Oct 23, 2007 - FÃ¶hringer Ring 6, Dâ80805 Munich, Germany. PACS: 02.10.Ud, 02.10.Yn, 02.60.Dc. Jacobi-type iterative algorithms for the eigenvalue d...

0 downloads 8 Views 118KB Size

Routines for the diagonalization of complex matrices T. Hahna a

Max-Planck-Institut f¨ ur Physik F¨ohringer Ring 6, D–80805 Munich, Germany

MPP-2006-85 PACS: 02.10.Ud, 02.10.Yn, 02.60.Dc

Jacobi-type iterative algorithms for the eigenvalue decomposition, singular value decomposition, and Takagi factorization of complex matrices are presented. They are implemented as compact Fortran 77 subroutines in a freely available library.

V ∈ Cn¯ ×m ,

1. Introduction

n ¯ = min(m, n) ,

This note describes a set of routines for the eigenvalue decomposition, singular value decomposition, and Takagi factorization of a complex matrix. Unlike many other implementations, the current ones are all based on the Jacobi algorithm, which makes the code very compact but suitable only for small to medium-sized problems. Although distributed as a library, the routines are self-contained and can easily be taken out of the library and included in own code, removing yet another installation prerequisite. Owing to the small size of the routines (each about 3 kBytes source code) it is possible, in fact quite straightforward, to adapt the diagonalization routine to one’s own conventions rather than vice versa.

2.3. Takagi Factorization The Takagi factorization [1,2] is a less known diagonalization method for complex symmetric matrices A = AT ∈ Cn×n , U ∗ A U † = diag(σ1 , . . . , σn ) , U

−1

†

=U ,

(3)

σi > 0 .

Although outwardly similar to the eigenvalue decomposition of a Hermitian matrix, it is really the special case of an SVD with V = W ∗ , as it applies even to singular matrices. Note also that the left and right factors, U ∗ and U † , are in general not inverses of each other. One might think that the Takagi factorization is merely a scaled SVD. For example, the matrix 1 2 A= (4) 2 1

2.1. Eigenvalue Decomposition The eigenvalue decomposition of a nonsingular matrix A ∈ Cn×n takes the form (1)

The eigenvalues σi and transformation matrix U can be further characterized if A possesses certain properties:

has the SVD V T diag(σ1 , σ2 )W = !T √1 √1 √1 3 0 2 2 2 0 1 √1 √1 √1 − − 2 2 2

• A = A† (Hermitian): U −1 = U † , σi ∈ IR, • A = AT (symmetric): U −1 = U T . 2.2. Singular Value Decomposition The singular value decomposition (SVD) can be applied to an arbitrary matrix A ∈ Cm×n : V ∗ A W † = diag(σ1 , . . . , σn¯ ) ,

σi > 0 .

V consists of orthonormal row vectors, i.e. is also unitary for m = n.

2. Mathematical Background

U A U −1 = diag(σ1 , . . . , σn ) , σi ∈ C .

W −1 = W † ∈ Cn×¯n ,

(5) √1 2 √1 2

!

which can indeed be scaled to yield U T diag(σ1 , σ2 )U =

(2) 1

(6)

2

T. Hahn √1 2 √i 2

√1 2 − √i2

!T 3 0

0 1

√1 2 √i 2

√1 2 − √i2

!

But consider the matrix 0 1 A= 1 0 which has 1 0

to medium-sized problems the Jacobi method is a strong competitor, in particular as it has the following advantages:

.

• It is conceptually very simple and thus very compact. (7)

the SVD T 0 1 0 0 1 1 0 1 1 0

(8)

whereas its Takagi factorization is √1 2 − √i2

√1 2 √i 2

!T 1 0

0 1

√1 2 − √i2

√1 2 √i 2

!

. (9)

Although occurring less frequently than the eigenvalue decomposition, the Takagi factorization does have real applications in physics, e.g. in the diagonalization of mass matrices of Majorana fermions. 3. Jacobi Algorithm The Jacobi algorithm [3] consists of iteratively applying a basic 2×2 diagonalization formula until the entire n × n matrix is diagonal. It works in several ‘sweeps’ until convergence is achieved. In each sweep it rotates away the non-zero offdiagonal elements using the 2 × 2 algorithm. Every such rotation of course creates other non-zero off-diagonal elements. It can be shown, however, that the sum of the absolute values of the offdiagonal elements is reduced in each sweep. More precisely, the Jacobi method has quadratic convergence [4]. Convergence is in most cases achieved in 5−10 sweeps, which for an n × n matrix translates into (10−20)n3 multiply–add operations to obtain the eigenvalues only and (15−30)n3 operations including the eigenvectors ([5], cf. also Sect. 10). This compares with 23 n3 +30n2 operations for the Householder/QL algorithm when just the eigenvalues are sought and 43 n3 + 3n3 when also the eigenvectors are needed. For large n, the Jacobi algorithm is thus not the most efficient method. Nevertheless, for small

• It delivers the eigenvectors at little extra cost. • The diagonal values are accurate to machine precision and, in cases where this is mathematically meaningful, the vectors of the transformation matrix are always orthogonal, almost to machine precision. For the various diagonalization problems discussed before, only the core 2 × 2 diagonalization formula changes, whereas the surrounding Jacobi algorithm stays essentially the same. The famous Linear Algebra Handbook gives an explicit implementation of the Jacobi algorithm for real symmetric matrices [4], taking particular care to minimize roundoff errors through mathematically equivalent variants of the rotation formulas. The present routines are closely patterned on this procedure. For the Takagi factorization, the use of the Jacobi algorithm was first advocated in two conference papers [6,7] which give only few details, however. The two-sided Jacobi version of the SVD is used in [8]. Ref. [9] gives a more thorough overview of literature on the Jacobi method. 4. The 2 × 2 Formulas 4.1. Eigenvalue decomposition Using the ansatz c1 t1 c 1 U= −t2 c2 c2

(10)

the equation U A = diag(σ1 , σ2 ) U becomes 1 A12 , t1

(11)

1 A21 = A22 − t2 A12 . t2

(12)

σ1 = A11 + t1 A21 = A22 + σ2 = A11 −

Solving for t1 and t2 yields t1 =

A12 , ∆+D

t2 =

A21 , ∆+D

(13)

3

Routines for the diagonalization of complex matrices 1 (A11 − A22 ) , 2p D = ± ∆2 + A12 A21 . ∆=

(14) (15)

For the numerical stability it is best to choose the sign of D which gives t1,2 the larger denominator. This corresponds to taking the smaller rotation angle (< π/4). The diagonal values are σ1 = A11 + δ , σ2 = A22 − δ ,

δ=

A12 A21 . ∆+D

(16)

In order that U smoothly becomes unitary as A becomes Hermitian, we choose 1 c1 = c2 = √ , 1 + t1 t2

(17)

which guarantees a unit determinant. 4.2. Takagi Factorization Substituting the unitary ansatz c t c eiϕ U= , −t c e−iϕ c 1 c= √ , t ∈ IR , 1 + t2

σ ˜2 = e

−iϕ

σ2 ,

A˜11 = eiϕ A11 , A˜22 = e−iϕ A22 ,

A12 , ˜ ˜ ∆+D ˜ = 1 (A˜11 − A˜22 ) , ∆ 2q

˜ =± ∆ ˜ 2 + A2 . D 12

σ2 = A22 − t A12 e .

(28)

The assumption t ∈ IR fixes the phase ϕ. It re˜ have the same phase, i.e. quires that A12 and ∆ ˜ = (real number) · A12 . Since both eiϕ and its ∆ ˜ we try the ansatz conjugate appear in ∆, eiϕ = αA12 + βA∗12

(29)

and choose coefficients to make the A∗12 term in ˜ ∝ (αA11 − β ∗ A22 )A12 + ∆ (βA11 − α∗ A22 )A∗12

(18)

(30)

A∗11 A12 + A22 A∗12 . |A∗11 A12 + A22 A∗12 |

(31)

5. Singular Value Decomposition (19)

(20) (21)

(22) (23)

Comparing with Eqs. (11) and (12), the solution can be read off easily: t=

(27)

iϕ

eiϕ =

we arrive at 1 σ ˜1 = A˜11 + tA12 = A˜22 + A12 , t 1 ˜ ˜ σ ˜2 = A11 − A12 = A22 − tA12 . t

σ1 = A11 + t A12 e−iϕ ,

vanish. This is achieved by α = A∗11 and β = A22 which also makes the coefficient of A12 real. Thus,

into U ∗ A = diag(σ1 , σ2 ) U and introducing σ ˜1 = eiϕ σ1 ,

Again it is best for numerical stability to choose ˜ which gives the larger denominator the sign of D for t. The diagonal values become

(24) (25) (26)

We insert unitary parameterizations for the left and right transformation matrices X = V, W , cX tX c X , (32) X= −t∗X cX cX 1 cX = p (33) , tX ∈ C , 1 + |tX |2

into V ∗ A = diag(σ1 , σ2 ) W and by eliminating σ1,2 arrive at A12 + A22 t∗V = (A11 + A21 t∗V )tW , A21 +

A22 t∗W

= (A11 +

A12 t∗W )tV

.

(34) (35)

The solutions are evidently related through exchange of the off-diagonal elements. Explicitly, MV , ∆V + DV MV = A∗11 A21 + A22 A∗12 , 1 ∆V = (|A11 |2 − |A22 |2 + 2 |A12 |2 − |A21 |2 ) ,

tV =

(36) (37)

(38)

4

T. Hahn DV = ±

q ∆2V + |MV |2 ,

tW = tV |A12 ↔A21 .

(39) (40)

DV and DW share the same sign, which is chosen to yield the larger set of denominators for better numerical stability. The singular values become cV (A11 + t∗V A21 ) , cW cV (A22 − tV A12 ) , σ2 = cW

σ1 =

(41) (42)

If A ∈ Cm×n is not a square matrix, we consider two cases: For m > n, we make A square by padding it with zero-columns at the right. For zero right column, A12 = A22 = 0, Eq. (42) guarantees that it is σ2 that vanishes. That is, the above Jacobi rotation never moves a singular value into the zeroextended part of the matrix. All singular values automatically end up as the first min(m, n) diagonal values of the Jacobi-rotated A. For m < n, we apply this algorithm to AT and at the end exchange V and W . This is the least involved solution, as A has to be copied to temporary storage for zero-extension anyway. 6. Installation The Diag package can be downloaded from the URL http://www.feynarts.de/diag. After unpacking the tar file, the library is built with ./configure make To compile also the Mathematica executable, one needs to issue “make all” instead of just “make.” The generated files are installed into a platformdependent directory with “make install” and at the end one can do a “make clean” to remove intermediate files. The routines in the Diag library allocate space for intermediate results according to a preprocessor variable MAXDIM, defined in diag.h. This effectively limits the size of the input and output matrices but is necessary because Fortran 77 offers no dynamic memory allocation. Since the Jacobi algorithm is not particularly suited for large

problems anyway, the default value of 16 should be sufficient for most purposes. 7. Description of the Fortran Routines The general convention is that each matrix is followed by its leading dimension in the argument list, i.e. the m in A(m,n). In this way it is possible to diagonalize submatrices with just a different invocation. Needless to add, the leading dimension must be at least as large as the corresponding matrix dimension. 7.1. Hermitian Eigenvalue Decomposition Hermitian matrices are diagonalized with subroutine HEigensystem(n, A,ldA, d, U,ldU, sort) integer n, ldA, ldU, sort double complex A(ldA,n), U(ldU,n) double precision d(n) The arguments are as follows: • n (input), the matrix dimension. • A (input), the matrix to be diagonalized. Only the upper triangle of A needs to be filled and it is further assumed that the diagonal elements are real. Attention: the contents of A are not preserved. • d (output), the eigenvalues. • U (output), the transformation matrix. • sort (input), a flag that determines sorting of the eigenvalues: 0 = do not sort, 1 = sort into ascending order, −1 = sort into descending order. The ‘natural’ (unsorted) order is determined by the choice of the smaller rotation angle in each Jacobi rotation. 7.2. Symmetric Eigenvalue Decomposition The second special case is that of a complex symmetric matrix:

5

Routines for the diagonalization of complex matrices subroutine SEigensystem(n, A,ldA, d, U,ldU, sort) integer n, ldA, ldU, sort double complex A(ldA,n), U(ldU,n) double complex d(n) The arguments have the same meaning as for HEigensystem, except that A’s diagonal elements are not assumed real and sorting occurs with respect to the real part only. 7.3. General Eigenvalue Decomposition The general case of the eigenvalue decomposition is implemented in subroutine CEigensystem(n, A,ldA, d, U,ldU, sort) integer n, ldA, ldU, sort double complex A(ldA,n), U(ldU,n) double complex d(n) The arguments are as before, except that A has to be filled completely. 7.4. Takagi Factorization The Takagi factorization is invoked in almost the same way as SEigensystem: subroutine TakagiFactor(n, A,ldA, d, U,ldU, sort) integer n, ldA, ldU, sort double complex A(ldA,n), U(ldU,n) double precision d(n) The arguments are as for SEigensystem. Also here only the upper triangle of A needs to be filled. 7.5. Singular Value Decomposition The SVD routine has the form subroutine SVD(m, n, A,ldA, d, V,ldV, W,ldW, sort) integer m, n, ldA, ldV, ldW, sort double complex A(ldA,n) double complex V(ldV,m), W(ldW,n) double precision d(min(m,n)) with the arguments • m, n (input), the dimensions of A. • A (input), the m × n matrix of which the SVD is sought.

• d (output), the singular values. • V (output), the min(m, n)×m left transformation matrix. • W (output), the min(m, n) × n right transformation matrix. • sort (input), the sorting flag with values as above. 8. Description of the C Routines The C version consists merely of an include file CDiag.h which sets up the correct interfacing code for using the Fortran routines. In particular the usual problem of transposition1 between Fortran and C is taken care of. The arguments are otherwise as for Fortran. To treat complex numbers uniformly in C and C++, CDiag.h introduces the new double_complex type which is equivalent to complex

uses row-major storage for matrices, whereas Fortran uses column-major storage, i.e. a matrix is a vector of row vectors in C, and a vector of column vectors in Fortran.

6

T. Hahn

fcc -Iinclude myprogram.c -Llib -ldiag The fcc script is configured and installed during the build process and automatically adds the necessary flags for linking with Fortran code. 9. The Mathematica Interface The Mathematica version may not seem as useful as the Fortran library since Mathematica already has perfectly functional eigen- and singular value decompositions. The Takagi factorization is not available in Mathematica, however, and moreover the interface is ideal for trying out, interactively using, and testing the Diag routines. The Mathematica executable is loaded with Install["Diag"] and makes the following functions available:

scaling behaviour of the Jacobi algorithm, hence the curves appear essentially linear. The absolute time values should be taken for orientation only, as they necessarily reflect the CPU speed. For reference, the numbers used in the figure above were obtained on an AMD64-X2 CPU running at 3 GHz. Each point is actually the average from diagonalizing 106 random matrices, to reduce quantization effects in the time measurement. t/ms SVD

0.7 0.6 0.5 0.4

SEigensystem CEigensystem

0.3

TakagiFactor HEigensystem

0.2 0.1

• HEigensystem[A] computes the eigenvalue decomposition {d, U} of the Hermitian matrix A such that U.A == DiagonalMatrix[d].U. PSfrag • SEigensystem[A] computes the eigenvalue decomposition {d, U} of the symmetric matrix A such that U.A == DiagonalMatrix[d].U.

23 53 73 83 93 103

113

123

133

143

153

n3

The average number of sweeps needed to diagonalize the 106 random matrices to machine precision is plotted in the next figure. # sweeps

• CEigensystem[A] computes the eigenvalue decomposition {d, U} of the general matrix A such that U.A == DiagonalMatrix[d].U.

8

7 6

• TakagiFactor[A] computes the Takagi factorization {d, U} of the symmetric matrix A such that Conjugate[U].A == DiagonalMatrix[d].U. • SVD[A] computes the singular value decomposition {V, d, W} of the matrix A such that Conjugate[V].A == DiagonalMatrix[d].W. Note that these routines do not check whether the given matrix fulfills the requirements, e.g. whether it is indeed Hermitian. 10. Timings The following plot shows the time for diagonalizing a random matrix of various dimensions. Note that the abscissa is divided in units of the dimension cubed; this accounts for the anticipated

SVD HEigensystem SEigensystem TakagiFactor CEigensystem

5 4 3 2

2

3

5

7

9

11

13

15

11. Summary The Diag package contains Fortran subroutines for the eigenvalue decomposition, singular value decomposition, and Takagi factorization of a complex matrix. The Fortran library is supplemented by interfacing code to access the routines from C/C++ and Mathematica. The routines are based on the Jacobi algorithm. They are self-contained and quite compact, thus it should be straightforward to use them outside

Routines for the diagonalization of complex matrices of the library. All routines are licensed under the LGPL. Acknowledgements TH thanks the National Center for Theoretical Studies, Hsinchu, Taiwan, for warm hospitality during the time this work was carried out. REFERENCES 1. T. Takagi, Japanese J. Math. 1 (1927) 83. 2. R.A. Horn, C.A. Johnson, Matrix Analysis, Cambridge University Press, 1990, p. 201 f. 3. C.G.J. Jacobi, Crelle J. 30 (1846) 51. 4. H. Rutishauser, Contribution II/1, in: Handbook for Automatic Computation, ed. J.H. Wilkinson, C. Reinsch, Springer, 1971. 5. W.H. Press et al., Numerical Recipes in Fortran, 2nd ed., Cambridge University Press, 1992, Chapter 11. 6. L. De Lathauwer, B. De Moor, EUSIPCO 2002 conference proceedings, https://www.ensieta.fr/e3i2/intranet/ Confs/Eusipco02/articles/paper323.html 7. X. Wang, S. Qiao, PDPTA 2002 conference proceedings Vol. I, http://www.dcss.mcmaster.ca/ ∼qiao/publications/pdpta02.ps.gz 8. R.P. Brent, F.T. Luk, C.F. van Loan, J. of VLSI and Computer Systems 3 (1983) 242. G.E. Forsythe, P. Henrici, Trans. Amer. Math. Soc. 94 (1960) 1. 9. G.H. Golub, C.F. van Loan, Matrix Computations, 3rd ed., The Johns Hopkins University Press, 1996, Sect. 8.4.

7