A Real Algorithm for the Hermitian Eigenvalue Decomposition

Haesun Park ^{∗} Computer Science Department University of Minnesota Minneapolis, MN 55455

Vjeran Hari^{†} Department of Mathematics University of Zagreb, P.O.Box 187 41000 Zagreb, Yugoslavia

April, 1992

Abstract

Modifying complex plane rotations, we derive a new Jacobi-type algorithm for the Hermitian eigendecomposition, which uses only real arithmetic. When the fast-scaled rotations are incorporated, the new algorithm brings a substantial reduction in computational costs. The new method has the same convergence properties and parallelism as the symmetric Jacobi algorithm. Computational test results show that the new method produces accurate eigenvalues and eigenvectors and achieves great reduction in computational time.

AMS categories : 65F15 Key-words : Complex rotation, Fast rotation, Hermitian eigenvalue decomposition, Jacobi method

∗The work of this author was supported in part by the National Science Foundation grant CCR-8813493 and by the University of Minnesota Army High Performance Computing Research Center contract DAAL 03-89-C-0038.

†The work of this author was supported in part by the University of Minnesota Army High Performance Computing Research Center contract DAAL 03-89-C-0038.

1 Introduction

Recent interest in Jacobi methods is due to their adaptability for vector and parallel computing environments and their ability to deliver the eigenvalues and eigenvectors with higher accuracy [2]. On the other hand, their disadvantage, compared to the QR methods, lies in the larger computational costs. In this paper, we show how to reduce the computational costs for the complex Hermitian Jacobi eigenreduction. It is well known that substantial reduction can be obtained by using fast rotations [13]. Additional reduction can be accomplished by replacing the two-sided method by its onesided variation [15]. We show that it is possible to further reduce the computational cost by 25% for two-sided methods and 17% for one-sided methods. The standard Jacobi method for a Hermitian matrix H ∈ C^{n×n} typically generates a new matrix in each step by H′ = Q^{∗}HQ with a complex rotation matrix Q. Our modification replaces each complex rotation Q by ΦR where Φ is the same as the identity matrix except for one unit complex diagonal element and R is a real rotation. Our modification preserves the convergence properties of the standard Hermitian Jacobi method.

The paper is organized as follows. In Section 2, we present basic facts about complex rotations and our modifications. The convergence properties of the Jacobi method based on the modified complex rotations are discussed. In Section 3, we derive the new algorithm for the Hermitian eigenvalue decomposition, which uses only real arithmetic. We also show how the fast-scaled rotations are incorporated to achieve further reduction in the computational cost. The results of numerical tests are presented in Section 4.

2 Modifying the Complex Rotations

If H = A + ıB = (h_{ij}) ∈ C^{n×n} is a Hermitian matrix, where A = (a_{ij}) ∈ R^{n×n} and B = (b_{ij}) ∈ R^{n×n}, then A = A^{T} and B = −B^{T}. Let W = U + ıV ∈ C^{n×n} be a unitary matrix that satisfies

W^{∗} H W = Λ = diag(λ_{1}, …, λ_{n}) ∈ R^{n×n}.  (1)

The equation (1) is equivalent to

W^{T} H W = [ \begin{matrix}Λ \\ & Λ \end{matrix} ] ∈ R^{2n×2n},

where H = [ \begin{matrix}A & −B \\ B & A \end{matrix} ] ∈ R^{2n×2n} is symmetric and W = [ \begin{matrix}U & −V \\ V & U \end{matrix} ] ∈ R^{2n×2n} is

orthogonal. Throughout this paper we denote 2n × 2n real matrices in boldface.

The standard Jacobi method for the eigendecomposition of the Hermitian matrix H uses the unitary similarity transformations with complex rotations. A typical complex rotation Q in an (l, m) plane, l < m, is the same as the identity matrix I_{n}, except for the four nontrivial elements in its pivot submatrix (or the (l, m)-restriction), which we denote with the subscript lm:

Q_{lm} = [ \begin{matrix}cos φ & e^{ıα} sin φ \\ −e^{−ıα} sin φ & cos φ \end{matrix} ] .  (2)

Note that Q can be factored into

Q = Φ(l, α) R(φ) Φ(l, −α)

or

Q = Φ(m, −α) R(φ) Φ(m, α),

where Φ(r, α) is the same as the identity matrix except for e^{ıα} on its (r, r)th position, and R is a real rotation where the (l, m)–restriction is

R_{lm} = [ \begin{matrix}c & s \\ −s & c \end{matrix} ] , \begin{matrix}c = cos φ, \\ s = sin φ. \end{matrix}

A typical Jacobi step performs the transformation H′ = Q^{∗} H Q, where Q_{lm} is chosen to maximally decrease off(H)^{2} = ∑_{1≤i<j≤n} |h_{ij}|^{2} for the given (l, m), i.e., to make Q_{lm}^{∗} H_{lm} Q_{lm} diagonal. The above requirement yields

e^{−ıα}h_{lm} = ±|h_{lm}|, tan 2φ = ±\frac{2|h_{lm}|}{h_{mm} − h_{ll}} .  (3)

We will choose the + sign in formulas (3) since both + and − signs produce the same matrix Q. With the choice of the “inner” angle φ, we have

\begin{aligned}e^{ıα} & = & \frac{a_{lm} + ıb_{lm}}{|h_{lm}|} , \\ tan 2φ & = & \frac{2 |h_{lm}|}{a_{mm} − a_{ll}} , − \frac{π}{4} ≤ φ ≤ \frac{π}{4} . \end{aligned}  (4)

From the formulas (4), we get

a′_{ll} = a_{ll} − |h_{lm}| tan φ, a′_{mm} = a_{mm} + |h_{lm}| tan φ,  (5)

where H′ = A′ + ıB′, A′ = (a′_{ij}), and B′ = (b′_{ij}), which shows that the diagonal elements a_{ll} and a_{mm} can be updated by simple arithmetic.

Our modification consists of using the transformation, G = Φ(l, α) R(φ) or G = Φ(m, −α) R(φ) instead of Q. As with Q, in each stage, G_{lm} is chosen so that

G_{lm}^{∗}H_{lm}G_{lm} = [ \begin{matrix}c & −s \\ s & c \end{matrix} ] [ \begin{matrix}e^{−ıα} & 0 \\ 0 & 1 \end{matrix} ] [ \begin{matrix}h_{ll} & h_{lm} \\ h_{ml} & h_{mm} \end{matrix} ] [ \begin{matrix}e^{ıα} & 0 \\ 0 & 1 \end{matrix} ] [ \begin{matrix}c & s \\ −s & c \end{matrix} ]

or

G_{lm}^{∗}H_{lm}G_{lm} = [ \begin{matrix}c & −s \\ s & c \end{matrix} ] [ \begin{matrix}1 & 0 \\ 0 & e^{ıα} \end{matrix} ] [ \begin{matrix}h_{ll} & h_{lm} \\ h_{ml} & h_{mm} \end{matrix} ] [ \begin{matrix}1 & 0 \\ 0 & e^{−ıα} \end{matrix} ] [ \begin{matrix}c & s \\ −s & c \end{matrix} ]

is diagonal. Thus, the pivot submatrix of the modified complex unitary matrix has the form

G_{lm} = [ \begin{matrix}e^{ıα} cos φ & e^{ıα} sin φ \\ − sin φ & cos φ \end{matrix} ]

or

G_{lm} = [ \begin{matrix}cos φ & sin φ \\ −e^{−ıα} sin φ & e^{−ıα} cos φ \end{matrix} ] .

Conceptually, we diagonalize H_{lm} in two stages: first, find Φ to make Φ_{lm}^{∗}H_{lm}Φ_{lm} real, and then find the real rotation R to make R_{lm}^{T}(Φ_{lm}^{∗}H_{lm}Φ_{lm})R_{lm} diagonal. In fact, these conditions produce exactly the same requirements (3) on α and φ as in the standard complex rotation Q. The advantage of using G instead of Q is in the savings of up to 25 % in the total operations as we shall show. We note that Petkov [12] considered a variation of the Jacobi method for H, but our approach and aims are different. A similar idea has been used in constructing complex CORDIC [14].

The standard and the modified methods share the same convergence properties. The global convergence with the row-cyclic and the column-cyclic orderings follows directly from the results in [6] since the transformation G is a special elementary unitary matrix. It is a matter of combinatorial argument [3, 9, 10] to prove the convergence of the modified Hermitian Jacobi method under other frequently used parallel cyclic orderings. The quadratic convergence also follows easily from the results of Wilkinson [16], originally obtained for the row-cyclic symmetric Jacobi method. We presume that all the existing global and quadratic convergence proofs for the Hermitian Jacobi method hold or could be modified for our modified method. This is because all the essential inequalities in these proofs require only the moduli of sin φ and of matrix elements. For example the recent result for the symmetric Jacobi method [7] which extends to the case of multiple or close eigenvalues also holds for our method.

3 Fast Real Algorithms for the Hermitian Eigenvalue Decomposition

We derive two fast algorithms that require only real arithmetic, by incorporating fast rotations [1, 4, 5, 13, 17] for the fast-scaled standard complex rotation and its modified versions. Recent results due to Anda and Park [1] make it possible to dynamically control the underflow and overflow in the fast rotations. In deriving the algorithms we shall make use of the complex matrices H ∈ C^{n×n} and W ∈ C^{n×n}, and real matrices H ∈ R^{2n×2n} and W ∈ R^{2n×2n}.

3.1 Standard Algorithm

To derive the fast-scaled algorithm, we assume that

H = D H̃ D, Q = F Δ,

where D = diag(d_{1}, …, d_{n}), H̃ = Ã + ıB̃, Δ = diag(1, …, 1, δ_{1}, 1, …, 1, δ_{2}, 1, …, 1) and δ_{1}, δ_{2} appear at l-th and m-th positions. Initially, D = I_{n} and H̃ = H. Note that

H′ = Q^{∗} H Q = (F Δ)^{∗} D H̃ D (F Δ) = (D Δ)^{∗} F̃^{∗} H̃F̃ (D Δ)  (6)

where

F̃ = D F D^{−1}.

If we update D and H̃ in (6) separately as

D′ = D Δ, H̃′ = F̃^{∗} H̃F̃,  (7)

we can save about 50% of the multiplications for real rotations and 33% for complex rotations compared to the standard rotations, e.g., by taking δ_{1} = δ_{2} = c, where c ≥ 1/\sqrt{2}, at any step. This implies that matrices Δ and D are nonsingular with positive diagonal elements. If the eigenvectors are required, we can update the eigenvector matrix W in the same manner. Assuming that W = W̃ D, where W = I_{n} initially, we have W′ = W Q = W̃ D F Δ = W̃ F̃ D′. Hence W̃ is updated as

W̃′ = W̃ F̃,  (8)

and the unit eigenvectors are ultimately obtained by normalizing the columns of W̃. The value of tan φ is computed from

t = \frac{λ}{| μ | +\sqrt{λ^{2} + μ^{2}}}

where μ = a_{mm} − a_{ll} and λ = 2 · d_{l} · d_{m} · |h̃_{lm}| · sign(μ). We may use the fast rotation [1] that dynamically controls the underflow in the diagonal matrix D :

F̃_{lm} = \{ \begin{matrix}[ \begin{matrix}1 & 0 \\ γ_{1} & 1 \end{matrix} ] · [ \begin{matrix}1 & β_{1} \\ 0 & 1 \end{matrix} ] , & d_{l} ≥ d_{m} \\ [ \begin{matrix}1 & β_{2} \\ 0 & 1 \end{matrix} ] · [ \begin{matrix}1 & 0 \\ γ_{2} & 1 \end{matrix} ] , & d_{l} < d_{m} \end{matrix}  (9)

where

\begin{matrix}γ_{1} = −e^{−ıα}t\frac{d_{m}}{d_{l}} , & β_{1} = e^{ıα}tc^{2} \frac{d_{l}}{d_{m}} \\ γ_{2} = −e^{−ıα}tc^{2} \frac{d_{m}}{d_{l}} , & β_{2} = e^{ıα}t\frac{d_{l}}{d_{m}} \end{matrix}  (10)

with e^{ıα} = \frac{h̃_{lm}}{|h̃_{lm}|}, and the diagonal matrix D is updated as

[ \begin{matrix}d′_{l} & 0 \\ 0 & d′_{m} \end{matrix} ] = \{ \begin{matrix}[ \begin{matrix}d_{l}c & 0 \\ 0 & d_{m}/c \end{matrix} ] , & d_{l} ≥ d_{m} \\ [ \begin{matrix}d_{l}/c & 0 \\ 0 & d_{m}c \end{matrix} ] , & d_{l} < d_{m} \end{matrix} .  (11)

The diagonal components of DÃD which contain the computed eigenvalues of the current stage are updated as

\begin{aligned}a′_{ll} & = & a_{ll} − td_{l}d_{m}|h̃_{lm}|, \\ a′_{mm} & = & a_{mm} + td_{l}d_{m}|h̃_{lm}|, \end{aligned}  (12)

and the pivot elements are updated as

ã′_{lm} = 0, b̃′_{lm} = 0.  (13)

The relations (6)-(13) are summarized in the fast-scaled standard algorithm, Algorithm FS.

Algorithm FS: Given a Hermitian matrix H = A + ıB ∈ C^{n×n}, this algorithm finds the eigenvalues and the corresponding eigenvectors of H by real arithmetic.

1. d = (1, ..., 1); U = I_{n}; V = 0 ;

2. Repeat until convergence

3. Repeat for one sweep /* current pivot pair is (l, m) if a_{lm} ≠ 0 or b_{lm} ≠ 0 3.1 Compute F̃ : (9) 3.2 Update the 2 × 2 submatrices of A and B : (12),(13) 3.3 Update d_{l} and d_{m} : (11) 3.4 Update the l-th and m-th rows and columns of A and B via F̃

except for the 2 × 2 submatrix updated in 3.2 (only the upper triangular parts) : (7)

3.5 Update the l-th and m-th columns of U and V via F̃ : (8) /* The eigenvalues are the diagonal elements of A. /* The eigenvectors are the normalized columns of U + ıV.

Algorithm FS requires about 16n multiplications and additions (8n if eigenvectors are not required) for each step that consists of 3.1 - 3.5, for the given pair (l, m). For one sweep of n(n−1)/2 steps this amounts to 8n^{3} (4n^{3}) multiplications which is four times larger than those for the fast-scaled standard symmetric Jacobi method. The storage requirement is about 3n^{2} when eigenvectors are needed since the real and imaginary parts of each off-diagonal element h̃_{lm} can be stored in the location X(l, m) and X(m, l) of a real array X and the diagonal elements of H can be stored in the diagonal part of X in any step.

3.2 Modified Algorithms

As we defined in Section 2, the modified algorithm uses the transformation G = ΦR. One step of the modified method consists of two substeps:

H′ = Φ^{∗} H Φ and H′′ = R^{T} H′ R.

Both transformations induce unitary similarity transformations on the real matrix H ∈ R^{2n×2n}. The corresponding transformations are

H′ = R^{T} H R,  (14)

where R = I_{2n}, the identity matrix of order 2n, except for R_{l,n+l} = [ \begin{matrix}c_{α} & −s_{α} \\ s_{α} & c_{α} \end{matrix} ] when

Φ = Φ(l, α) and R_{m,n+m} = [ \begin{matrix}c_{α} & s_{α} \\ −s_{α} & c_{α} \end{matrix} ] when Φ = Φ(m, −α), and

H′′ = [ \begin{matrix}R^{T} \\ & R^{T} \end{matrix} ] H′ [ \begin{matrix}R \\ & R \end{matrix} ] ,  (15)

where R = I_{n} except for R_{lm} = [ \begin{matrix}c_{φ} & s_{φ} \\ −s_{φ} & c_{φ} \end{matrix} ], and c_{x} = cos x, s_{x} = sin x for real x.

For updating eigenvectors, we compute

W′ = W R  (16)

and

W′′ = W′ [ \begin{matrix}R \\ & R \end{matrix} ] .  (17)

To derive the fast-scaled algorithm we assume

\begin{matrix}H & = & D H̃ D, & W & = & W̃ D, \\ H′ & = & D′ H̃′ D′, & W′ & = & W̃′ D′, \\ H′′ & = & D′′ H̃′′ D′′, & W′′ & = & W̃′′ D′′, \end{matrix} \}  (18)

where

D = [ \begin{matrix}D \\ & D \end{matrix} ] , D′ = [ \begin{matrix}D′ \\ & D′ \end{matrix} ] , D′′ = [ \begin{matrix}D′′ \\ & D′′ \end{matrix} ] .

Thus, X̃ denotes a scaled factor of X for any matrix X. Also, we assume that

R = F Δ and R = F Δ  (19)

where F and F are fast rotations, and Δ and Δ are diagonal matrices.

First, we consider the transformations (14) and (16). The relations (18) and (19) imply that

H′ = R^{T} H R = Δ F^{T} D H̃ D F Δ = (D Δ) F̃^{T} H̃ F̃ (D Δ)

and W′ = W R = W̃ D F Δ = W̃ F̃ (D Δ) , where F̃ = D F D^{−1}. Since D_{l,n+l} = d_{l}I_{2} or D_{m,n+m} = d_{m}I_{2} we have F̃ = F. With D′ = D Δ and (18), we obtain

H̃′ = F^{T} H̃ F  (20)

and

W̃′ = W̃ F.  (21)

Initially, H̃ = H = [ \begin{matrix}A & −B \\ B & A \end{matrix} ] and the structure of H̃ is preserved after each step.

Thus, the transformation (20) has only to update the left 2n×n part [ Ã ; B̃ ] of H̃. The

parameter tan α is computed as

tan α = \frac{b_{lm}}{a_{lm}} = \frac{b̃_{lm}}{ã_{lm}} ,

where unlike φ, the angle α is not bounded. The fast rotation F, when Φ = Φ(l, α), is defined as

F_{l,n+l} = \{ \begin{matrix}[ \begin{matrix}1 & −t_{α} \\ t_{α} & 1 \end{matrix} ] , & |α| ≤ \frac{π}{4} \\ [ \begin{matrix}1/t_{α} & −1 \\ 1 & 1/t_{α} \end{matrix} ] , & \frac{π}{4} < |α| ≤ \frac{π}{2} \end{matrix}  (22)

and the diagonal matrix D is updated as

d′_{l} = \{ \begin{matrix}d_{l}c_{α}, & |α| ≤ \frac{π}{4} \\ d_{l}s_{α}, & \frac{π}{4} < |α| ≤ \frac{π}{2}. \end{matrix}  (23)

When Φ = Φ(m, −α), F_{m,n+m} can be defined and D is updated in a way which is similar to (22) and (23), respe ctively. The choice for Φ

Φ = \{ \begin{matrix}Φ(l, α) & if |d_{l}| ≥ |d_{m}| \\ Φ(m, −α) & otherwise \end{matrix}  (24)

will keep the deviation in the size of elements d_{j}’s small. After the first transformation, we have

ã′_{lm} = \{ \begin{matrix}\sqrt{ã_{lm}^{2} + b̃_{lm}^{2}}/c_{α} & , & |α| ≤ \frac{π}{4} \\ \sqrt{ã_{lm}^{2} + b̃_{lm}^{2}}/s_{α} & , & |α| > \frac{π}{4} \end{matrix} ,

b̃′_{lm} = 0,
a′_{ll} = a_{ll},
b′_{ll} = b_{ll}.

For the second transformation, from the relations (15),(18) and (19),

H′′ = [ \begin{matrix}R^{T} \\ & R^{T} \end{matrix} ] H′ [ \begin{matrix}R \\ & R \end{matrix} ]
= [ \begin{matrix}ΔF^{T} \\ & ΔF^{T} \end{matrix} ] D′ H̃′ D′ [ \begin{matrix}F Δ \\ & F Δ \end{matrix} ]
= ([ \begin{matrix}Δ \\ & Δ \end{matrix} ] D′) F̃^{T} H̃′ F̃ (D′ [ \begin{matrix}Δ \\ & Δ \end{matrix} ]) ,

where

F̃ ≡ D′ [ \begin{matrix}F \\ & F \end{matrix} ] D′^{−1} = [ F̃ ; F̃ ] .

Similarly, the relations (17),(18) and (19) imply

W′′ = W̃′ D′ [ \begin{matrix}F Δ \\ & F Δ \end{matrix} ] = W̃′ F̃ (D′ [ \begin{matrix}Δ \\ & Δ \end{matrix} ]) .

With

D′′ = D′ [ \begin{matrix}Δ \\ & Δ \end{matrix} ] ,  (25)

we obtain

H̃′′ = F̃^{T} H̃′ F̃, W̃′′ = W̃′ F̃.  (26)

Since

H̃′ = [ \begin{matrix}Ã′ & −B̃′ \\ B̃′ & Ã′ \end{matrix} ] , W̃′ = [ \begin{matrix}Ũ′ & −Ṽ′ \\ Ṽ′ & Ũ′ \end{matrix} ] ,

and

H̃′′ = [ \begin{matrix}Ã′′ & −B̃′′ \\ B̃′′ & Ã′′ \end{matrix} ] , W̃ = [ \begin{matrix}Ũ′′ & −Ṽ′′ \\ Ṽ′′ & Ũ′′ \end{matrix} ] ,

the relations (25) and (26) can be simplified to

\begin{matrix}D′′ & = & D′ Δ, \\ Ã′′ & = & F̃^{T} Ã′ F̃, \\ B̃′′ & = & F̃^{T} B̃′ F̃, \\ Ũ′′ & = & Ũ′ F̃, \\ Ṽ′′ & = & Ṽ′ F̃ \end{matrix} \}  (27)

where F̃ is given by (9) and (10) with α = 0. The parameter t_{φ} is computed from

t_{φ} = tan φ = \frac{λ}{| μ | +\sqrt{λ^{2} + μ^{2}}}  (28)

where μ = a_{mm} − a_{ll} and λ = 2 · d_{l} · d_{m} · \sqrt{ã_{lm}^{2} + b̃_{lm}^{2}} · sign(μ). The 2 × 2 submatrix in the (l, m) plane of A is updated as

\begin{matrix}a′′_{ll} & = & a_{ll} − t · d_{l} · d_{m} · \sqrt{ã_{lm}^{2} + b̃_{lm}^{2}} \\ a′′_{mm} & = & a_{mm} + t · d_{l} · d_{m} · \sqrt{ã_{lm}^{2} + b̃_{lm}^{2}} \\ ã′′_{lm} & = & 0 \end{matrix} \}  (29)

and the corresponding 2 × 2 block of B remains as zero after the first transformation. The fast-scaled modified real algorithm, Algorithm FM, also operates on the upper triangular parts of the arrays A and B and on the real arrays U and V for eigenvector updates. The matrices A and B are stored in the same way as in Algorithm FS.

Algorithm FM: Given a Hermitian matrix H = A + ıB ∈ C^{n×n}, this algorithm finds the eigenvalues and the corresponding eigenvectors of H by real arithmetic.

1 d = (1, ..., 1); U = I_{n}; V = 0 ; 2 Repeat until convergence

3 Repeat for one sweep /* current pivot pair is (l, m) if a_{lm} ≠ 0 or b_{lm} ≠ 0 3.1 Compute F : (22),(24) 3.2 Update d_{l} or d_{m} : (23),(24) 3.3 Update the l-th or m-th row and column of A and B via F

except for the diagonal and pivot elements (only the upper triangular parts) : (20)

3.4 Update the l-th or m-th columns of U and V via F : (21)

3.5 Compute F̃ : (28),(9),(10) with α = 0

3.6 Update the 2 × 2 submatrices of A and B : (29)

3.7 Update d_{l} and d_{m} : (11)

3.8 Update the l-th and m-th rows and columns of A and B via F̃ : (26) except for the 2 × 2 submatrix updated in 3.6 (only the upper triangular parts) : (26), (27)

3.9 Update the l-th and m-th columns of U and V via F̃ : (26), (27)

/* The eigenvalues are the diagonal components of A. /* The eigenvectors are the normalized columns of U + ıV.

Each step represented in 3.1-3.9 for a given pair (l, m) of the modified algorithm requires around 12n (6n if eigenvectors are not required) multiplications and the same number of additions. This amounts to 6n^{3} (3n^{3}) multiplications per cycle which is a reduction of 25% compared to the fast-scaled standard algorithm, Algorithm FS. Both algorithms can be arranged to operate on three real n × n arrays.

If the matrix H is Hermitian positive definite, then the one-sided algorithm that operates on the Cholesky factor L of H, where H = LL^{∗}, can further reduce the computational complexity [15]. At the expense of storing an additional vector of size n for the real diagonal elements, the diagonal elements can be updated simply (see the relation (5)), and only one off diagonal element has to be computed via the four real dot products for each pair (i, m) for the computation of rotation parameters for H. The one-sided version of Algorithm FS requires about 12n multiplications and 12n additions for each step. This is a 25% reduction compared to Algorithm FS. Also the storage requirement is less by one third since the complex matrix L and an additional vector of size n will provide the information for both eigenvectors and eigenvalues. The one-sided modified algorithm requires about 10n multiplications and 10n additions for each step. This is a 17% reduction compared to Algorithm FM. Details for the onesided Jacobi algorithm for positive definite symmetric matrices can be found in Veselić and Hari [15].

4 Computational Results

Computational tests were performed on an IBM compatible 386 DX machine with IN-TEL 80387 coprocessor. The TURBO PASCAL 6.0 compiler with two floating point arithmetics (the real and the extended precisions) were used. The real precision works with 12 decimal digit operands and the extended precision works with 18 decimal digit operands ^{1}, with the smallest and largest positive value around 10^{−38} and 10^{38} (10^{−4096} and 10^{4096}), respectively. The extended arithmetic was used for determining the errors in the computed eigenvalues.

The test matrices were generated in the following three ways. The real and imaginary parts of each element h_{lm} of a Hermitian matrix H are stored in the locations X(l, m) and X(m, l) of the real array X, respectively.

1machine precision ep≐ 9.95 · 10^{−13} (1.08 · 10^{−19})

Group 1: The elements X(i, j) of the array X are uniformly distributed integers in the interval [−100, 100]. This choice of X(i, j) ensures that the values, when stored as real and extended data types, are equal. The tests were performed for n = 10, 12, 14, 16, 18, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70.

Group 2: The elements X(i, j) are generated with uniformly distributed real numbers in the interval [−1, 1]. They are then scaled so that X(1, 1), X(2, 2), …, X(n, n) appear in decreasing order. The quotients |a_{i+1,i+1}/a_{ii}| are 0.1. The tests were performed for n = 10, 12, 14, 16, 18, 20.

Group 3: The matrices in this group are obtained from those of Group 2 using transformation P^{T} XP where P = [e_{n}, e_{n−1}, … e_{1}]. Thus, the quotients |a_{i+1,i+1}/a_{ii}| are 10.

The matrices from Group 2 and Group 3 are graded. In all procedures we used the row-cyclic ordering. Generic names of procedures are listed below.

S: This procedure is based on the standard Hermitian Jacobi algorithm, which uses the complex rotations Q_{lm} shown in (2).

M: This procedure is based on modified Hermitian Jacobi algorithm. It uses modified transformations G = Φ R.

FS: This procedure is based on Algorithm FS.

FM: This procedure is based on Algorithm FM.

For the eigenvalue tests we used the real and the extended precision arithmetics; for the eigenvector tests we used the extended precision arithmetic only.

4.1 Results for Eigenvalues

In each table, the following values are shown:

Relative error: max_{1≤i≤n} | \frac{λ̃_{i}−λ_{i}}{λ̃_{i}} | where λ_{i} (λ̃_{i}) is computed in real (extended) precision arithmetic.

d_{min}, d_{max} : d_{min} = min_{1≤i≤n} |d_{i}| and d_{max} = max_{1≤i≤n} |d_{i}|, where d_{i}’s are the elements in the final matrix D in the fast scaled procedures FS and FM.

Cycle time: the time required for the first sweep. The first sweep was chosen in order to minimize the number of rotations skipped.

For each case, the iteration was performed for 6 sweeps. Tables 1 and 2 show typical cases for matrices of order 18 and 70. The average of all 16 matrices from Group 1 is summarized in Table 3, where d_{min} and d_{max} are the minimum and maximum taken over all steps and over all matrices. The extremal eigenvalues for the whole group are 4.27266E+00, and 1.29756E+03.

name relative error

d_{min} d_{max}

cycle time

S  9.615E − 11  —  —  1.204
M  4.767E − 11  —  —  1.106
FS  7.352E − 11 8.043E − 01 1.193E + 00  0.930
FM  6.140E − 11 1.778E − 03 3.290E − 03  0.901

Table 1. n = 18, λ_{max} = 6.02561E + 02, λ_{min} = 8.25112E + 00,

name relative error

d_{min} d_{max}

cycle time

S  1.306E − 10  —  —  68.94
M  2.038E − 10  —  —  57.11
FS  4.847E − 11 7.075E − 01 1.443E + 00  46.67
FM  7.105E − 11 5.208E − 15 2.178E − 14  40.65

Table 2. n = 70, λ_{max} = 1.26662E + 03, λ_{min} = 1.58898E + 01

name

average relative error

\begin{matrix}minimum \\ d_{min} \end{matrix}

\begin{matrix}maximum \\ d_{max} \end{matrix}

average cycle time

S  1.583E − 10  —  —  2059
M  2.720E − 10  —  —  1719
FS  2.392E − 10 7.075E − 01 1.452E + 00  1443
FM  3.072E − 10 5.208E − 15 7.528E − 02  1414

Table 3. Statistics for matrices from Group 1

In Algorithm FM, d_{min} and d_{max} can become very small due to the first part of the algorithm that updates only one of d_{l} and d_{m} for a given pivot pair (l, m). But the choice of Φ and F based on the magnitude of d_{l} and d_{m} as shown in (3.19) and (3.6) makes the elements d_{j}’s stay in a small range.

Tables 4 and 5 show the results of matrices of order 20 from Groups 2 and 3, respectively. name relative error

d_{min} d_{max}

cycle time

S  4.672E − 11  —  —  1.602
M  8.719E − 11  —  —  1.446
FS  3.611E − 11 7.836E − 01 1.360E + 00  1.248
FM  3.611E − 11 5.613E − 02 9.489E − 02  1.201

Table 4. n = 20, λ_{max} = 3.79752E + 10, λ_{min} = 1.02180E − 08

name relative error

d_{min} d_{max}

cycle time

S  4.520E − 10  —  —  1.631
M  6.464E − 10  —  —  1.474
FS  3.473E − 10 6.875E − 01 1.317E + 00  1.257
FM  4.970E − 10 3.188E − 02 8.155E − 02  1.200

Table 5. n = 20, λ_{max} = 3.79752E + 10, λ_{min} = 1.02180E − 08

4.2 Results for Eigenvectors

We used the stopping criterion off(A) ≤ 100 ep ∥A∥. Since CPU time ratios remained almost the same as above we display only measures for accuracy of computed eigenvectors:

Maxdot: The maximum of eigenvector dot products max_{1≤i<j≤n} | (w_{i} | w_{j}) |,

Maxres: The maximum of residual norms max_{1≤i≤n} ∥A w − λ_{i} w∥.

For all matrices, but especially for those of smaller orders, accuracy of output data was within limits placed by rounding errors. The results are summarized in Tables 6-8.

name  maxdot  maxres
S  9.08730E − 19 7.07101E − 15
M  1.31694E − 18 8.92388E − 15
FS  8.25433E − 19 6.35780E − 15
FM  1.05848E − 18 7.77509E − 15

Table 6. Matrix presented in Table 2, n = 70, number of sweeps = 9

name  maxdot  maxres
S  1.70013E − 19 4.81035E − 09
M  3.37583E − 19 9.16319E − 09
FS  1.92113E − 19 7.46300E − 09
FM  4.06576E − 19 3.18438E − 09

Table 7. Matrix presented in Table 4, n = 20, number of sweeps = 3

name  maxdot  maxres
S  6.65694E − 19 1.42884E − 08
M  9.70510E − 19 1.69044E − 08
FS  8.24911E − 19 8.75591E − 09
FM  8.84352E − 19 1.96058E − 08

Table 8. Matrix presented in Table 5, n = 20, number of sweeps = 11

As the numerical tests suggest, the modified Jacobi algorithms compute the eigenvalues and eigenvectors of Hermitian matrices with approximately the same accuracy as the standard ones. Gain in computational effort becomes more significant for matrices of larger dimension, but it also depends on the computer.

Acknowledgment: The authors are indebted to Sanja Singer, Univ. of Zagreb, for the computational tests.

References

[1] A.A. Anda and H. Park: Fast plane rotations with dynamic scaling. To appear in SIAM J. Matrix Anal. Appl.

[2] J. Demmel and K. Veselić: Jacobi’s Method is more accurate than QR. Courant Institute of Mathematical Sciences, Department of Computer Science, New York University, 1989.

[3] P.J. Eberlein and H. Park: Efficient implementation of Jacobi algorithms and Jacobi sets on distributed memory architectures. Journal of Parallel and Distributed Computing, special issue on Algorithms for Hypercube Computers, 8, pp. 358-366, 1990.

[4] W.M. Gentleman: Least squares computations by Givens transformations without square roots. J. Inst. Maths. Applics, 12, 329-336, 1973.

[5] S. Hammarling: A note on modifications to the Givens plane rotation. J. Inst. Maths. Applics., 13, 215-218, 1974.

[6] V. Hari: On the convergence of cyclic Jacobi- like processes. Linear Algebra and its Appl. 81 (1986) 105-127.

[7] V. Hari: On sharp quadratic convergence bounds for the serial Jacobi methods. Tech. Report, University of Zagreb, 1990.

[8] M.R. Hestenes: Inversion of matrices by biorthogonalization and related results. J. SIAM 6 (1958) 51-90.

[9] F.T. Luk and H. Park: On parallel Jacobi orderings. SIAM J. Sci. Statist. Comput. 10 (1989) 18-26.

[10] F.T. Luk and H. Park: A Proof of convergence for two parallel Jacobi SVD algorithms. IEEE Trans. Comput. 38 (1989) 806-811.

[11] H. Park and P.J. Eberlein: Factored Jacobi-like algorithms for eigensystem computations on vector processors. Tech. Report TR90-11, Comp. Scie. Dept., Univ. of Minn., Feb., 1990.

[12] M. Petkov: Numerical Solution of Linear Algebraic Problems in the Complex Case. (Bulgarian language) Annuarie de ĺ Université de Sofia “Kliment Ohridski”, Tome 71, 2^{e} Partie, Faculté de Mathematiques et Mecanique 1976/77.

[13] W. Rath: Fast Givens rotations for orthogonal similarity transformations. Numer. Math. 40 (1982) 47-56

[14] A. van der Veen and E.F. Deprettere: Parallel VLSI matrix pencil algorithm for high resolution direction finding. IEEE Tran. on Signal Processing 39 (1991) 383-394.

[15] K. Veselić and V. Hari: A note on the one-sided Jacobi algorithm. Numer.Math. 56 (1989) 627-633.

[16] J.H. Wilkinson: Note on the quadratic convergence of the cyclic Jacobi process. Numer. Math. 4 (1962) 296-300.

[17] J.H. Wilkinson: Some Recent Advances in Numerical Linear Algebra. The State of Art in Numerical Analysis, ed. D.A.H. Jacobs, Academic press, New York, pp. 1-53.