Jacobi eigenvalue algorithm: Difference between revisions

Revision as of 21:45, 7 January 2011

Latest revision as of 04:54, 22 August 2025

(92 intermediate revisions by 66 users not shown)

Line 1:

{{short description|Numerical linear algebra algorithm}}

In [[numerical linear algebra]], the '''Jacobi eigenvalue algorithm''' is an [[iterative method]] for the calculation of the [[eigenvalue]]s and [[eigenvector]]s of a [[real number|real]] [[symmetric matrix]]. It is named after [[Carl Gustav Jacob Jacobi]], who first proposed the method in 1846,<ref>{{cite journal

In [[numerical linear algebra]], the '''Jacobi eigenvalue algorithm''' is an [[iterative method]] for the calculation of the [[eigenvalue]]s and [[eigenvector]]s of a [[real number|real]] [[symmetric matrix]] (a process known as [[Matrix diagonalization#Diagonalization|diagonalization]]). It is named after [[Carl Gustav Jacob Jacobi]], who first proposed the method in 1846,<ref>{{cite journal

|last=Jacobi |first=C.G.J. |authorlink=Carl Gustav Jacob Jacobi

|url=http://resolvergdz.sub.uni-goettingen.de/purldms/load/img/?PID=GDZPPN002144522

|title=Über ein leichtes Verfahren, die in der Theorie der Säkularstörungen vorkommenden Gleichungen numerisch aufzulösen

|journal=[[Crelle's Journal]]

|volume=1846 |issue=30 |year=1846 |pages=51–94

|language=German

|doi=10.1515/crll.1846.30.51 |s2cid=199546177 |url-access=subscription }}</ref> but it only became widely used in the 1950s with the advent of computers.<ref>{{cite journal

|firstfirst1=G.H. |lastlast1=Golub |author1-link=Gene H. Golub

|first2=H.A. |last2=van der Vorst|author2-link=Henk van der Vorst

|title=Eigenvalue computation in the 20th century

|journal=Journal of Computational and Applied Mathematics

|volume=123 |issue=1-21–2 |year=2000 |pages=35&ndash;65

|doi=10.1016/S0377-0427(00)00413-1

|doi-access=free}}</ref>

This algorithm is inherently a [[dense matrix]] algorithm: it draws little or no advantage from being applied to a sparse matrix, and it will destroy sparseness by creating fill-in. Similarly, it will not preserve structures such as being [[band matrix|banded]] of the matrix on which it operates.

== Description ==

Let ''A''<math>S</math> be a symmetric matrix, and ''<math>G''&nbsp;=&nbsp;''G''(''i'',''j'',''&\theta;'')</math> be a [[Givens rotation|Givens rotation matrix]]. Then:

:<math>AS'=G^\top AS G \, </math>

is symmetric and [[similar matrix(linear algebra)|similar]] to ''A''<math>S</math>.

Furthermore, ''A&<math>S^\prime;''</math> has entries:

:<math>\begin{align}

AS'_{ii} &= c^2\, A_S_{ii} - 2\, s c \,A_S_{ij} + s^2\, A_S_{jj} \\

AS'_{jj} &= s^2 \,A_S_{ii} + 2 s c\, A_S_{ij} + c^2 \, A_S_{jj} \\

AS'_{ij} &= AS'_{ji} = (c^2 - s^2 ) \, A_S_{ij} + s c \, (A_S_{ii} - A_S_{jj} ) \\

AS'_{ik} &= AS'_{ki} = c \, A_S_{ik} - s \, A_S_{jk} & k \ne i,j \\

AS'_{jk} &= AS'_{kj} = s \, A_S_{ik} + c \, A_S_{jk} & k \ne i,j \\

AS'_{kl} &= A_S_{kl} &k,l \ne i,j

\end{align}</math>

where ''<math>s''&nbsp;=&nbsp;\sin(''&\theta;'')</math> and ''<math>c''&nbsp;=&nbsp;\cos(''&\theta;'')</math>.

AsSince they<math>G</math> areis similarorthogonal, ''A''<math>S</math> and ''A''&<math>S^\prime;</math> have the same [[Frobenius norm]] <math>||·\cdot||<sub>F_F</submath> (the square-root sum of squares of all components), however we can choose ''&<math>\theta;''</math> such that ''A''&prime;<submath>''S^\prime_{ij''}=0</submath>&nbsp;=&nbsp;0, in which case ''A''&<math>S^\prime;</math> has a larger sum of squares on the diagonal:

:<math> AS'_{ij} = \cos(2\theta) A_S_{ij} + \tfrac{1}{2} \sin(2\theta) (A_S_{ii} - A_S_{jj}) </math>

Set this equal to 0, and rearrange:

:<math> \tan(2\theta) = \frac{2 A_S_{ij}}{A_S_{jj} - A_S_{ii}} </math>

if <math> A_S_{jj} = A_S_{ii} </math>

:<math> \theta = \frac{\pi} {4} </math>

In order to optimize this effect, ''AS''_''ij'' should be the largest [[off-diagonal componentelement]] with the largest [[absolute value]], called the ''pivot''.

The Jacobi eigenvalue method repeatedly [[Jacobi rotation|performs rotations]] until the matrix becomes almost diagonal. Then the elements in the diagonal are approximations of the (real) eigenvalues of ''AS''.

== Convergence ==

If <math> p = S_{kl} </math> is a [[pivot element]], then by definition <math> |S_{ij} | \le |p| </math> for <math> 1 \le i, j \le n, i \ne j</math> . SinceLet ''<math>\Gamma(S)^2</math> ''denote the sum of squares of all off-diagonal entries of <math>S</math>. Since <math>S</math> has exactly 2''<math> N2N '':= ''n ''( ''n '' - 1) </math> off-diagdiagonal elements, we have <math> p^2 \le \Gamma(S )^2 \le 2 N p^2 </math> or <math> 2 p^2 \ge \Gamma(S )^2 / N </math> . Now <math>\Gamma(S^J)^2=\Gamma(S)^2-2p^2</math>. This implies

<math> \Gamma(S^J )^2 \le (1 - 1 / N ) \Gamma (S )^2 </math> or <math> \Gamma (S^ J ) \le (1 - 1 / N )^{1 / 2} \Gamma(S ) </math> ,;

i.e.that is, the sequence of Jacobi rotations converges at least linearly by a factor <math> (1 - 1 / N )^{1 / 2} </math> to a [[diagonal matrix]].

A number of ''<math> N ''</math> Jacobi rotations is called a sweep; let <math> S^{\sigma} </math> denote the result. The previous estimate yields

: <math> \Gamma(S^{\sigma} ) \le \left(1 - \frac{1 / }{N} \right)^{N / 2} \Gamma(S ) </math>,;

i.e.that is, the sequence of sweeps converges at least linearly with a factor ≈ <math> e ^{1 / 2}</math> .

However the following result of [[Arnold Schönhage|Schönhage]]<ref>{{cite journal

|last=Schönhage |first=A.|authorlink=Arnold Schönhage

|title=Zur quadratischen Konvergenz des Jacobi-Verfahrens

|journal=Numerische Mathematik

|volume=6 |issue=1 |year=1964 |pages=410–412

|language=German

|doi=10.1007/BF01386091 |idmr={{MR|174171}}

|s2cid=118301078

}}</ref> yields locally quadratic convergence. To this end let ''S'' have ''m'' distinct eigenvalues <math> \lambda_1, ... , \lambda_m </math> with multiplicities <math> \nu_1, ... , \nu_m </math> and let ''d'' > 0 be the smallest distance of two different eigenvalues. Let us call a number of

}}</ref> yields locally quadratic convergence. To this end let ''S'' have ''m'' distinct eigenvalues <math> \lambda_1, ... , \lambda_m </math> with multiplicities <math> \nu_1, ... , \nu_m </math> and let ''d'' > 0 be the smallest distance of two different eigenvalues. Let us call a number of

: <math> N_S := \frac{1}{2} n (n - 1)}{2} - \sum_{\mu = 1}^{m} \frac{1}{2} \nu_{\mu} (\nu_{\mu} - 1) \le N </math>

Jacobi rotations a Schönhage-sweep. If <math> S^ s </math> denotes the result then

: <math> \Gamma(S^ s ) \le\sqrt{\frac{n}{2} - 1} \left(\frac{\gamma^2}{d - 2\gamma}\right), \quad \gamma := \Gamma(S ) </math> .

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

== Cost ==

Each JacobiGivens rotation can be done in ''<math> O(n'') </math> steps when the pivot element ''p'' is known. However the search for ''p'' requires inspection of all ''N''&nbsp;≈&nbsp;½{{sfrac|1|2}}&nbsp;''n''² off-diagdiagonal elements., We canwhich reducemeans this tosearch ''n''dominates stepsthe toooverall ifcomplexity weand introduce an additional index array <math> m_1, \, \dots \, , \, m_{n - 1} </math> withpushes the propertycomputational that <math> m_i </math> is the indexcomplexity of thea largest elementsweep in row ''i'', (''i'' = 1, … , ''n''&nbsp;&minus;&nbsp;1) of the currentclassical ''S''.Jacobi Thenalgorithm (''k'', ''l'') must be one of the pairsto <math> O(i, m_in^4) </math> . SinceCompeting onlyalgorithms columns ''k'' and ''l'' change, onlyattain <math> m_k \mbox{ and } m_l </math> must be updated, which again can be done in ''n'' steps. Thus each rotation has O(''n'') cost and one sweep has O(''n''^{^3}) cost which is equivalent to one matrix multiplication. Additionally the <math> m_i </math> complexity mustfor bea initializedfull before the process starts, this can be done in ''n''² stepsdiagonalisation.

=== Caching row maximums ===

We can reduce the complexity of finding the pivot element from O(''N'') to O(''n'') if we introduce an additional index array <math> m_1, \, \dots \, , \, m_{n - 1} </math> with the property that <math> m_i </math> is the index of the largest element in row ''i'', (''i'' = 1, ..., ''n''&nbsp;&minus;&nbsp;1) of the current ''S''. Then the indices of the pivot (''k'', ''l'') must be one of the pairs <math> (i, m_i) </math>. Also the updating of the index array can be done in O(''n'') [[average-case complexity]]: First, the maximum entry in the updated rows ''k'' and ''l'' can be found in O(''n'') steps. In the other rows ''i'', only the entries in columns ''k'' and ''l'' change. Looping over these rows, if <math> m_i </math> is neither ''k'' nor ''l'', it suffices to compare the old maximum at <math> m_i </math> to the new entries and update <math> m_i </math> if necessary. If <math> m_i </math> should be equal to ''k'' or ''l'' and the corresponding entry decreased during the update, the maximum over row ''i'' has to be found from scratch in O(''n'') complexity. However, this will happen on average only once per rotation. Thus, each rotation has O(''n'') and one sweep O(''n''³) average-case complexity, which is equivalent to one [[matrix multiplication]]. Additionally the <math> m_i </math> must be initialized before the process starts, which 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 <math>N_S < N</math>.

=== Cyclic and parallel Jacobi ===

An alternative approach is to forego the search entirely, and simply have each sweep pivot every off-diagonal element once, in some predetermined order. It has been shown that this ''cyclic Jacobi'' attains quadratic convergence,<ref>{{cite journal |last=Wilkinson |first=J.H. |authorlink=James H. Wilkinson |title=Note on the Quadratic Convergence of the Cyclic Jacobi Process |journal=Numerische Mathematik |date=1962 |volume=6 |pages=296–300|doi=10.1007/BF01386321 }}</ref><ref>{{cite journal |last1=van Kempen |first1=H.P.M. |title=On Quadratic Convergence of the Special Cyclic Jacobi Method |journal=Numerische Mathematik |date=1966 |volume=9 |pages=19–22|doi=10.1007/BF02165225 }}</ref> just like the classical Jacobi.

The opportunity for parallelisation that is particular to Jacobi is based on combining cyclic Jacobi with the observation that Givens rotations for [[disjoint sets]] of indices commute, so that several can be applied in parallel. Concretely, if <math> G_1 </math> pivots between indices <math> i_1, j_1 </math> and <math> G_2 </math> pivots between indices <math> i_2, j_2 </math>, then from <math> \{i_1,j_1\} \cap \{i_2,j_2\} = \varnothing </math> follows <math> G_1 G_2 = G_2 G_1 </math> because in computing <math> G_1 G_2 A </math> or <math> G_2 G_1 A </math> the <math> G_1 </math> rotation only needs to access rows <math> i_1, j_1 </math> and the <math> G_2 </math> rotation only needs to access rows <math> i_2, j_2 </math>. Two processors can perform both rotations in parallel, because no matrix element is accessed for both.

Partitioning the set of index pairs of a sweep into classes that are pairwise disjoint is equivalent to partitioning the edge set of a [[complete graph]] into [[Matching (graph theory)|matching]]s, which is the same thing as [[edge colouring]] it; each colour class then becomes a round within the sweep. The minimal number of rounds is the [[chromatic index]] of the complete graph, and equals <math>n</math> for odd <math>n</math> but <math>n-1</math> for even <math>n</math>. A simple rule for odd <math>n</math> is to handle the pairs <math> \{i_1,j_1\} </math> and <math> \{i_2,j_2\} </math> in the same round if <math> i_1+j_1 \equiv i_2+j_2 \textstyle\pmod{n} </math>. For even <math>n</math> one may create <math> n-1 </math> rounds <math> k = 0, 1, \dotsc, n-2 </math> where a pair <math> \{i,j\} </math> for <math> 1 \leqslant i < j \leqslant n-1 </math> goes into round <math> (i+j) \bmod (n-1) </math> and additionally a pair <math> \{i,n\} </math> for <math> 1 \leqslant i \leqslant n-1 </math> goes into round <math> 2i \bmod (n-1) </math>. This brings the time complexity of a sweep down from <math> O(n^3) </math> to <math> O(n^2) </math>, if <math> n/2 </math> processors are available.

A round would consist of each processor first calculating <math>(c,s)</math> for its rotation, and then applying the rotation from the left (rotating between rows). Next, the processors [[synchronise]] before applying the transpose rotation from the right (rotating between columns), and finally synchronising again. A matrix element may be accessed by two processors during a round, but not by both during the same half of this round.

Further parallelisation is possible by dividing the work for a single rotation between several processors, but that might be getting too fine-grained to be practical.

== 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; that is, i.e. <math> e_i </math> is an eigenvalue and the column <math> E_i </math> an orthonormal eigenvector for <math> e_i </math> , ''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;''d'', ''r'' ∈ '''R'''

''ind'' &isin;∈ '''N'''^''n''

''changed'' &isin;∈ '''L'''^''n''

'''function''' maxind(''k'' &isin;∈ '''N''') &isin;∈ '''N''' ! ''index of largest off-diagonal element in row k''

''m'' := ''k''+1

'''for''' ''i'' := ''k''+2 '''to''' ''n'' '''do'''

Line 103 ⟶ 122:

'''endfunc'''

'''procedure''' update(''k'' &isin;∈ '''N'''; ''t'' &isin;∈ '''R''') ! ''update e_k 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 S_ij, S_kl''

┌ ┐ ┌ ┐┌ ┐

│''S''_''kl''│ │''c'' −''s''││''S''_''kl''│

Line 121 ⟶ 140:

''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≠0 '''do''' ! ''next rotation''

''m'' := 1 ! ''find index (k,l) of pivot p''

'''for''' ''k'' := 2 '''to''' ''n''−1 '''do'''

'''if''' │''S''_{''k''&nbsp;''ind''}_{_''k''}│ > │''S''_{''m''&nbsp;''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; ''td'' := │''y''│+&radic;√(''p''²+''y''²)

''sr'' := &radic;√(''p''²+''td''²); ''c'' := ''td''/''sr''; ''s'' := ''p''/''sr''; ''t'' := ''p''²/''td''

'''if''' ''y''<0 '''then''' ''s'' := −''s''; ''t'' := −''t'' '''endif'''

''S''_''kl'' := 0.0; update(''k'',−''t''); update(''l'',''t'')

Line 139 ⟶ 158:

'''for''' ''i'' := 1 '''to''' ''n'' '''do'''

┌ ┐ ┌ ┐┌ ┐

│''E''_''kiik''│ │''c'' −''s''││''E''_''kiik''│

│ │ := │ ││ │

│''E''_''liil''│ │''s'' ''c''││''E''_''liil''│

└ ┘ └ ┘└ ┘

'''endfor'''

! ''rowsupdate k,all l havepotentially changed, update rows ind_k, ind_li''

''ind'for'_{''k ''i''} := maxind(1 ''k'to''' ''n'' '''do'''); ''ind''_''li'' := maxind(''li'') '''endfor'''

'''loop'''

'''endproc'''

Line 151 ⟶ 170:

=== Notes ===

1. The logical array ''changed'' holds the status of each eigenvalue. If the numerical value of <math> e_k </math> or <math> e_l </math> 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 <math> e_1,\, ...\, , e_n </math> 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 flointingfloating 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'''

Line 162 ⟶ 183:

'''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'''

Line 178 ⟶ 204:

Let

<math>

S = \begin{pmatrix} 4 & -30 & 60 & -35 \\ -30 & 300 & -675 & 420 \\ 60 & -675 & 1620 & -1050 \\ -35 & 420 & -1050 & 700 \end{pmatrix}

[[Image:smat.png]]

</math>

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

Line 206 ⟶ 233:

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

</math>

<math>

Line 215 ⟶ 241:

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

</math>

these link can explain the example

http://www.scribd.com/doc/24953275/Jacobi-Transformation-Used-for-Eigenproblems-for-Symmetric-Matrix

== Applications for real symmetric matrices ==

Line 226 ⟶ 248:

;Singular values

:The singular values of a (square) matrix ''<math>A''</math> are the square roots of the (non -negative) eigenvalues of <math> A^T A </math>. In case of a symmetric matrix ''<math>S''</math> we have of <math> S^T S = S^2 </math>, hence the singular values of ''<math>S''</math> are the absolute values of the eigenvalues of ''<math>S''</math>.

;2-Normnorm and spectral radius

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

;Condition number

:The condition number of a nonsingular matrix ''<math>A''</math> is defined as <math> \mbox{cond} (A) = \| A \|_2 \| A^{-1}\|_2 </math>. 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 matrix|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&nbsp;&times;×&nbsp;10⁸.

;Rank

:A matrix ''<math>A''</math> has rank ''<math>r''</math> if it has ''<math>r''</math> columns that are linearly independent while the remaining columns are linearilylinearly dependent on these. Equivalently, ''<math>r''</math> is the dimension of the range of&nbsp;''<math>A''</math>. 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.

;Pseudo-inverse

;Pseudoinverse

:The pseudo inverse of a matrix ''<math>A''</math> is the unique matrix <math> X = A^+ </math> for which ''<math>AX''</math> and ''<math>XA''</math> are symmetric and for which ''<math>AXA = A, XAX = X'' </math> holds. If ''<math>A''</math> is nonsingular, then '<math> A^+ = A^{-1} </math>.

:When procedure jacobi (S, e, E) is called, then the relation <math> S = E^T \mbox{Diag} (e) E </math> holds where Diag(''e'') denotes the diagonal matrix with vector ''e'' on the diagonal. Let <math> e^+ </math> denote the vector where <math> e_i </math> is replaced by <math> 1/e_i </math> if <math> e_i \le 0 </math> and by 0 if <math> e_i </math> is (numerically close to) zero. Since matrix ''E'' is orthogonal, it follows that the pseudo-inverse of S is given by <math> S^+ = E^T \mbox{Diag} (e^+) E </math>.

;Least squares solution

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

;Matrix exponential

:From <math> S = E^T \mbox{Diag} (e) E </math> one finds <math> \exp S = E^T \mbox{Diag} (\exp e) E </math> where exp&nbsp;''<math>e''</math> is the vector where <math> e_i </math> is replaced by <math> \exp e_i </math>. In the same way, ''<math>f''(''S'')</math> can be calculated in an obvious way for any (analytic) function ''<math>f''</math>.

;Linear differential equations

:The differential equation ''<math>x'&nbsp;'' =&nbsp;''Ax'', ''x''(0) = ''a''</math> has the solution ''<math>x''(''t'') = \exp(''t&nbsp;A''tA)&nbsp;''a''</math>. For a symmetric matrix ''<math>S'' </math>, it follows that <math> x(t) = E^T \mbox{Diag} (\exp t e) E a </math>. If <math> a = \sum_{i = 1}^n a_i E_i </math> is the expansion of ''<math>a''</math> by the eigenvectors of ''<math>S''</math>, then <math> x(t) = \sum_{i = 1}^n a_i \exp(t e_i) E_i </math>.

:Let <math> W^s </math> be the vector space spanned by the eigenvectors of ''<math>S''</math> which correspond to a negative eigenvalue and <math> W^u </math> analogously for the positive eigenvalues. If <math> a \in W^s </math> then <math> \mbox{lim}_{t \rightarrow \infty} x(t) = 0 </math>; i.e.that is, the equilibrium point 0 is attractive to ''<math>x''(''t'')</math>. If <math> a \in W^u </math> then <math> \mbox{lim}_{t \rightarrow \infty} x(t) = \infty </math>; that is, i.e. 0 is repulsive to ''<math>x''(''t'')</math>. <math> W^s </math> and <math> W^u </math> are called ''stable'' and ''unstable'' manifolds for ''<math>S''</math>. If ''<math>a''</math> has components in both manifolds, then one component is attracted and one component is repelled. Hence ''<math>x''(''t'')</math> approaches <math> W^u </math> as <math> t \to \infty </math>.

== Julia implementation ==

The following code is a straight-forward implementation of the mathematical description of the Jacobi eigenvalue algorithm in the [[Julia (programming language)|Julia programming language]].

<syntaxhighlight lang="julia">

using LinearAlgebra, Test

function find_pivot(Sprime)

n = size(Sprime,1)

pivot_i = pivot_j = 0

pivot = 0.0

for j = 1:n

for i = 1:(j-1)

if abs(Sprime[i,j]) > pivot

pivot_i = i

pivot_j = j

pivot = abs(Sprime[i,j])

end

end

end

return (pivot_i, pivot_j, pivot)

end

# in practice one should not instantiate explicitly the Givens rotation matrix

function givens_rotation_matrix(n,i,j,θ)

G = Matrix{Float64}(I,(n,n))

G[i,i] = G[j,j] = cos(θ)

G[i,j] = sin(θ)

G[j,i] = -sin(θ)

return G

end

# S is a symmetric n by n matrix

n = 4

sqrtS = randn(n,n);

S = sqrtS*sqrtS';

# the largest allowed off-diagonal element of U' * S * U

# where U are the eigenvectors

tol = 1e-14

Sprime = copy(S)

U = Matrix{Float64}(I,(n,n))

while true

(pivot_i, pivot_j, pivot) = find_pivot(Sprime)

if pivot < tol

break

end

θ = atan(2*Sprime[pivot_i,pivot_j]/(Sprime[pivot_j,pivot_j] - Sprime[pivot_i,pivot_i] )) / 2

G = givens_rotation_matrix(n,pivot_i,pivot_j,θ)

# update Sprime and U

Sprime .= G'*Sprime*G

U .= U * G

end

# Sprime is now (almost) a diagonal matrix

# extract eigenvalues

λ = diag(Sprime)

# sort eigenvalues (and corresponding eigenvectors U) by increasing values

i = sortperm(λ)

λ = λ[i]

U = U[:,i]

# S should be equal to U * diagm(λ) * U'

@test S ≈ U * diagm(λ) * U'

</syntaxhighlight>

== Generalizations ==

Line 257 ⟶ 354:

Since singular values of a real matrix are the square roots of the eigenvalues of the symmetric matrix <math> S = A^T A</math> 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 error]]s. Note that <math> J S J^T = J A^T A J^T = J A^T J^T J A J^T = B^T B </math> with <math> B \, := J A J^T </math> .

The Jacobi Method is also well suited for parallelism.

== References ==

Line 265 ⟶ 360:

== Further reading ==

{{refbegin}}

*{{Citation | last1=Press | first1=WH | last2=Teukolsky | first2=SA | last3=Vetterling | first3=WT | last4=Flannery | first4=BP | year=2007 | title=Numerical Recipes: The Art of Scientific Computing | edition=3rd | publisher=Cambridge University Press | ___location=New York | isbn=978-0-521-88068-8 | chapter=Section 11.1. Jacobi Transformations of a Symmetric Matrix | chapter-url=http://apps.nrbook.com/empanel/index.html#pg=570 | access-date=2011-08-13 | archive-date=2011-08-11 | archive-url=https://web.archive.org/web/20110811154417/http://apps.nrbook.com/empanel/index.html#pg=570 | url-status=dead }}

* {{cite journal

|last=Rutishauser |first=H.

Line 270 ⟶ 366:

|journal=Numerische Mathematik

|volume=9 |issue=1 |year=1966 |pages=1–10

|doi=10.1007/BF02165223 |idmr={{MR|1553948}}

|s2cid=120520713

}}

}}

* {{cite journal

|last=Sameh |first=A.H.

Line 277 ⟶ 374:

|journal=[[Mathematics of Computation]]

|volume=25 |issue=115 |year=1971 |pages=579–590

|idmr={{MR|297131}} {{JSTOR| jstor = 2005221}} | doi = 10.1090/s0025-5718-1971-0297131-6

|doi-access=free }}

}}

* {{cite journal

|last=Shroff |first=Gautam M.

|title=A parallel algorithm for the eigenvalues and eigenvectors of a general complex matrix

|journal=Numerische Mathematik

|volume=58 |issue=1 |year=1991 |pages=779–805

|doi=10.1007/BF01385654 |mr=1098865

|citeseerx=10.1.1.134.3566

|s2cid=13904356

}}

* {{cite journal

|last=Veselić |first=K.

Line 284 ⟶ 390:

|journal=Numerische Mathematik

|volume=33 |issue=2 |year=1979 |pages=157–172

|doi=10.1007/BF01399551 |idmr={{MR|549446}}

|s2cid=119919630

}}

}}

* {{cite journal

|lastlast1=Veselić |firstfirst1=K.

|last2=Wenzel |first2=H. J.

|title=A quadratically convergent Jacobi-like method for real matrices with complex eigenvalues

|journal=Numerische Mathematik

|volume=33 |issue=4 |year=1979 |pages=425–435

|doi=10.1007/BF01399324 |idmr={{MR|553351}}

|s2cid=119554420

}}

}}

* {{cite journal

* Yousef Saad: "Revisiting the (block) Jacobi subspace rotation method for the symmetric eigenvalue problem", Numerical Algorithms, vol.92 (2023), pp.917-944. https://doi.org/10.1007/s11075-022-01377-w .

|last=Shroff |first=Gautam M.

|title=A parallel algorithm for the eigenvalues and eigenvectors of a general complex matrix

|journal=Numerische Mathematik

|volume=58 |issue=1 |year=1991 |pages=779–805

|doi=10.1007/BF01385654 |id={{MR|1098865}}

}}

{{refend}}

== External links ==

*[https://groups.google.com/group/sci.math.num-analysis/msg/8282d0d412f72d2e Matlab implementation of Jacobi algorithm that avoids trigonometric functions]

*[http://www.maths.lancs.ac.uk/~gilbert/m306c/node17.html Jacobi's Method]

*[https://github.com/jewettaij/jacobi_pd C++11 implementation]

*[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]

*[http://code.google.com/p/parallel-jacobi/ C++ implementation of the Jacobi algorithm in parallel]

{{Numerical linear algebra}}

Line 313 ⟶ 413:

[[Category:Numerical linear algebra]]

[[Category:Articles with example pseudocode]]

[[de:Jacobi-Verfahren (Eigenwerte)]]

[[mhr:Шкеыҥжым муашлан Йакобин алгоритмже]]

[[nl:Jacobi-eigenwaarde algoritme]]

[[uk:Метод обертання Якобі]]