The eigen-structures of real (skew) circulant matrices with some applications

THE EIGEN-STRUCTURES OF REAL (SKEW) CIRCULANT MATRICES WITH SOME APPLICATIONS ZHONGYUN LIU∗, SIHENG CHEN∗, WEIJIN XU∗,AND YULIN ZHANG† Abstract. The circulant matrices and skew-circulant matrices are two special classes of Toeplitz matrices and play vital roles in the computation of Toeplitz matrices. In this paper, we focus on real circulant and skewcirculant matrices. We first investigate their real Schur forms, which are closely related to the family of discrete cosine transform (DCT) and discrete sine transform (DST). Using those real Schur forms, we then develop some fast algorithms for computing real circulant, skew-circulant and Toeplitz matrix-real vector multiplications. Also, we develop a DCT-DST version of circulant and skew-circulant splitting (CSCS) iteration for real positive definite Toeplitz systems. Compared with the fast Fourier transform (FFT) version of CSCS iteration, the DCTDST version is more efficient and saves a half storage. Numerical experiments are presented to illustrate the effectiveness of our method. Key words. Real Schur form, real circulant matrices, real skew-circulant matrices, real Toeplitz matrices, CSCS iteration. AMS subject classifications. 15A23, 65F10, 65F15. 1. Introduction. Recall that a matrix T= (tjk)n−1 j,k=0 is said to be Toeplitz if tjk =tj−k; a matrix C= (cjk)n−1 j,k=0 is said to be circulant if cjk =cj−kand c−l=cn−lfor 1 ≤l≤n−1; and a matrix S= (sjk)n−1 j,k=0 is said to be skew-circulant if sjk =sj−kand s−l=−sn−lfor 1≤l≤n−1. Toeplitz matrices arise in a variety of applications in mathematics, scientific computing and engineering, for instance, signal processing, algebraic differential equation, time series and control theory, see e.g. [3] and a large literature therein. Those applications have motivated both mathematicians and engineers to develop specific algorithms for solving Toeplitz systems for instance [3, 6, 13] and references therein. The discrete Fourier transform (DFT) matrix F= (Fjk) is defined by Fjk =1 √nω−jk, j, k = 0,1,··· , n −1,where ω= exp(2π ni), i =√−1.(1.1) It is known that any circulant matrix Cand skew-circulant matrix Spossess the following Schur canonical forms [3, 10, 5], respectively, C=FΛF∗and S =˜ F˜ Λ˜ F∗,(1.2) ∗School of Mathematics and Statistics, Changsha University of Science and Technology, Changsha 410076, P. R. China ([email protected], [email protected], [email protected]). †Centro de Matem´atica, Universidade do Minho, 4710-057 Braga, Portugal ([email protected]). 1 2Z. Y. Liu, et al. where ˜ F=DF∗is an unitary matrix with D= diag1, eπ ni,··· , e(n−1)π ni, Λ and ˜ Λ are diagonal matrices, holding the eigenvalues of Cand Srespectively. Moreover, Λ and ˜ Λ can be obtained in O(nlog n) operations by using two FFTs of the first rows of Cand S, respectively. Due to the Schur canonical forms (1.2) of Cand S, the products Cxor Sxfor any vector xcan be computed by 3FFTs ( 1 FFT for computing eigenvalues) in O(nlog n) operations. Very often, circulant and skew circulant matrices are used to deal with Toeplitz issues. An important property is that a Toeplitz matrix Tcan be split into the following circulant and skew-circulant splitting (CSCS)[10] T=C+S(1.3) where C= (cjk) is a circulant matrix and S= (sjk) is a skew-circulant matrix, which are defined as follows. cjk =(1 2t0,if j=k, (tj−k+tj−k−n) 2,otherwise, and sjk =(1 2t0,if j=k, tj−k−tj−k−n 2,otherwise. Actually, due to Tx=Cx+Sx,Txcan be computed by 6 FFTs of n-vector. Also, any linear system of equations Cx=b(Sx=b) that contains circulant matrices (skew-circulant matrices) may be quickly solved by using the FFT. However, all operations, due to FFTs, are involved into complex arithmetics, even if C(S) and bare real. Now, one may ask when C and Sare real, could we find an analogue of (1.2) for Cand Sto avoid complex arithmetics in matrix-vector multiplication? and/or when C(S) and bare real, could we develop an algorithm which only involves real arithmetics for solving Cx=b(Sx=b)? This is the main motivation of this paper. The organization of this paper is as follows. In the next section, by exploring the eigenstructures of Cand S, we will give the real Schur forms of the circulant matrix Cand the skew-circulant matrix S. In Sections 3 and 4, with the real schur forms, we will develop a real method to fast calculate Toeplitz matrix-vector multiplication and an algorithm based on the CSCS iteration in [10] to solve Tx=bby real arithmetics. Numerical experiments are presented in Section 5 to show the effectiveness of our method. A brief conclusion and the acknowledgements are finally followed. 2. The Real Schur Forms of Real (Skew) Circulant Matrices. In the section, making use of the eigen-structures of a real circulant matrix Cand a real skew-circulant matrix S, we develop their corresponding real Schur forms. 2.1. Preliminaries. Let’s begin with some basic definitions. For convenience, throughout the paper, we define Jnthe permutation matrix of order nwith ones on the cross diagonal (bottom left to top right) and zeros elsewhere, and Ppq is the (q−p)×nrestriction matrix satisfying Ppq (xj)n−1 j=0 = (xj)q−1 j=p, q > p. Definition 2.1. [8] A vector x∈Rnis said to be symmetric if Jnx=xand skewsymmetric if Jnx=−x. The real eigen-structure 3 Now, let’s recall the definitions of DCTs and DSTs. The family of discrete trigonometric transforms consists of 8 versions of DCTs and corresponding 8 versions of DSTs [15, 12, 11]. In this paper, we only need four versions of them which will be used in the sequel. Definition 2.2. The DCT-I, DCT-II, DCT-V and DCT-VI matrices are defined as follows. CI n+1 =rn 2τjτkcos jkπ nn j,k=0 ,CV n=2 √2n−1τjτkcos 2jkπ 2n−1n−1 j,k=0 , CII n=rn 2τjcos j(2k+ 1)π 2nn−1 j,k=0 ,CVI n=2 √2n−1τjιkcos j(2k+ 1)π 2n−1n−1 j,k=0 , where τl(l=j, k)=(1,if l6= 0 and l6=n 1 √2,if l= 0 or l=n, and ιk=(1,if k6=n−1 1 √2,if k=n−1. Definition 2.3. The DST-I, DST-II, DST-V and DST-VI matrices are defined as follows. SI n−1=r2 nsin jkπ nn−1 j,k=1 ,SV n−1=2 √2n−1sin 2jkπ 2n−1n−1 j,k=1 , SII n=r2 nτjsin j(2k−1)π 2nn j,k=1 ,SVI n−1=2 √2n−1sin j(2k−1)π 2n−1n−1 j,k=1 , where τjis defined as in Definition 2.2. Note that all those transform matrices are all orthogonal. 2.2. Real Circulant Matrices. Let us first to investigate the real eigen-structure of a real circulant matrix C. It is shown in [5, 9] that the eigenvalues of a real circulant matrix can be arranged in the following order 1. λ= [λ0, λ1,··· , λm−1, λm,¯ λm−1,··· ,¯ λ1]T, where λ0, λm∈Rand n= 2m, 2. λ= [λ0, λ1,··· , λm,¯ λm,··· ,¯ λ1]T, where λ0∈Rand n= 2m+ 1. Partitioning F∗= [ f0,··· ,fn−1], we have fn−k=¯ fk. For any eigenvalue λk,Cfk= λkfkmeans C(fk+¯ fk) = λkfk+¯ λk¯ fkand C(fk−¯ fk) = λkfk−¯ λk¯ fk. If we denote λk=αk+iβkand fk=ˆ ck+iˆ sk, then we have C[ˆ ck,ˆ sk] = [ ˆ ck,ˆ sk]αkβk −βkαk, where αkand βk,ˆ ckand ˆ skare the real and pure imaginary parts of λkand fk, for k= 0,··· , n −1. Now, we show how to construct an orthogonal matrix Uwhich transforms Cinto its real Schur form. Notice that ˆ ck=ˆ cn−kand ˆ sk=−ˆ sn−k, for k= 0,··· , n −1, so we need only 4Z. Y. Liu, et al. normalize the first half of the vectors ˆ ckand ˆ skto get an orthogonal U. Namely, Ucan be chosen as follows, U=(ˆ c0,√2ˆ c1,··· ,√2ˆ cm−1,ˆ cm,√2ˆ sm−1,··· ,√2ˆ s1, n = 2m, ˆ c0,√2ˆ c1,··· ,√2ˆ cm,√2ˆ sm,··· ,√2ˆ s1, n = 2m+ 1.(2.1) A straightforward calculation shows that UTCU ≡Ω is real and has the following structure: Ω2m=                α0 α1β1 ... ... αm−1βm−1 αm −βm−1αm−1 ... ... −β1α1                (2.2) or Ω2m+1 =              α0 α1β1 ... ... αmβm −βmαm ... ... −β1α1              ,(2.3) which can be transformed into the real Schur canonical form by a permutation. Therefore we also refer to (2.2) or (2.3) as the real Schur form of C. This leads to the following theorem. Theorem 2.4 (Real Schur form of real circulant matrices). Let Ube defined as in (2.1). If the circulant matrix Cis real, then UTCU = Ω is the real Schur form of C. Proof. The proof can be directly given by the above analysis and thus omitted. 2.3. Real Skew-Circulant Matrices. Similarly, the eigenvalues of a real skew-circulant matrix Scan be arranged in the following order 1. ˜ λ= [˜ λ0,··· ,˜ λm−1,¯ ˜ λm−1,··· ,¯ ˜ λ0]T, where n= 2m, 2. ˜ λ= [˜ λ0,··· ,˜ λm−1,˜ λm,¯ ˜ λm−1,··· ,¯ ˜ λ0]T, where ˜ λm∈Rand n= 2m+ 1. The next procedure is very much like section 2.2. Let ˜ F∗= [ ˜ f0,··· ,˜ fn−1]. Analogously, denoting ˜ λk= ˜αk+i˜ βkand ˜ fk=˜ ck+i˜ sk, then we have S[˜ ck,˜ sk]=[˜ ck,˜ sk]˜αk˜ βk −˜ βk˜αk. Due to ˜ ck=˜ cn−k−1and ˜ sk=−˜ sn−k−1, for k= 0,··· , n−1, we can choose an orthogonal The real eigen-structure 5 matrix ˜ Uas follows, ˜ U=(√2˜ c0,··· ,√2˜ cm−1,√2˜ sm−1,··· ,√2˜ s0,if n= 2m, √2˜ c0,··· ,√2˜ cm−1,˜ cm,√2˜ sm−1,··· ,√2˜ s0,if n= 2m+ 1,(2.4) which makes ˜ UTS˜ U≡Σ being of the following structure. Σ2m=            ˜α0˜ β0 ... ... ˜αm−1˜ βm−1 −˜ βm−1˜αm−1 ... ... −˜ β0˜α0            (2.5) or Σ2m+1 =              ˜α0˜ β0 ... ... ˜αm−1˜ βm−1 ˜αm −˜ βm−1˜αm−1 ... ... −˜ β0˜α0              .(2.6) Similarly, we refer to (2.5) or (2.6) as the real Schur form of S. Based on the above analysis, we conclude the following theorem. Theorem 2.5 (Real Schur form of real skew-circulant matrices). Let ˜ Ube defined as in (2.4). If the skew-circulant matrix Sis real, then ˜ UTS˜ U= Σ is the real Schur form of S. We remark here that different from C=FΛF∗(i.e., Cis factorized into the product of 3 complex matrices), Theorems 2.4 - 2.5 tell us that both Cand Scan be factorized into the products of 3 real matrices, respectively. This fact allows us to fast calculate matrix-vector multiplication and solve Cx=band Sx=bby only real operations. In the next two sections, we will derive this strategy. 3. Fast Matrix-vector Multiplication. In this section, we first reduce the matrices Uand ˜ Uinto simpler forms by exploiting the structures of Uand ˜ U, then show how to fast compute Ω in (2.2) or (2.3) and Σ in (2.5) or (2.6), and finally develop fast algorithms for computing Cx,Sxand Tx. Recall Definition 2.1, P1nˆ ck(remove the first entry of ˆ ck) is a symmetric vector and P1nˆ sk (remove the first entry of ˆ sk) is a skew-symmetric vector. Notice that if we delete the first row of U, then the first m+ 1 columns of the submatrix are all symmetric, and the last columns are all skew-symmetric. Therefore, the Uof (2.1) can be partitioned into the following form 6Z. Y. Liu, et al. U2m=         σ1qT m+1 0 ˆ C−SI m−1Jm−1 σ1vT m+1 0 Jm−1ˆ CJm−1SI m−1Jm−1         and U2m+1 =   σ1pT m+1 0 ˆ C−SV mJm Jmˆ CJmSV mJm    where σ1=q2 n,σ2=1 √2,pm+1 = ( 1 √2,1,··· ,1)T,qm+1 = ( 1 √2,1,··· ,1,1 √2)T,vm+1 = (1 √2,−1,··· ,(−1)m−1,(−1)m √2)T, and ˆ C=(σ2P1,mCI m+1 ∈R(m−1)×(m+1),if n= 2m, σ2P1,m+1CV m+1 ∈Rm×(m+1),if n= 2m+ 1. Similarly, the ˜ Uof (2.4) can be partitioned into the following form ˜ U2m=       σ1eT0 ˜ C S II mJm 0σ1uT mJm −Jm−1˜ CJm−1SII mJm        and ˜ U2m+1 =    σ1pT m+1 0 ˜ C S VI mJm −Jm˜ CJmSVI mJm     where um= (1,−1,··· ,(−1)m−1)Tand ˜ C=(σ2P1mCII m∈R(m−1)×m,if n= 2m, σ2P1mCVI m+1 ∈Rm×(m+1),if n= 2m+ 1. Now we construct an orthogonal matrix Qof the form Q=                        1 √2    √2 Im−1Jm−1 √2 −Jm−1Im−1     ,if n= 2m, 1 √2  √2 ImJm −JmIm ,if n= 2m+ 1. (3.1) Then we can get the following simple formulas. Theorem 3.1. Let U,˜ Uand Qbe defined as in (2.1),(2.4) and (3.1), respectively. Then we have QU =           CI m+1 JmSI m−1Jm,if n= 2m, CV m+1 JmSV mJm,if n= 2m+ 1. (3.2) The real eigen-structure 7 and QT˜ U=           CII m JmSII mJm,if n= 2m. CVI m+1 JmSVI mJm,if n= 2m+ 1. (3.3) Proof. The proof of this theorem can be completed by a tedious straightforward calculation and thus omitted. The Theorem 3.1 provides us a fast way to the calculation of Ux, it can be obtained by 1 DCT-I of (m+ 1)-vector and 1 DST-I of (m−1)-vector if n= 2m, and 1 DCT-V of (m+ 1)- vector and 1 DST-V of m-vector if n= 2m+ 1. The product ˜ Uxcan be obtained by a similar mode by employing the second and sixth versions of DCT and DST. The entries of Ω and Σ can be computed by ΩUTe1= (QU)TQCe1and Σ ˜ UTe1= (QT˜ U)TQTSe1,(3.4) where e1= (1,0,··· ,0)T. The left-hand sides of (3.4) are as follows, respectively, √nΩUTe1=   (α0 √2, α1,··· , αm−1,αm √2,−βm−1,··· ,−β1)T, n = 2m, (α0 √2, α1,··· , αm,−βm,··· ,−β1)T, n = 2m+ 1,(3.5) and √nΣ˜ UTe1=   (˜α0,··· ,˜αm−1,−˜ βm−1,··· ,−˜ β0)T, n = 2m, (˜α0,··· ,˜αm−1,˜αm √2,−˜ βm−1,··· ,−˜ β0)T, n = 2m+ 1.(3.6) This means we only need 1 DCT and 1 DST of about n 2-vector to get Ω or Σ. Now, we show the calculations of Cxand Sxfor any real n-vector xusing DCT and DST. According to Theorem 2.4 and Theorem 3.1, Cxcan be easily obtained by three DSTs and three DCTs (version I or V) of about n 2-vector. As for the storage required, we need one temporary n-vector and an extra n-vector for storing Ω. In fact, we don’t need to compute and store Uand Ω explicitly. It can be written as the following Algorithm 1. If we compute Cxby (1.2), it requires three FFTs of n-vector, and one temporary complex n-vector and an extra complex n-vector for storing Λ, equivalently, two temporary real n-vectors and two extra real n-vectors for storing Λ. Similarly, according to Theorem 2.5 and Theorem 3.1, we develop the following Algorithm 2 for computing the product Sx, which can be obtained by three DSTs and three DCTs (version II or VI) of about n 2-vector. From (1.3), we have that a Toeplitz matrix-vector multiplication Tx=Cx+Sxcan be fast calculated by employing the Algorithms 1 - 2. 8Z. Y. Liu, et al. Algorithm 1 To calculate Cx 1: Compute v=Qc1directly. 2: Compute ˆ v= (QU)Tvby DCT and DST. 3: Form Ω. 4: Compute y1=Qxdirectly. 5: Compute y2= (QU)Ty1by DCT and DST. 6: Compute y3= Ωy2directly. 7: Compute y4= (QU)y3by DCT and DST. 8: Compute QTy4, i.e., Cx. Algorithm 2 To calculate Sx 1: Compute u=QTs1directly. 2: Compute ˆ u= (QT˜ U)Tuby DCT and DST. 3: Form Σ. 4: Compute z1=QTxdirectly. 5: Compute z2= (QT˜ U)Tz1by DCT and DST. 6: Compute z3= Σz2directly. 7: Compute z4= (QT˜ U)y3by DCT and DST. 8: Compute Qz4, i.e., Sx. 4. Solving Tx=bby the CSCS iteration. Consider the iterative solution to a large scale system of linear equations Tx=b,(4.1) where T∈Rn×nis a Toeplitz matrix and b∈Rn. Based on the splitting (1.3), Ng proposed in [10] the following CSCS iteration for solving (4.1). The CSCS iteration: Given an initial guess x(0), for k= 0,1,··· ,until {x(k)}converges, compute    (θI +C)x(k+1 2)= (θI −S)x(k)+b, (θI +S)x(k+1) = (θI −C)x(k+1 2)+b, (4.2) where θis a given positive constant. It is shown in [10] that the CSCS iteration converges unconditionally, if both Cand Sare positive definite. Then applying (1.2) to (4.2), we get The FFT version of CSCS iteration: Given an initial guess x(0), for k= 0,1,···, The real eigen-structure 9 until {x(k)}converges, compute    F(θI + Λ)F∗x(k+1 2)=˜ F(θI −˜ Λ) ˜ F∗x(k)+b, ˜ F(θI +˜ Λ) ˜ F∗x(k+1) =F(θI −Λ)F∗x(k+1 2)+b. (4.3) where θis a given positive constant. In the preparatory stage, two FFTs of n-vector for computing Λ and ˜ Λ are required. In the iterative stage, six FFTs of n-vector for solving (4.3) are needed. Therefore the computational complexity is O(nlog n) complex flops at each iteration. However, all operations, due to FFTs, are involved into complex arithmetics, even if Tand bare real. In this section, we develop the DCT-DST version of (4.2) based on the DCT and DST. Also, we compare the computational cost of our version with the FFT version of the CSCS iteration (4.3). Note by Theorem 2.4 and Theorem 2.5 that the CSCS iteration (4.2) can be reformulated as the following form.    U(θI + Ω)UTx(k+1 2)=˜ U(θI −Σ) ˜ UTx(k)+b, ˜ U(θI + Σ) ˜ UTx(k+1) =U(θI −Ω)UTx(k+1 2)+b. (4.4) The equation (4.4) can be further reduced into a simpler form due to Theorem 3.1. For example, consider the case n= 2m(The odd case is similar to the even case), we have the following version, The DCT-DST version of CSCS iteration: Given an initial guess x(0) ∈Rn, compute x(k), for k= 0,1,···, until {x(k)}converges:                                QTCI m+1 SI m−1(θI + Ω) CI m+1 SI m−1(Qx(k+1 2)) =QCII m SII m(θI −Σ) CII m SII mT (QTx(k)) + b, QCII m SII m(θI + Σ) CII m SII mT (QTx(k+1)) =QTCI m+1 SI m−1(θI −Ω) CI m+1 SI m−1(Qx(k+1 2)) + b, (4.5) where θis a given constant. We emphasize here that our version has the same convergence rate and optimal parameter as the CSCS iteration does. The computational complexity. The iteration (4.5) consists of the preparatory stage and computational stage. In preparatory stage, we only need to calculate the matrices Ω and Σ which can be obtained by two DCTs and DSTs of about n/2-vector, see (3.4), (3.5) and (3.6).