Content deleted Content added
→Top: Dense algorithm |
WikiBuddy2 (talk | contribs) Link suggestions feature: 3 links added. |
||
| (15 intermediate revisions by 8 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]] (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
Line 6 ⟶ 7:
|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
|first1=G.H. |last1=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
Line 49 ⟶ 50:
:<math> \theta = \frac{\pi} {4} </math>
In order to optimize this effect, ''S''<sub>''ij''</sub> should be the [[off-diagonal element]] 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 ''S''.
Line 55 ⟶ 56:
== 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> . Let <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 <math> 2N := n(n-1) </math> off-diagonal 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>;
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
Line 64 ⟶ 65:
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
Line 83 ⟶ 84:
== Cost ==
Each Givens rotation can be done in <math> O(
=== 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'' − 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''<sup>3</sup>) 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''<sup>2</sup> 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 |
The
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 ==
Line 241 ⟶ 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-norm and spectral radius
Line 347 ⟶ 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> .
== References ==
| |||