# Real Fast Structure-Preserving Algorithm for Eigenproblem of Complex Hermitian Matrices

Authors: Jiangzhou Lai; Linzhang Lu
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

It is well known that the flops for complex operations are usually 4 times of real cases. In the paper, using real operations instead of complex, a real fast structure-preserving algorithm for eigenproblem of complex Hermitian matrices is given. We make use of the
real symmetric and skew-Hamiltonian structure transformed by Wilkinson's way, focus on symplectic orthogonal similarity transformations and their structure-preserving property, and then reduce it into a two-by-two block tridiagonal symmetric matrix. Finally a real algorithm can be quickly obtained for eigenvalue problems of the original Hermitian matrix. Numerical experiments show that the fast algorithm can solve real complex Hermitian matrix efficiently, stably, and with high precision.


---

## Body

## 1. Introduction

We consider eigenvalue problems of the complex Hermitian matrix of the form(1)Hx=λx,H∈ℂn×n,x∈ℂn,
where λ is a real scalar and H∈ℂn×n is a complex Hermitian matrix of the form
(2)H=A+iB,
where A is a real symmetric matrix; that is, AT=A∈ℝn and B is a real skew-symmetric nonzero matrix, that is, B=-BT≠0∈ℝn, given by
(3)A=HT+H2,B=H-HT2i,
and x∈ℂn can be written as x=u+iv, where u, v∈ℝn.We know that the property of the complex Hermitian matrix is also possessed by the symmetry matrix; therefore little work has been done on complex Hermitian matrix. However, from the view of numerical computation, the biggest difference between these two types of matrices is a complex Hermitian matrix that involves complex operations, while the latter does not. It is well known that a complex Hermitian matrix is unitary similar to a real diagonal matrix; that is, it has real eigenvalue but complex eigenvector. Because the flops for the complex operation are usually 4 times of the real case, so how to reduce the computation cost of the eigenvalue problem of complex Hermitian matrices is very meaningful and worthy of research.In 1965, the famous numerical analysis expert Wilkinson [1] gave a way to transform a complex Hermitian matrix into real case. To be more precise, if (λ,u+iv) is the eigenpair of H, then
(4)Hx=(A+iB)(u+iv)=λ(u+iv).Since the eigenvalues ofH are real, the previous equation can be written as
(5)[A-BBA][uv]=λ[uv].So the eigenproblem (1) is transformed into the eigenproblem of the real symmetry matrix
(6)S=[A-BBA]∈ℝ2n×2n.Note that the matrixS in (6) is not only real symmetry but also skew-Hamiltonian. Therefore, in the following discussion, inspired by the algorithms for the eigenproblem of Hamiltonian (skew-Hamiltonian) matrix [2–11], we devise an algorithm by taking full advantage of the special structure of (6).In this paper we first prove that orthogonal symplectic similarity transformations not only preserve the symmetry and skew-Hamiltonian structure of (6), but also eigenvectors of (6), since the eigenvectors of (6) have special structure; that is, it is also orthogonal symplectic. Therefore, by these transformations, a real stable, accurate, and fast method is devised to reduce (6) into a block symmetry tridiagonal matrix, and then we get the eigenproblem of the original Hermitian matrix. When implementing the method, although its dimension doubles the one of H, we can make full use of the special structure of S to avoid superfluous calculation and storage.The remainder of the paper is organized as follows. In Section2, we show the basic orthogonal symplectic matrix and its algorithm. In Section 3, we present the overall real algorithm for eigenproblem of complex Hermitian matrix. Some implementation aspects are discussed in Section 4. Finally, in Section 5, the results of numerical experiments are discussed to end the paper.


## 2. Basic Orthogonal Symplectic Matrix and Its Algorithm

It is well known that by choosing a proper unitary vectorw, Householder transformation can zero some entries of a given vector; likewise, we defined Householder symplectic transformation which has the following form.Definition 1.
A real2n×2n matrix S is a symplectic matrix if SJST=J, where ST is the transpose of S  and J is the orthogonal matrix:
(7)J=[0In-In0]∈ℝ2n×2n.
Here In is the identity n×n matrix and 0 is the zero n×n matrix.Definition 2.
Lettingk∈{1,2,…,n}, a matrix H(k,w)∈ℝ2n×2n is called Householder symplectic matrix, if it has the form
(8)H(k,w)=[diag(Ik-1,P)00diag(Ik-1,P)]∈ℝ2n×2n,
where
(9)P=I-2wwTwTw,w∈ℝn-k+1.
Obviously, H=diag(Ik-1,P) is Householder transformation of order n, when w=0, H(k,w)=I2n.Note that Householder symplectic matrixH(k,w) is just a direct sum of two “ordinary” n-by-n Householder matrices, specially when w=0, H(k,w)=I2n. So based on the algorithm of “ordinary” Householder transformation, we obtain the following algorithm for Householder symplectic transformation H(k,w) of order 2n.Algorithm 3 (Householder symplectic transformation [8]).
Givenk(1≤k<n) and y,z∈ℝn, the following algorithm determines w=(wk,…,wn)T, such that if
(10)H(k,w)(yz)=(vx),
then xi=0,i=k+1,…,n:
(11)σ:=(zk2+zk+12+⋯+zn2)1/2,wk:=zk+sign(zk)σ,wi:=zi
end.H
(
k
,
w
) is completely determined by kth to nth components of vector z; it makes x satisfy xi=0 for i=k+1,…,n. Note that by interchanging the roles of y and z, the previous algorithm can be used to determine H(k,w) such that vi=0 for i=k+1,…,n.Householder symplectic transformation can be used to zero large portions of a vector as Householder transformation; likewise, we can define Givens symplectic transformation, which can be used to zero single entries.Definition 4.
Givenk∈{1,2,…,n}, a matrix J(k,θ)∈ℝ2n×2n is called Givens symplectic matrix, if it has the form
(12)J(k,θ)=[CS-SC]=[Ik-10k-1cosθsinθIn-k0n-k0k-1Ik-1-sinθcosθ0n-kIn-k],
where
(13)C=diag(Ik-1,cosθ,In-k),S=diag(0k-1,sinθ,0n-k).
J
(
k
,
θ
) is an “ordinary” 2n-by-2n Givens rotation that rotates in planes k and n+k. The corresponding algorithm is as follows.Algorithm 5 (Givens symplectic transformation [8]).
Givenk(1≤k<n) and y,z∈ℝn, the following algorithm determines c=cosθ, s=sinθ such that if(14)J(k,w)(yz)=(vx),
then xk=0:
(15)σ≔(yk2+zk2)1/2,ifσ=0thenc:=1ands:=0elsec:=yk/σands:=zk/σend.


## 3. The Reduction of the Symmetry and Skew-Hamiltonian MatrixS

Now we consider some properties of the symmetry and skew-Hamiltonian matrixS of (6).Lemma 6.
If(λ,[uv]) is the eigenpair of
(16)S=[A-BBA],
then (λ,[-vu]) is also its eigenpair.Proof.
If(17)[A-BBA][uv]=λ[uv],
then we have
(18)[A-BBA][-vu]=λ[-vu].According to Lemma6, there are two different eigenvectors corresponding to the same eigenvalue, so it can be easy to prove that
(19)[A-BBA][U-VVU]=Λ[U-VVU],
where Λ={λ1,…,λn} is the set of eigenvalues of (6) and
(20)[U-VVU]∈ℝ2n×2n
is the eigenvectors of (6). Note that the matrix in (20) is symplectic orthogonal. The following lemma tells us that the product of symplectic orthogonal matrices is also symplectic orthogonal.Lemma 7.
Supposing that2n-by-2n real matrices of the form
(21)Q=[Q1Q2-Q2Q1],Q~=[Q~1Q~2-Q~2Q~1]
are both orthogonal and symplectic, then their product QQ~ is both orthogonal and symplectic.Lemma7 shows that the special structure of eigenvectors of S can be preserved by the symplectic orthogonal similarity transformations. The following theorem is the foundation of the fast and stable structure-preserving algorithm given in the following discussion.Theorem 8.
SupposeAT=A∈Rn×n,BT=-B≠0∈Rn×n, and
(22)Q=[Q1Q2-Q2Q1]
is an orthogonal symplectic matrix. If
(23)[A~11A~12A~21A~22]
satisfy
(24)[A~11A~12A~21A~22]=[Q1Q2-Q2Q1]T[A-BBA][Q1Q2-Q2Q1],
then
(25)A~11=A~22,A~21=-A~12,A~11T=A~11,A~12T=-A~12.Proof.
It is easy to obtain(26)A~11=Q1TAQ1+Q1TBQ2-Q2TBQ1+Q2TAQ2,A~22=Q2TAQ2-Q2TBQ1+Q1TAQ1+Q1TBQ2,A~12=Q1TAQ2-Q1TBQ1-Q2BQ2-Q2TAQ1,A~21=Q2TAQ1+Q2TBQ2+Q1TBQ1-Q1TAQ2.
Then
(27)A~11=A~22,A~21=-A~12.
MoreoverAT=A,BT=-B, hence we have
(28)A~11T=A~11,A~12T=-A~12.The previous theorem shows that the symplectic similarity transformation preserves the symmetry and skew-Hamiltonian structure of(29)S=[A-BBA],(AT=A,BT=-B).This structure-preserving property provides a way to transform the matrix into a special skew-Hamiltonian matrix:(30)[T00T]∈ℝ2n×2n,
where T∈ℝn×n is a real symmetry tridiagonal matrix. So we have the following results, which provide a theoretical basis for the real fast structure-preserving algorithm for eigenproblems of a Hermitian matrix.Theorem 9.
LettingS=[A-BBA]∈ℝ2n×2n, where AT=A,BT=-B, then there exists an orthogonal symplectic matrix Q=[Q1Q2-Q2Q1]∈ℝ2n×2n such that
(31)QTSQ=[T00T]∈ℝ2n×2n,
where T is a symmetric tridiagonal matrix. Suppose that T has the eigenvalue decomposition
(32)Q~1TTQ~1=Λ,
that is
(33)[Q~100Q~1]T[T00T][Q~100Q~1]=[Λ00Λ].
Then eigenvalues of S are
(34)λ(S)=λ(T)⊔λ(T),
and eigenvectors are the columns of
(35)QQ~=[Q1Q2-Q2Q1][Q~100Q~1]=[Q1Q~1Q2Q~1-Q2Q~1Q1Q~1].We now consider how an orthogonal symplecticQ∈ℝ2n×2n can be determined to reduce S=[A-BBA] (AT=A,BT=-B) to
(36)QTSQ=[T00T],
where T is a symmetry tridiagonal matrix. It is worth illustrating a few steps of the algorithm before presenting it in full generality. We depict an 8-by-8 symmetry and skew-Hamiltonian matrix (n=4) as follows:
(37)S=[A-BBA]=[xxxx0xxxxxxxx0xxxxxxxx0xxxxxxxx00xxxxxxxx0xxxxxxxx0xxxxxxxx0xxxx]S=[A-BBA]xxxxxxxX(n=4).
The zero diagonals in B follow from skew symmetry. It is easy to verify that the (2, 1) and (1, 2) block of this matrix can be zeroed by applying a sequence of Householder and Givens symplectic similarity transformations.The first step is to zerob31 and b41 using a Householder symplectic H1=H(2,w). This matrix can be determined by executing Algorithm 3 with k=2, y=Ae1, and z=Be1. Note that after performing, only the rows and columns 2, 3, and 4 of A, B, and -B are affected, the result being
(38)S=[A-BBA]:=H1SH1T=[xxxx0x00xxxxx0xxxxxx0x0xxxxx0xx00x00xxxxx0xxxxxx0x0xxxxx0xx0xxxx]=[A-BBA]:=H1Sxx(n=4).
Then (1,3) and (1,4) positions of B are zeroed; meanwhile the same positions of -B are zeroed because Householder symplectic similarity transformation preserves the structure of S.The next step is to zerob21 using Givens similarity J1=J(2,θ), which can be achieved by applying Algorithm 5 with k=2, y=Ae1, and z=Be1. And notice that only the second row and column of A, B and -B be affected. This implies that
(39)S=[A-BBA]:=J1SJ1T=[xxxx0000xxxx00xxxxxx0x0xxxxx0xx00000xxxx00xxxxxx0x0xxxxx0xx0xxxx]S=[A-BBA]:==(n=4).Next we compute a Householder symplecticH(2,w) to zero a31 and a41, which can be achieved by applying Algorithm 3 with k=2, y=Be1, and z=Ae1. Denote G1=H(2,w), the result being
(40)S=[A-BBA]:=G1SG1T=[xx000000xxxx00xx0xxx0x0x0xxx0xx00000xx0000xxxxxx0x0x0xxx0xx00xxx]S=[A-BBA]=xxx(n=4).Note that the(1,3) and (1,4) positions of A are zeroed. This completes the zeroing in the first column of S. Likewise we can zero the 2th of S with the property that
(41)S=[A-BBA]=[xx000000xxx000000xxx000x00xx00x00000xx000000xxx0000x0xxx00x000xx]S=[A-BBA]=xxx(n=4).Notice that only the(3,4) and (4,3) positions of B are nonzero; they are zeroed only by applying Givens symplectic similarity transformation G3=J(4,θ) to S and we get
(42)S=G3SG3T=[T00T]=[xx000000xxx000000xxx000000xx00000000xx000000xxx000000xxx000000xx]S=[T00T]=mm(n=4).This completes the zeros ofS with n=4. The similar method can be applied to zero the general matrix of order 2n. Overall we have the following procedure “Reduction to block symmetry tridiagonal matrix.” Given matrix
(43)S=[A-BBA]withAT=A,BT=-B.The following algorithm updatesA to the symmetry tridiagonal matrix T; meanwhile B and -B are reduced into an n×n zero matrix; moreover we have λ(S)=λ(T)⊔λ(T).Algorithm 10 (reduction to block symmetry tridiagonal matrix).
Fork=1,…,n-2 do:
if k≤n-2 do:
(1)  Apply Algorithm3 with y=Aek and z=Bek to determine Hk=H(k+1,w). Update is as follows:
(44)[A-BBA]:=Hk[A-BBA]HkT.
(2)  Apply Algorithm5 with y=Aek and z=Bek to determine Jk=J(k+1,θ).  Updates are as follows:
(45)[A-BBA]:=Jk[A-BBA]JkT.
(3)  Apply Algorithm3 with y=Bek,z=Aek to determine Gk=H(k+1,w). Updates are as follows:
(46)[A-BBA]:=Gk[A-BBA]GkT.
Ifk=n-1 do the folowing.
(4) Apply Algorithm5 with y=Aen-1 and z=Ben-1 to determine Jn-1=J(n,θ). Update is as follows:
(47)[A-BBA]:=Jn-1[A-BBA]Jn-1T.After completion, we get the block symmetry tridiagonal matrix[T00T], where T is a symmetry tridiagonal matrix. The orthogonal matrix QT is given by
(48)QT=Jn-1(Gn-2Jn-2Hn-2)⋯(G2J2H2)(G1J1H1).The method of reductionS into block symmetry tridiagonal matrix is given by Algorithm 10; in the former n-2 steps, each step we premultiply by one Givens symplectic similarity and two Householder symplectic similarity and by postmultiply their transpose. Finally in the n-1 step, only the Givens symplectic similarity is done. According to (48), it is clear that QT is the product of 3(n-2)+1 orthogonal matrices, which is backward stability.Algorithm 11 (the overall algorithm).
(1) Transform Hermitian matrixH into the symmetry and skew-Hamiltonian matrix:
(49)S=[A-BBA]whereAT=A,BT=-B.
(2)  Apply Algorithm10 to reduce S into a block symmetry tridiagonal matrix:
(50)QTHQ=[T00T].
(3)  Compute eigenpairs (Λ, X) of T via QR.
(4)  Setting(51)[Q11Q21]=[Q1Q2-Q2Q1][X0]=[Q1X-Q2X],
then Λ is eigenvalue of H, and Q11+iQ21 is the eigenvector of H.


## 4. Implementation Aspects

When implementing Algorithm10, the matrix S must be taken full advantage of in order to avoid superfluous calculation and storage. In this section there are two parts; the first is to consider the computational cost and storage for the reduction to this block symmetry tridiagonal matrix, and the second is the computational cost and storage of the eigenvector QT.Due to the special structureA=AT and BT=-B, shared by A and B, only n2 workspace is required for matrix
(52)S=[A-BBA]whereAT=A,BT=-B.Now we consider the computation cost for reducingS to the block symmetry tridiagonal form, for Householder symplectic similarity, let
(53)Hk=H(k+1,w)=diag(Ik,P,Ik,P),(54)S=[A11A21B11B21A21TA22-B21TB22B11T-B21TA11A12TB21TB22TA12A22]kn-kkn-kkn-kkn-k,
then
(55)HkSHk=[A11PA12TB11PB21A12PPA22P-B21TPPB22PB11T-PB21TA11PA12TB21TPPB22TA12PA22P]kn-kkn-kkn-kkn-k.Inkth step, updating S to HkSHk is equal to two ordinary Householder similarities for n×n matrix, one for symmetry matrix A and one for skew-symmetry matrix B, so we only need to compute PA12T, PA22P, PB21, and PB22P. Likewise, we can update S to GkSGk by the same way.Algorithm 12 (Householder symplectic similarity).
For
k
=
1
:
n
-
2
(56)
[
v
,
β
]
=
house
(
B
(
k
+
1
:
n
,
k
)
)
,
w
=
β
B
(
k
+
1
:
n
,
k
:
n
)
T
v
,
B
(
k
+
1
:
n
,
k
:
n
)
=
B
(
k
+
1
:
n
,
k
:
n
)
-
v
w
T
,
w
=
β
B
(
k
:
n
,
k
+
1
:
n
)
v
T
,
B
(
k
:
n
,
k
+
1
:
n
)
=
B
(
k
:
n
,
k
+
1
:
n
)
-
w
v
T
,
w
=
β
A
(
k
+
1
:
n
,
k
:
n
)
T
v
,
A
(
p
+
1
:
n
,
p
:
n
)
=
A
(
k
+
1
:
n
,
k
:
n
)
-
v
w
T
,
w
=
β
A
(
k
:
n
,
k
+
1
:
n
)
v
,
A
(
k
:
n
,
k
+
1
:
n
)
=
A
(
k
:
n
,
k
+
1
:
n
)
-
w
v
T
end.Updating by Givens symplectic similarity is equally simple; settingc=cos(θ), s=sin(θ)  then (57)Jk=J(k+1,θ)=[Ik0kcsIn-k-10n-k-10kIk-sc0n-k-1In-k-1],S=[A11A21A31B11B21B31A12A22A32B12B22B32A13A23A33B13B23B33-B11-B21-B31A11A21A31-B12-B22-B32A12A22A32-B13-B23-B33A13A23A33]k+11n-k-2k+11n-k-2k+11n-k-2k1n-k-2.DenotingSJ=(JkSJkT), then we have
(58)SJ(1:2n,1:n)=[A11cA21+sB21A31B11cB21-sA21B31cA12-sB12A22cA32-sB32cB12+sA12B22cB32+sA32A31TA32TA33B13cB23-SA23B33]k+11n-k-2k+11n-k-2k+11n-k-2.Note that updatingS into JkSJkT such that B(k+1,k)=0 affects the rows  (resp. the columns) k+1, n+k+1 of S, that is, the rows(columns) k+1, n+k+1 of A and B. According to (58), we only need to compute cA21+sB21, cA32-sB32, cB21-sA21 and cB32+sA32 because SJ is also symmetry and skew-Hamiltonian matrix. Therefore, only 2(n-k) flops are required in kth step, altogether ∑k=1n-12(n-k)=n2-n flops. Bellow is the algorithm for JkSJkT in kth step.Algorithm 13 (Givens symplectic similarity).
Consider(59)a=A(k+1,k),b=B(k+1,k),[c,s]=givens(a,b),A(k+1,k)=cA(k+1,k)+sB(k+1,k),A(k,k+1)=A(k+1,k),fori=k+2:nS(i,k+1)=cA(i,k+1)-sB(i,k+1),B(i,k+1)=cB(i,k+1)+sA(i,k+1),B(k+1,i)=-B(i,k+1),A(i,k+1)=S(i,k+1),A(k+1,i)=A(i,k+1)endB(k+1,k)=0,B(k,k+1)=0.Overall(4/3)n3×4+𝒪(n2)=(16/3)n3+𝒪(n2) flops are required for reducing S to [T00T].This reduction includes 4 times computational cost of Householder tridiagonalized forn×n symmetry matrix and n2-n flops for Givens similarity transformation. Since one complex operation is four times of real operation, (16/3)n3+𝒪(n2) flops are required for complex Hermitian tridiagonalization.Now we consider the computation and storage problem for eigenvectors(60)QT=Jn-1(Gn-2Jn-2Hn-2)⋯(G2J2H2)(G1J1H1).To obtainQT, we can premultiply only by n×n matrices by using the special structure of Jk, Gk, and Hk though they are 2n-by-2n matrices. When premultiplying by Householder symplectics Hk, let a 2n×2n matrix QT be of the form
(61)QT=[QN×N-UN×NUN×NQN×N]=[Q11Q21U11U21Q12Q22U12U22-U11-U21Q11Q21-U12-U22Q12Q22]kn-kkn-kkn-kkn-k.So twon2 workspace is required for QT, one for QN×N and one for UN×N. Premultiplying by Hk=diag(Ik+1,PA,Ik+1,PA) in kth step, we have
(62)HkQT=[Q11PAQ21U11PAU21Q12PAQ22U12PAU22-U11-PAU21Q11PAQ21-U12-PAU22Q12PAQ22]kn-kkn-kkn-kkn-k.
Therefore, only Q21, Q22 and the real part of QT are updated to PAQ21, PAQ22 and also imaginary part U21, U22 updated to PAU21, PAU22. The main work is computing the product of Householder matrix H=I-βvvT∈ℝ(n-k)×(n-k) and [Q21Q22] ∈ℝ(n-k)×(n), and the product of H and [U21U22] ∈ℝ(n-k)×(n), where the first columns of Q21 and U21 are zero. So 4(n-k)(n-k-1) flops are required each step, altogether ∑k=1n-28(n-k)(n-k-1)=(8/3)n3+𝒪(n2) flops. Likewise premultiplying by Gk requires other (8/3)n3+𝒪(n2) flops.Algorithm 14 (computation ofHkQT).
For k=1:n-2(63)[v,β]=house(B(k+1:n,k)),w=βQ1(k+1:n,2:n)Tv,Q1(k+1:n,2:n)=Q1(k+1:n,2:n)-vwT,w=βQ2(k+1:n,2:n)Tv,Q2(k+1:n,2:n)=Q2(k+1:n,2:n)-vwTend.For Givens symplectic matrixJk, set(64)QT=[Q11Q21Q31U11U21U31Q12Q22Q32U12U22U32Q13Q23Q33U13U23U33-U11-U21-U31Q11Q21Q31-U12-U22-U32Q12Q22Q32-U13-U23-U33Q13Q23Q33]k1n-k-1k1n-k-1k1n-k-1k1n-k-1.Premultiplying byJk, denoting QJT=JkQT, we have
(65)QJT(1:2n,1:n)=[Q11cQ21+sU21Q310cU21-sQ210Q12cQ22+sU22Q320cU22-sQ220Q13cQ23+sU23Q330cU23-sQ230]k1n-k-1k1n-k-1k1n-k-1.
Therefore, each step updating the k+1th and n+k+1th row of QT requires 4(n-1) flops, altogether 4(n-1)2 flops. Below is the algorithm.Algorithm 15 (computation ofJkQT).
Fork=1:n-2(66)T1=Q1,T2=Q2,Q1(k+1,2:n)=cQ1(k+1,2:n),Q2(k+1,2:n)=sQ2(k+1,2:n),Q2(k+2:n,2:n)=0,Q1(k+1:n,2:n)=Q1(k+1:n,2:n)+Q2(k+1:n,2:n),T1(k+1,2:n)=-sT1(k+1,2:n),T1(k+2:n,2:n)=0,T2(k+1,2:n)=cT2(k+1,2:n),Q2(k+1:n,2:n)=T1(k+1:n,2:n)+T2(k+1:n,2:n)
end.
Overall in each step we only need to premultiply twon×n Householder matrix and one n×n Givens matrix in the procedure of computing the eigenvector.


## 5. Numerical Experiments

All codes were run in MATLAB 7.04 in double precision and all experiments were performed on a personal computer with 1.41 Ghz central processing unit (AMD Sempron(tm) Processor 2500+), 520 MB memory, and Windows XP system. Machine epsiloneps≈2.22×10-16.Example 16.
To illustrate Algorithm11 more clearly, we give a simple example and show the processes and results of Algorithm 11 of H1:
(67)H1=[1-3i-4i3i14i4i4i1].
We first transformH1 into symmetry and skew-Hamiltonian matrix S:
(68)H→S=[100034010-304001-4-400-3-410030-4010440001].
Then
(69)S→HTSH=[100050010-50-40010400-501005040100-40001]→JTHTSHJ=[150000514000041000000150000514000041]=[T00T],
where
(70)H=[10000000.60.800000.8-0.600000010000000.60.800000.8-0.6],J=[1000000000-10001000000100010000000001].
Denoting that(D,V) is the eigenpair of T(71)V=[0.5522-0.62470.5522-0.70710.00000.70710.44170.78090.4417],D=[-5.40310001.00000007.4031],
then
(72)Q[V0]=H1J[V0]=[0.5522-0.62470.55220.35340.62470.3534-0.2650-0.4685-0.2650000-0.42430.00000.4243-0.56570.00000.5657]
is the first half of eigenvectors of S. The computing eigenproblem of H1 is showed in Table 1.
According to Tables1 and 2, it can be seen that the eigenvectors are different except the one corresponding to 1. In fact, by Table 1, the eigenvector corresponding to 7.4031 is v(1)=(0.5522,0.3534+0.4243i,-0.2650+0.5657i)T, and v(2)=(-0.2343-0.5000i,0.2343-0.5000i,0.6247)T by Table 2. It is easy to prove that rank(v(1),v(2))=1, that is, v(1)=eiθv(2). So the complex eigenvector after normalizing is not unique, which is the biggest difference from real eigenvector.Table 1
Example16—real fast algorithm for H1 (Algorithm 11).
Eigenvalue
Eigenvector
7.4031
0.5522
0.3534
+
0.4243
i
-
0.2650
+
0.5657
i
−5.4031
0.5522
0.3534
-
0.4243
i
-
0.2650
-
0.5657
i
1.0000
−0.6247
0.6247
-
0.4685Table 2
Example16—complex algorithm (eig) for eigenproblem H1.
Eigenvalue
Eigenvector
7.4031
-
0.2343
-
0.5000
i
0.2343
-
0.5000
i
0.6247
−5.4031
0.2343
-
0.5000
i
-
0.2343
-
0.5000
i
-
0.6247
1.0000
-
0.6247
-
0.0000
i
0.6247
+
0.0000
i
-
0.4685Example 17.
Now we give a little bigger matrixH2; we depict its eigenproblem in Tables 3 and 4, respectively:
(73)H2=[12+i3+i4+i5+i2-i23-1i463-i3+1i3154-i4131+i5-i651-i3].
According to tables, we can see that the eigenvectors are complex and it is easy to prove that the two eigenvectors corresponding to the same eigenvalue are linear dependent, so the eigenvectors, which are computed by the real fast algorithm (Algorithm11) span the same space as complex algorithm (eig).Table 3
Example17—the real fast algorithm for eigenproblem of H2.
Eigenvalue
Eigenvector
16.6334
0.4207
0.4417
+
0.1637
i
0.4177
+
0.1008
i
0.3343
+
0.0604
i
0.5220
+
0.1688
i
-
6.0472
0.4928
0.4942
-
0.0179
i
0.0460
+
0.0178
i
-
0.3757
-
0.1146
i
-
0.5966
+
0.0018
i
3.3730
-
0.1369
-
0.1327
+
0.0555
i
0.4347
+
0.1567
i
-
0.7694
-
0.2026
i
0.3370
+
0.0198
i
-
1.8356
-
0.5240
0.2953
-
0.2007
i
0.4425
+
0.4636
i
0.1982
+
0.1590
i
-
0.2902
-
0.1959
i
-
0.1236
0.5356
-
0.5467
-
0.2943
i
0.1736
+
0.3981
i
0.0803
+
0.1618
i
-
0.0590
-
0.3210
iTable 4
Example17—the real fast algorithm for eigenproblem of H2.
Eigenvalue
Eigenvector
16.6334
0.4003
-
0.1295
i
0.4706
+
0.0198
i
0.4285
-
0.0326
i
0.3367
-
0.0454
i
0.5486
-
6.0472
0.4928
+
0.0015
i
0.4942
-
0.0164
i
0.0459
+
0.0179
i
-
0.3753
-
0.1158
i
-
0.5966
3.3730
0.1367
-
0.0080
i
0.1293
-
0.0631
i
-
0.4431
-
0.1310
i
0.7800
+
0.1572
i
-
0.3376
-
1.8356
0.4343
-
0.2932
i
-
0.1325
+
0.3315
i
-
0.6262
-
0.1367
i
-
0.2532
-
0.0208
i
0.3502
-
0.1236
-
0.0968
+
0.5268
i
0.3882
-
0.4845
i
-
0.4229
+
0.0988
i
-
0.1736
+
0.0497
i
0.3263Example 18.
In this example we construct complex Hermitian matrices of the form(74)H3=diag(1,2,…,2*n)+[0-12-120-12-12⋱⋱⋱⋱-12-120]2n×2n+i*[0B-B0],
where B=(1/8)diag(eye(n)).  In tables we only give a few eigenvalues with largest modulus and compare it with complex algorithm (eig). Let n=10, 100; that is, the order of matrix H3 is 20, 200, respectively. The eigenvalues and errors, which are computed respectively by real fast algorithm and complex algorithm (eig), are showed in Tables 5 and 6. All the eigenvalues of matrix of order 20 are given, while only the 20th largest modulus eigenvalues of matrix order 200 are present. It is easy to prove that eigenvectors corresponding to the same eigenvalue are linear dependent. According to the tables, we have the same result as the complex algorithm (eig), but only the real operation is used.Table 5
Example17—the real fast algorithm for eigenproblem of H2.
Real fast algorithm
Complex algorithm (eig)
Error
2.022700124819288
e
+
001
2.022700124819288
e
+
001
7.105427357601002
e
-
015
1.902502989496266
e
+
001
1.902502989496268
e
+
001
1.421085471520200
e
-
014
1.800263600389322
e
+
001
1.800263600389323
e
+
001
1.065814103640150
e
-
014
1.700158594984117
e
+
001
1.700158594984119
e
+
001
2.131628207280301
e
-
014
1.600156256127496
e
+
001
1.600156256127496
e
+
001
0
1.500156225818778
e
+
001
1.500156225818779
e
+
001
8.881784197001252
e
-
015
1.400156222655317
e
+
001
1.400156222655317
e
+
001
3.552713678800501
e
-
015
1.300156061869106
e
+
001
1.300156061869103
e
+
001
2.842170943040401
e
-
014
1.200151108208256
e
+
001
1.200151108208256
e
+
001
3.552713678800501
e
-
015
1.100086192964745
e
+
001
1.100086192964745
e
+
001
0
9.999138070352572
e
+
000
9.999138070352565
e
+
000
7.105427357601002
e
-
015
8.998488917917468
e
+
000
8.998488917917442
e
+
000
2.664535259100376
e
-
014
7.998439381308970
e
+
000
7.998439381308969
e
+
000
8.881784197001252
e
-
016
6.998437773446823
e
+
000
6.998437773446838
e
+
000
1.509903313490213
e
-
014
5.998437741812211
e
+
000
5.998437741812216
e
+
000
5.329070518200751
e
-
015
4.998437438725036
e
+
000
4.998437438725052
e
+
000
1.598721155460225
e
-
014
3.998414050158813
e
+
000
3.998414050158812
e
+
000
8.881784197001252
e
-
016
2.997363996106763
e
+
000
2.997363996106775
e
+
000
1.199040866595169
e
-
014
1.974970105037331
e
+
000
1.974970105037334
e
+
000
3.774758283725532
e
-
015
7.729987518071178
e
-
001
7.729987518071133
e
-
001
4.440892098500626
e
-
015Table 6
Example17—the real fast algorithm for eigenproblem of H2.
Real fast algorithm
Complex algorithm (eig)
Error
2.002255917401552
e
+
002
2.002255917401532
e
+
002
2.074784788419493
e
-
012
1.990236230837739
e
+
002
1.990236230837733
e
+
002
6.252776074688882
e
-
013
1.980012299298954
e
+
002
1.980012299298937
e
+
002
1.676880856393836
e
-
012
1.970001799412475
e
+
002
1.970001799412465
e
+
002
9.663381206337363
e
-
013
1.960001565550514
e
+
002
1.960001565550502
e
+
002
1.193711796076968
e
-
012
1.950001562523487
e
+
002
1.950001562523486
e
+
002
1.136868377216160
e
-
013
1.940001562497741
e
+
002
1.940001562497715
e
+
002
2.586375558166765
e
-
012
1.930001562497571
e
+
002
1.930001562497559
e
+
002
1.165290086646564
e
-
012
1.920001562497562
e
+
002
1.920001562497558
e
+
002
4.547473508864641
e
-
013
1.910001562497567
e
+
002
1.910001562497556
e
+
002
1.193711796076968
e
-
012
1.900001562497573
e
+
002
1.900001562497560
e
+
002
1.364242052659392
e
-
012
1.890001562497554
e
+
002
1.890001562497561
e
+
002
7.389644451905042
e
-
013
1.880001562497553
e
+
002
1.880001562497555
e
+
002
1.989519660128281
e
-
013
1.870001562497555
e
+
002
1.870001562497559
e
+
002
3.694822225952521
e
-
013
1.860001562497558
e
+
002
1.860001562497558
e
+
002
8.526512829121202
e
-
014
1.850001562497552
e
+
002
1.850001562497558
e
+
002
5.400124791776761
e
-
013
1.840001562497563
e
+
002
1.840001562497560
e
+
002
2.842170943040401
e
-
013
1.830001562497566
e
+
002
1.830001562497558
e
+
002
8.526512829121202
e
-
013
1.820001562497562
e
+
002
1.820001562497559
e
+
002
3.410605131648481
e
-
013
1.810001562497546
e
+
002
1.810001562497558
e
+
002
1.193711796076968
e
-
012Example 19.
Complex Hermitian matrices are constructed by Hilbert matrix, and the MATLAB code is as follows:(75)n=10;(76)H=hilb(n);(77)H4=H+(tril(H)-diag(diag(H)))′i-(tril(H)-diag(diag(H))).
When n=10, the 6 largest modulus eigenvalues are present and 3 eigenvectors corresponding to the 3 largest modulus eigenvalues are given here:
(78)λ(H4)eig={2.137981825207990.51461933549659-0.49873888113430-0.188658336083330.17747332069964-0.08500068573387},λ(H4)Algorithm11={2.137981825207990.51461933549659-0.49873888113430-0.188658336083330.17747332069964-0.08500068573387},(79)V(eig)=[-0.1267-0.0857i-0.0156-0.1526i0.1174+0.0100i0.0529+0.2493i-0.0069+0.2192i-0.1508-0.1187i0.0647-0.2646i-0.2177-0.2291i0.1777+0.1521i-0.3229+0.1315i0.3057-0.0144i-0.0986-0.2811i0.2230+0.2709i-0.1554+0.3121i-0.0529+0.2981i0.2089-0.2644i-0.2960-0.2141i0.2578-0.2379i-0.2801-0.2083i0.2142-0.2740i-0.3615-0.0351i-0.2635+0.2411i0.2822+0.1932i0.2113+0.3073i0.1389+0.3183i-0.1273+0.3242i0.1640-0.3632i0.3304-0.3474-0.4016],(80)Q=[-0.6616-0.4938-0.3499-0.3850-0.2226i0.1121-0.3048i-0.0667+0.5036i-0.2449-0.2266i0.2977+0.0431i0.2830+0.3297i-0.1682-0.2102i0.1960+0.2453i0.3451+0.0933i-0.1203-0.1921i0.0573+0.3070i0.2918-0.0624i-0.0880-0.1754i-0.0550+0.2958i0.2086-0.1463i-0.0649-0.1608i-0.1329+0.2524i0.1273-0.1825i-0.0477-0.1480i-0.1816+0.1980i0.0581-0.1897i-0.0346-0.1368i-0.2082+0.1428i0.0028-0.1802i-0.0243-0.1270i-0.2190+0.0916i-0.0394-0.1614i].
We prove thatrank([V(:,k)Q(:,k)])=1 for k=1,2,3. It is very interesting that the eigenvectors computed from our algorithm are unique but not the complex algorithm (eig).Example 20.
At last we construct some rand Hermitian matrices of orderN; its MATLAB code is as follows.
(81)A=rand(n),SA=A′+A,(82)B=rand(n),SB=B-B′,(83)H=SA+i*SB.
Letting2*N=200:100:1000, denote
(84)max(relerr)=maxj|λjeig-λjSympl||λjSympl|,(85)min(relerr)=minj|λjeig-λjSympl||λjSympl|.
All the relative errors we compute are present by a factor of order𝒪(10d). The average relative errors for all examples of each dimension we computed for each example in Table 7 are denoted by average (relerr). The max, min, and average errors of r=Hvi-viλ are also considered in Table 8. According to two tables, all the “relerr” and “residual” are order 𝒪(10-13) and very small. Hence it can be concluded that eigenproblems computed by the real fast algorithm have high precision, and the method always converges.Table 7
Example20—relative error.
2
*
N
Max (relerr)
Min (relerr)
Average (relerr)
200
0.2645
e
-
013
0
0.0288
e
-
013
300
0.7872
e
-
013
0
0.0447
e
-
013
400
0.1949
e
-
013
0
0.0242
e
-
013
500
0.3536
e
-
013
0
0.0287
e
-
013
600
0.3297
e
-
012
0
0.0058
e
-
012
700
0.2149
e
-
013
0
0.0292
e
-
013
800
0.1773
e
-
011
0
0.0007
e
-
011
900
0.2841
e
-
012
0
0.0037
e
-
012
1000
0.6865
e
-
012
0
0.0046
e
-
012Table 8
Example20—residual r=Hvi-viλ.
2
*
N
ma
x
i
=
1
,
…
,
N
(
H
v
i
-
v
i
λ
)
mi
n
i
=
1
,
…
,
N
(
H
v
i
-
v
i
λ
)
Average
(
H
v
i
-
v
i
λ
)
200
0.1140
e
-
012
0.0113
e
-
012
0.0282
e
-
012
300
0.2110
e
-
012
0.0168
e
-
012
0.0420
e
-
012
400
0.2275
e
-
012
0.0225
e
-
012
0.0556
e
-
012
500
0.3608
e
-
012
0.0271
e
-
012
0.0682
e
-
012
600
0.6141
e
-
012
0.0324
e
-
012
0.0819
e
-
012
700
0.9271
e
-
012
0.0396
e
-
012
0.0954
e
-
012
800
0.8451
e
-
012
0.0452
e
-
012
0.1080
e
-
012
900
0.6920
e
-
012
0.0496
e
-
012
0.1175
e
-
012
1000
0.1218
e
-
011
0.0056
e
-
011
0.0132
e
-
011



---

*Source: 438320-2013-03-27.xml*