Jacobi eigenvalue algorithm

Jacobi eigenvalue algorithm

The Jacobi eigenvalue algorithm is a numerical procedure for the calculation of all eigenvalues and eigenvectors of a real symmetric matrix.

Description

Let varphi in mathbb{R}, , 1 le k < l le n and let J(varphi, k, l) denote the n imes n matrix with components

: J_{kk} = cos varphi, ,!: J_{kl} = - sin varphi, ,!: J_{lk} = sin varphi, ,!: J_{ll} = cos varphi, ,!and: J_{ii} = 1 for the remaining diagonal elements: J_{ij} = 0 for the remaining off-diagonal elements

J(varphi, k, l) describes a rotation with angle φ in the plane spanned by the "k". and "l". unit vector.We have J(varphi, k, l)^T = J(varphi, k, l)^{- 1} , i.e. J(varphi, k, l) is orthogonal.

For a real symmetric matrix "S" let S(varphi, k, l) := J(varphi, k, l)^{-1}, S , J(varphi, k, l) . This matrix differs from "S" only in rows and columns "k" and "l" where ( if S ' = S(varphi, k, l), , c = cos varphi, , s = sin varphi for short) :: S_{kj}' = S_{jk}' = c , S_{kj} + s , S_{jl} , , S_{lj}' = S_{jl}' = - s ,S_{kj} + c, S_{jl} quad mbox{for} quad j e k, l

and

: S_{kk}' = c^2, S_{kk} + 2, s c ,S_{kl} + s^2, S_{ll} : S_{ll}' = s^2 ,S_{kk} - 2 s c, S_{kl} + c^2 , S_{ll} : S_{kl}' = S_{lk}' = (c^2 - s^2 ) , S_{kl} + s, c , (S_{ll} - S_{kk} ).

Define

: ||S ||_F^2 := mbox{trace} (S^ TS ) (the sum of squares of all components),

: sigma^2 := sum_{i = 1}^n S_{ii}^2 (the sum of squares of all diagonal components),

: Gamma(S )^2 := ||S ||_F^2 - sigma^2 (the sum of squares of all off-diagonal components).

||S||_F is the Frobenius norm of "S" and Γ ("S" ) is called off-diag norm since Gamma (S ) = 0 if and only if "S" is diagonal.

The Frobenius norm does not change under unitary rotations, so we have ||S'||_F = ||S||_F . If we choose φ according to

: an 2varphi = frac{2 S_{kl{S_{kk} - S_{ll

then S_{kl}' = 0 and one finds Gamma(S ')^2 = Gamma (S )^2 - 2 S_{kl}^ 2. Thus the off-norm decreases and hence "S' " becomes "more" diagonal than "S".In order to optimize this effect, S_{kl} should be the largest off-diagonal component. Such an element is called a pivot and the rotated matrix S^ J is called a Jacobi rotation.

The Jacobi eigenvalue method repeatedly performs Jacobi rotations until the matrix becomes almost diagonal. Then the elements in the diagonal are approximations of the (real) eigenvalues of "S".

Convergence

If p = S_{kl} is a pivot element, then by definition |S_{ij} | le |p| for 1 le i, j le n, i e j . Since "S " has exactly 2" N ":= "n "( "n " - 1) off-diag elements, we have p^2 le Gamma(S )^2 le 2 N p^2 or 2 p^2 ge Gamma(S )^2 / N . This implies Gamma(S^J )^2 le (1 - 1 / N ) Gamma (S )^2 or Gamma (S^ J ) le (1 - 1 / N )^{1 / 2} Gamma(S ) ,i.e. the sequence of Jacobi rotations converges at least linearly by a factor (1 - 1 / N )^{1 / 2} to a diagonal matrix.

A number of "N " Jacobi rotations is called a sweep; let S^{sigma} denote the result. The previous estimate yields: Gamma(S^{sigma} ) le (1 - 1 / N )^{N / 2} Gamma(S ) ,i.e. the sequence of sweeps converges at least linearly with a factor ≈ e ^{1 / 2} .

However the following result of Schönhage yields locally quadratic convergence. To this end let "S" have "m" distinct eigenvalues lambda_1, ... , lambda_m with multiplicities u_1, ... , u_m and let "d" > 0 be the smallest distance of two different eigenvalues. Let us call a number of

: N_S := frac{1}{2} n (n - 1) - sum_{mu = 1}^{m} frac{1}{2} u_{mu} ( u_{mu} - 1) le N

Jacobi rotations a Schönhage-sweep. If S^ s denotes the result then: Gamma(S^ s ) lesqrt{frac{n}{2} - 1} frac{gamma^2}{d - 2gamma}, quad gamma := Gamma(S ) .

Thus convergence becomes quadratic as soon as Gamma(S ) < d / (2 + sqrt{frac{n}{2} - 1}) .

Cost

Each Jacobi rotation can be done in "n" steps when the pivot element "p" is known. However the search for "p" requires inspection of all "N" ≈ ½ "n"² off-diag elements. We can reduce this to "n" steps too if we introduce an additional index array m_1, , ... , , , m_{n - 1} with the property that m_i is the index of the largest element in row "i", ("i" = 1, … , "n" - 1) of the current "S". Then ("k", "l") must be one of the pairs (i, m_i) . Since only columns "k" and "l" change, only m_k mbox{ and } m_l must be updated, which again can be done in "n" steps. Thus each rotation has O("n") cost and one sweep has O("n"³) cost which is equivalent to one matrix multiplication. Additionally the m_i must be initialized before the process starts, this can be done in "n"² steps.

Typically the Jacobi method converges within numerical precision after a small number of sweeps. Note that multiple eigenvalues reduce the number of iterations since N_S < N.

Algorithm

The following algorithm is a description of the Jacobi method in math-like notation.It calculates a vector "e" which contains the eigenvalues and a matrix "E" which contains the corresponding eigenvectors, i.e. e_i is an eigenvalue and the column E_i an orthonormal eigenvector for e_i , "i" = 1, … , "n".

procedure jacobi("S" &isin; R"n"×"n"; out "e" &isin; R"n"; out "E" &isin; R"n"×"n") var "i", "k", "l", "m", "state" &isin; N "s", "c", "t", "p", "y" &isin; R "ind" &isin; N"n" "changed" &isin; L"n" function maxind("k" &isin; N) &isin; N ! "index of largest element in row k" "m" := "k"+1 for "i" := "k"+2 to "n" do if │"S""ki"│ > │"S""km"then "m" := "i" endif endfor return "m" endfunc procedure update("k" &isin; N; "t" &isin; R) ! "update ek and its status" "y" := "e""k"; "e""k" := "y"+"t" if "changed""k" and ("y"="e""k") then "changed""k" := false; "state" := "state"−1 elsif (not "changed""k") and ("y"&ne;"e""k") then "changed""k" := true; "state" := "state"+1 endif endproc procedure rotate("k","l","i","j" &isin; N) ! "perform rotation of Sij, Skl" ┌ ┐ ┌ ┐┌ ┐ │"S""kl"│ │"c" −"s"││"S""kl"│ │ │ := │ ││ │ │"S""ij"│ │"s" "c"││"S""ij"│ └ ┘ └ ┘└ endproc ! "init e, E, and arrays ind, changed" "E" := "I"; "state" := "n" for "k" := 1 to "n" do "ind""k" := maxind("k"); "e""k" := "S""kk"; "changed""k" := true endfor while "state"&ne;0 do ! "next rotation" "m" := 1 ! "find index (k,l) of pivot p" for "k" := 2 to "n"−1 do if │"S""k" "ind""k"│ > │"S""m" "ind""m"then "m" := "k" endif endfor "k" := "m"; "l" := "ind""m"; "p" := "S""kl" ! "calculate c = cos &phi;, s = sin &phi;" "y" := ("e""l"−"e""k")/2; "t" := │"y"│+&radic;("p"2+"y"2) "s" := &radic;("p"2+"t"2); "c" := "t"/"s"; "s" := "p"/"s"; "t" := "p"2/"t" if "y"<0 then "s" := −"s"; "t" := −"t" endif "S""kl" := 0.0; update("k",−"t"); update("l","t") ! "rotate rows and columns k and l for "i" := 1 to "k"−1 do rotate("i","k","i","l") endfor for "i" := "k"+1 to "l"−1 do rotate("k","i","i","l") endfor for "i" := "l"+1 to "n" do rotate("k","i","l","i") endfor ! "rotate eigenvectors" for "i" := 1 to "n" do ┐ ┌ ┐┌ ┐ │"E""ki"│ │"c" −"s"││"E""ki"│ │ │ := │ ││ │ │"E""li"│ │"s" "c"││"E""li"│ └ ┘ └ ┘└ endfor ! "rows k, l have changed, update rows indk, indl" "ind""k" := maxind("k"); "ind""l" := maxind("l") loop endproc

Notes

1. The logical array "changed" holds the status of each eigenvalue. If the numerical value of e_k or e_l changes during an iteration, the corresponding component of "changed" is set to "true", otherwise to "false". The integer "state" counts the number of components of "changed" which have the value "true". Iteration stops as soon as "state" = 0. This means that none of the approximations e_1,, ..., , e_n has recently changed its value and thus it is not very likely that this will happen if iteration continues. Here it is assumed that floating point operations are optimally rounded to the nearest flointing point number.

2. The upper triangle of the matrix "S" is destroyed while the lower triangle and the diagonal are unchanged. Thus it is possible to restore "S" if necessary according to

for "k" := 1 to "n"−1 do ! "restore matrix S" for "l" := "k"+1 to "n" do "S""kl" := "S""lk" endfor endfor

3. The eigenvalues are not necessarily in descending order. This can be achieved by a simple sorting algorithm.

for "k" := 1 to "n"−1 do "m" := "k" for "l" := "k"+1 to "n" do if "e""l" > "e""m" then "m" := "l" endif endfor if "k" &ne; "m" then swap "e""m","e""k"; swap "E""m","E""k" endif endfor

4. The algorithm is written using matrix notation (1 based arrays instead of 0 based).

5. When implementing the algorithm, the part specified using matrix notation must be performed simultaneously.

Example

Let

Then "jacobi" produces the following eigenvalues and eigenvectors after 3 sweeps (19 iterations) :

e_1 = 2585.25381092892231

E_1 = egin{pmatrix}0.0291933231647860588\ -0.328712055763188997\ 0.791411145833126331\ -0.514552749997152907end{pmatrix}

e_2 = 37.1014913651276582

E_2 = egin{pmatrix}-0.179186290535454826\ 0.741917790628453435\ -0.100228136947192199\ -0.638282528193614892end{pmatrix}

e_3 = 1.4780548447781369

E_3 = egin{pmatrix}-0.582075699497237650\ 0.370502185067093058\ 0.509578634501799626\ 0.514048272222164294end{pmatrix}

e_4 = 0.1666428611718905

E_4 = egin{pmatrix}0.792608291163763585\ 0.451923120901599794\ 0.322416398581824992\ 0.252161169688241933end{pmatrix}

Applications for real symmetric matrices

When the eigenvalues (and eigenvectors) of a symmetric matrix are known, the followingvalues are easily calculated.

;Singular values:The singular values of a (square) matrix "A" are the square roots of the (non negative) eigenvalues of A^T A . In case of a symmetric matrix "S" we have of S^T S = S^2 , hence the singular values of "S" are the absolute values of the eigenvalues of "S"

;2-Norm and spectral radius:The 2-norm of a matrix "A" is the norm based on the euclidian vectornorm, i.e. the largest value || A x||_2 when x runs through all vectors with ||x||_2 = 1 . It is the largest singular value of "A". In case of a symmetric matrix it is largest absolute value of its eigenvectors and thus equal to its spectral radius.

;Condition number:The condition number of a nonsingular matrix "A" is defined as mbox{cond} (A) = || A ||_2 * || A^{-1} ||_2 . In case of a symmetric matrix it is the absolute value of the quotient of the largest and smallest eigenvalue. Matrices with large condition numbers can cause numerically unstable results : small perturbation can result in large errors. Hilbert matrices are the most famous ill-conditioned matrices. For example the fourth order Hilbert matrix has a condition of 15514, while for order 8 it is 2.7*10^8 .

;Rank:A matrix "A" has rank r if it has r columns which are linearily independent while the remaining columns are linearily dependent on these. Equivalently r is the dimension of the range of A. Furthermore it is the number of nonzero singular values.:In case of a symmetric matrix r is the number of nonzero eigenvalues. Unfortunately because of rounding errors numerical approximations of zero eigenvalues may not be zero (it may also happen that a numerical approximation is zero while the true value is not). Thus one can only calculate the "numerical" rank by making a decision which of the eigenvalues are close enough to zero.

;Pseudoinverse:The pseudo inverse of a matrix "A" is the unique matrix X = A^+ for which "AX" and "XA" are symmetric and for which "AXA = A, XAX = X" holds. If "A" is nonsingular, then ' A^+ = A^{-1} .:When procedure jacobi (S, e, E) is called, then the relation S = E^T mbox{Diag} (e) E holds where Diag("e") denotes the diagonal matrix with vector "e" on the diagonal. Let e^+ denote the vector where e_i is replaced by 1/e_i if e_i le 0 and by 0 if e_i is (numerically close to) zero. Since matrix "E" is orthogonal, it follows that the pseudo-inverse of S is given by S^+ = E^T mbox{Diag} (e^+) E .

;Least squares solution:If matrix "A" does not have full rank, there may not be a solution of the linear system "Ax = b". However one can look for a vector x for which || Ax - b ||_2 is minimal. The solution is x = A^+ b . In case of a symmetric matrix "S" as before, one has x = S^+ b = E^T mbox{Diag} (e^+) E b .

;Matrix exponential:From S = E^T mbox{Diag} (e) E one finds mbox{exp} S = E^T mbox{Diag} (mbox{exp} e) E where exp "e" is the vector where e_i is replaced by exp e_i . In the same way "f(S)" can be calculated in an obvious way for any (analytic) function "f".

;Linear Differential Equations:The differential equation " x' = A x', x(0) = a " has the solution " x(t) " = exp ("t A") "a". For a symmetric matrix "S" it follows that x(t) = E^T mbox{Diag} (mbox{exp} t e) E a . If a = sum_{i = 1}^n a_i E_i is the expansion of "a" by the eigenvectors of "S", then x(t) = sum_{i = 1}^n a_i exp(t e_i) E_i .:Let W^s be the vector space spanned by the eigenvectors of "S" which correspond to a negative eigenvalue and W^u analogously for the positive eigenvalues. If a in W^s then mbox{lim}_{t infty} x(t) = 0 i.e. the equilibrium point 0 is attractive to "x(t)". If a in W^u then mbox{lim}_{t infty} x(t) = infty , i.e. 0 is repulsive to "x(t)". W^s and W^u are called "stable" and "unstable" manifolds for "S". If "a" has components in both manifolds, then one component is attracted and one component is repelled. Hence "x(t)" approaches W^u as t infty .

Generalizations

The Jacobi Method has been generalized to complex hermitian matrices, general nonsymmetric real and complex matrices as well as block matrices.

Since singular values of a real matrix are the squares of the eigenvalues of the symmetric matrix S = A^T A it can also be used for the calculation of these values. For this case, the method is modified in such a way that "S" must not be explicitly calculated which reduces the danger of round-off errors. Note that J S J^T = J A^T A J^T = J A^T J^T J A J^T = B^T B with B , := J A J^T .

The Jacobi Method is also well suited for parallelism.

References

* Jacobi, " [http://resolver.sub.uni-goettingen.de/purl/?GDZPPN002144522 Über ein leichtes Verfahren, die in der Theorie der Säkularstörungen vorkommenden Gleichungen numerisch aufzulösen] ", Crelle's Journal 30 (1846), 51–94
* Schönhage, "Zur quadratischen Konvergenz des Jacobi-Verfahrens", Numer. Math. 6 (1964), 410–412
* Rutishauser, "The Jacobi Method for Real Symmetric Matrices", Numer. Math. 9 (1966), 1–10
* Sameh, "On Jacobi and Jacobi-like algorithms for a parallel computer", Math. Comp. 25 (1971), 579–590
* Veselic, "On a class of Jacobi-like procedures for diagonalizing arbitrary real matrices", Numer. Math. 33 (1979), 157–172
* Veselic, "A quadratically convergent Jacobi-like method for real matrices with complex conjugate eigenvalues", Numer. Math. 33 (1979), 425–435
* Shroff, "A parallel algorithm for the eigenvalues and eigenvectors of a general complex matrix", Numer. Math. 58 (1991), 779–805

External links

* [http://www.maths.lancs.ac.uk/~gilbert/m306c/node17.html Jacobi's Method]
* [http://math.fullerton.edu/mathews/n2003/JacobiMethodMod.html Jacobi Iteration for Eigenvectors]
* [http://groups.google.com/group/sci.math.num-analysis/msg/8282d0d412f72d2e Matlab implementation of Jacobi algorithm that avoids trigonometric functions]


Wikimedia Foundation. 2010.

Игры ⚽ Нужно сделать НИР?

Look at other dictionaries:

  • Eigenvalue algorithm — In linear algebra, one of the most important problems is designing efficient and stable algorithms for finding the eigenvalues of a matrix. These eigenvalue algorithms may also find eigenvectors. Contents 1 Characteristic polynomial 2 Power… …   Wikipedia

  • Jacobi rotation — In numerical linear algebra, a Jacobi rotation is a rotation, Q k ℓ, of a 2 dimensional linear subspace of an n dimensional inner product space, chosen to zero a symmetric pair of off diagonal entries of an n × n real symmetric matrix, A , when… …   Wikipedia

  • Carl Gustav Jacob Jacobi — Carl Jacobi Carl Gustav Jacob Jacobi Born December 10, 1804(1804 …   Wikipedia

  • List of numerical analysis topics — This is a list of numerical analysis topics, by Wikipedia page. Contents 1 General 2 Error 3 Elementary and special functions 4 Numerical linear algebra …   Wikipedia

  • Singular value decomposition — Visualization of the SVD of a 2 dimensional, real shearing matrix M. First, we see the unit disc in blue together with the two canonical unit vectors. We then see the action of M, which distorts the disk to an ellipse. The SVD decomposes M into… …   Wikipedia

  • List of mathematics articles (J) — NOTOC J J homomorphism J integral J invariant J. H. Wilkinson Prize for Numerical Software Jaccard index Jack function Jacket matrix Jackson integral Jackson network Jackson s dimensional theorem Jackson s inequality Jackson s theorem Jackson s… …   Wikipedia

  • Pidgin code — In computer programming, pidgin code is a mixture of several programming languages in the same program, or pseudocode that is a mixture of a programming language with natural language descriptions. Hence the name: the mixture is a programming… …   Wikipedia

  • List of mathematics articles (E) — NOTOC E E₇ E (mathematical constant) E function E₈ lattice E₈ manifold E∞ operad E7½ E8 investigation tool Earley parser Early stopping Earnshaw s theorem Earth mover s distance East Journal on Approximations Eastern Arabic numerals Easton s… …   Wikipedia

  • Nonlinear eigenproblem — A nonlinear eigenproblem is a generalization of an ordinary eigenproblem to equations that depend nonlinearly on the eigenvalue. Specifically, it refers to equations of the form: where x is a vector (the nonlinear eigenvector ) and A is a matrix… …   Wikipedia

  • Maurice Clint — was professor of Computer Science at Queen s University Belfast, Northern Ireland. His research interests include parallel algorithms for numerical linear algebra and the use of formal methods in the development of parallel and distributed… …   Wikipedia

Share the article and excerpts

Direct link
Do a right-click on the link above
and select “Copy Link”