5.1_Jacobi_迭代法

5.1 Jacobi迭代法

该算法的基本思想是通过一系列的 Jacobi 旋转 JkJ_{k} , 将 AA 正交相似于一个对角矩阵, 即

A(0)=A,A(k+1)=JkA(k)JkT,k=0,1,,A ^ {(0)} = A, \quad A ^ {(k + 1)} = J _ {k} A ^ {(k)} J _ {k} ^ {\mathsf {T}}, \quad k = 0, 1, \ldots ,

A(k)A^{(k)} 收敛到一个对角矩阵,其中 JkJ_{k} 为Jacobi旋转,通常选取 JkJ_{k} 为Givens变换,即

Jk=G(ik,jk,θk)=[11cosθksinθksinθkcosθk11]J _ {k} = G \left(i _ {k}, j _ {k}, \theta_ {k}\right) = \left[ \begin{array}{c c c c c c c c c} 1 & & & & & & & & \\ & \ddots & & & & & & & \\ & & 1 & & & & & & \\ & & & \cos \theta_ {k} & & \sin \theta_ {k} & & & \\ & & & & \ddots & & & & \\ & & & - \sin \theta_ {k} & & \cos \theta_ {k} & & & \\ & & & & & & 1 & & \\ & & & & & & & \ddots & \\ & & & & & & & & 1 \end{array} \right]

易知, 在 A(k)A^{(k)} 两边分别左乘 JkJ_{k} 和右乘 JkTJ_{k}^{\mathsf{T}} 时, 只会修改 A(k)A^{(k)} 的第 iki_{k} 和第 jkj_{k} 行, 以及第 iki_{k} 和第 jkj_{k} 列.

由于 A(k)A^{(k)} 是对称矩阵, 由下面的引理可知, 通过选取适当的 θk\theta_{k} , 可以将 A(k)(ik,jk)A^{(k)}(i_k,j_k)A(k)(jk,ik)A^{(k)}(j_k,i_k) 同时化为0.

引理5.1设 AR2×2A\in \mathbb{R}^{2\times 2} 是对称矩阵,则存在Givens变换 GR2×2G\in \mathbb{R}^{2\times 2} ,使得 GAGGAG^{\top} 为对角矩阵

(板书)

证明. 设

A=[abbc],G=[cosθsinθsinθcosθ],A = \left[ \begin{array}{c c} a & b \\ b & c \end{array} \right], \quad G = \left[ \begin{array}{c c} \cos \theta & \sin \theta \\ - \sin \theta & \cos \theta \end{array} \right],

GAGT=[cosθsinθsinθcosθ][abbc][cosθsinθsinθcosθ]T=[acos2θ+csin2θ+bsin2θ12(ca)sin2θ+bcos2θ12(ca)sin2θ+bcos2θasin2θ+ccos2θbsin2θ]\begin{array}{l} G A G ^ {\mathsf {T}} = \left[ \begin{array}{l l} \cos \theta & \sin \theta \\ - \sin \theta & \cos \theta \end{array} \right] \left[ \begin{array}{l l} a & b \\ b & c \end{array} \right] \left[ \begin{array}{l l} \cos \theta & \sin \theta \\ - \sin \theta & \cos \theta \end{array} \right] ^ {\mathsf {T}} \\ = \left[ \begin{array}{l l} a \cos^ {2} \theta + c \sin^ {2} \theta + b \sin 2 \theta & \frac {1}{2} (c - a) \sin 2 \theta + b \cos 2 \theta \\ \frac {1}{2} (c - a) \sin 2 \theta + b \cos 2 \theta & a \sin^ {2} \theta + c \cos^ {2} \theta - b \sin 2 \theta \end{array} \right] \\ \end{array}

12(ca)sin2θ+bcos2θ=0\frac{1}{2} (c - a)\sin 2\theta +b\cos 2\theta = 0 ,可得

ac2b=cot2θ=1tan2θ2tanθ.\frac {a - c}{2 b} = \cot 2 \theta = \frac {1 - \tan^ {2} \theta}{2 \tan \theta}.

解得

tanθ=sign(τ)τ+1+τ2,τ=ac2b.\tan \theta = \frac {\operatorname {s i g n} (\tau)}{| \tau | + \sqrt {1 + \tau^ {2}}}, \qquad \tau = \frac {a - c}{2 b}.

故引理结论成立.

为了使得 A(k)A^{(k)} 收敛到一个对角矩阵, 其非对角线元素必须趋向于 0. 记 off(A)\operatorname{off}(A) 为所有非对角线元素的平方和, 即

off(A)=ijaij2=AF2i=1naii2,\operatorname {o f f} (A) = \sum_ {i \neq j} a _ {i j} ^ {2} = \| A \| _ {F} ^ {2} - \sum_ {i = 1} ^ {n} a _ {i i} ^ {2},

我们的目标就是使得 off(A)\operatorname{off}(A) 尽快趋于 0.

引理5.2 设 A=[aij]Rn×nA = [a_{ij}] \in \mathbb{R}^{n \times n} 是对称矩阵, A^=[a^ij]=JAJ,J=G(i,j,θ)\hat{A} = [\hat{a}_{ij}] = JAJ^{\top}, J = G(i,j,\theta) , 其中 θ\theta 的选取使得 a^ij=a^ji=0\hat{a}_{ij} = \hat{a}_{ji} = 0 , 则

off(A^)=off(A)2aij2,ij.\mathrm {o f f} (\hat {A}) = \mathrm {o f f} (A) - 2 a _ {i j} ^ {2}, \quad i \neq j.

(板书)

证明. 记 A=[a1,a2,,an]A = [a_{1}, a_{2}, \ldots, a_{n}] . 令 A~=JA=[a~ij]n×n\tilde{A} = JA = [\tilde{a}_{ij}]_{n \times n} . 由于 JJ 是正交阵, 故

Jak2=ak2,k=1,2,,n.\| J a _ {k} \| _ {2} = \| a _ {k} \| _ {2}, \quad k = 1, 2, \dots , n.

JJ 左乘 aka_{k} 时, 只影响其第 ii 和第 jj 个元素的值, 故由 Jai2=ai2\| Ja_i\|_2 = \| a_i\|_2Jaj2=aj2\| Ja_j\|_2 = \| a_j\|_2 可得

a~ii2+a~ji2=aii2+aji2,a~ij2+a~jj2=aij2+ajj2.(5.1)\tilde {a} _ {i i} ^ {2} + \tilde {a} _ {j i} ^ {2} = a _ {i i} ^ {2} + a _ {j i} ^ {2}, \quad \tilde {a} _ {i j} ^ {2} + \tilde {a} _ {j j} ^ {2} = a _ {i j} ^ {2} + a _ {j j} ^ {2}. \tag {5.1}

同理, 由 A^=A~JT\hat{A} = \tilde{A} J^{\mathsf{T}} 可得

a^ii2+a^ij2=a~ii2+a~ij2,a^ji2+a^jj2=a~ji2+a~jj2.(5.2)\hat {a} _ {i i} ^ {2} + \hat {a} _ {i j} ^ {2} = \tilde {a} _ {i i} ^ {2} + \tilde {a} _ {i j} ^ {2}, \quad \hat {a} _ {j i} ^ {2} + \hat {a} _ {j j} ^ {2} = \tilde {a} _ {j i} ^ {2} + \tilde {a} _ {j j} ^ {2}. \tag {5.2}

a^ij=a^ji=0\hat{a}_{ij} = \hat{a}_{ji} = 0 ,故

a^ii2+a^jj2=aii2+ajj2+aij2+aji2=aii2+ajj2+2aij2.\hat {a} _ {i i} ^ {2} + \hat {a} _ {j j} ^ {2} = a _ {i i} ^ {2} + a _ {j j} ^ {2} + a _ {i j} ^ {2} + a _ {j i} ^ {2} = a _ {i i} ^ {2} + a _ {j j} ^ {2} + 2 a _ {i j} ^ {2}.

由于 JAJTJAJ^{\mathsf{T}} 只影响 AA 的第 i,ji,j 行和第 i,ji,j 列, 故对角线元素中只有 aiia_{ii}ajja_{jj} 受影响. 所以

k=1na^kk2=k=1nakk2+2aij2,\sum_ {k = 1} ^ {n} \hat {a} _ {k k} ^ {2} = \sum_ {k = 1} ^ {n} a _ {k k} ^ {2} + 2 a _ {i j} ^ {2},

off(A^)=A^22k=1na^kk2=A22k=1nakk22aij2=off(A)2aij2,\mathrm {o f f} (\hat {A}) = \| \hat {A} \| _ {2} ^ {2} - \sum_ {k = 1} ^ {n} \hat {a} _ {k k} ^ {2} = \| A \| _ {2} ^ {2} - \sum_ {k = 1} ^ {n} a _ {k k} ^ {2} - 2 a _ {i j} ^ {2} = \mathrm {o f f} (A) - 2 a _ {i j} ^ {2},

即引理结论成立.

由此可知, off(A(k))\operatorname{off}(A^{(k)}) 总是不断减小的. 下面给出 Jacobi 迭代算法

算法5.1.Jacobi迭代算法

1: Given a symmetric matrix ARn×nA \in \mathbb{R}^{n \times n}
2: if eigenvectors are desired then
3: set J=IJ = I and flag=1flag = 1
4: end if
5: while not converge do
6: choose an index pair (i,j)(i,j) such that aij0a_{ij} \neq 0
7: τ=(aiiajj)/(2aij)\tau = (a_{ii} - a_{jj}) / (2a_{ij})
8: t=sign(τ)/(τ+1+τ2)t = \operatorname{sign}(\tau) / (|\tau| + \sqrt{1 + \tau^2}) %\% 计算 tanθ\tan \theta
9: c=1/1+t2,s=ctc = 1 / \sqrt{1 + t^2}, \quad s = ct % 计算 cosθ\cos \thetasinθ\sin \theta
10: A=G(i,j,θ)AG(i,j,θ)TA = G(i,j,\theta)AG(i,j,\theta)^{\mathsf{T}}
11: if flag=1flag = 1 then
12: J=G(i,j,θ)JJ = G(i,j,\theta)J

13: end if
14: end while

该算法涉及到 aija_{ij} 的选取问题, 一种直观的选取方法就是使得 aija_{ij} 为所有非对角线元素中绝对值最大的一个, 这就是经典 Jacobi 迭代算法.

算法5.2.经典Jacobi迭代算法

1: Given a symmetric matrix ARn×nA \in \mathbb{R}^{n \times n}
2: if eigenvectors are desired then
3: set J=IJ = I and flag=1flag = 1
4: end if
5: while off(A)>tol\operatorname{off}(A) > \operatorname{tol} do
6: choose (i,j)(i,j) such that aij=maxklakl|a_{ij}| = \max_{k\neq l}|a_{kl}| % 选取绝对值最大的元素
7: τ=(aiiajj)/(2aij)\tau = (a_{ii} - a_{jj}) / (2a_{ij})
8: t=sign(τ)/(τ+1+τ2)t = \mathrm{sign}(\tau) / (|\tau| + \sqrt{1 + \tau^2})
9: c=1/1+t2,s=ctc = 1 / \sqrt{1 + t^2}, \quad s = ct % 计算 cosθ\cos \thetasinθ\sin \theta
10: A=G(i,j,θ)AG(i,j,θ)TA = G(i,j,\theta)AG(i,j,\theta)^{\mathsf{T}}
11: if flag=1flag = 1 then
12: J=G(i,j,θ)JJ = G(i,j,\theta)J
13: end if
14: end while

可以证明,经典Jacobi算法至少是线性收敛的

定理5.3 对于经典Jacobi算法5.2,有

off(A(k+1))(11N)off(A(k)),N=n(n1)2.\mathrm {o f f} (A ^ {(k + 1)}) \leq \left(1 - \frac {1}{N}\right) \mathrm {o f f} (A ^ {(k)}), \quad N = \frac {n (n - 1)}{2}.

kk 步迭代后,有

off(A(k))(11N)koff(A(0))=(11N)koff(A).\operatorname {o f f} (A ^ {(k)}) \leq \left(1 - \frac {1}{N}\right) ^ {k} \operatorname {o f f} (A ^ {(0)}) = \left(1 - \frac {1}{N}\right) ^ {k} \operatorname {o f f} (A).

(板书)

证明. 由于在经典 Jacobi 算法 5.2 中, aij=maxklakl|a_{ij}| = \max_{k \neq l} |a_{kl}| , 故 off(A(k))n(n+1)(aij(k))2\operatorname{off}(A^{(k)}) \leq n(n + 1) \left(a_{ij}^{(k)}\right)^2 , 即

2(aij(k))21Noff(A(k)),N=n(n1)2.2 \left(a _ {i j} ^ {(k)}\right) ^ {2} \geq \frac {1}{N} \mathrm {o f f} (A ^ {(k)}), \quad N = \frac {n (n - 1)}{2}.

所以由引理5.2可知

off(A(k+1))=off(A(k))(aij(k))2(11N)off(A(k)).\operatorname {o f f} (A ^ {(k + 1)}) = \operatorname {o f f} (A ^ {(k)}) - \left(a _ {i j} ^ {(k)}\right) ^ {2} \leq \left(1 - \frac {1}{N}\right) \operatorname {o f f} (A ^ {(k)}).

事实上, 经典 Jacobi 算法最终是 (渐进) 二次收敛的 [30, 100]

定理5.4经典Jacobi算法5.2是 NN 步(渐进)二次收敛的,即对足够大的 kk ,有

off(A(k+N))=O(off2(A(k))).\operatorname {o f f} \left(A ^ {(k + N)}\right) = O \left(\operatorname {o f f} ^ {2} \left(A ^ {(k)}\right)\right).

由于在经典 Jacobi 算法中, 每一步都要寻找绝对值最大的非对角元, 比较费时, 因此实用性较差. 我们可以通过逐行扫描来选取 (i,j)(i,j) , 这就是循环 Jacobi 迭代算法.

算法5.3. 循环Jacobi迭代算法(逐行扫描)
1: Given a symmetric matrix ARn×nA \in \mathbb{R}^{n \times n}
2: if eigenvectors are desired then
3: set J=IJ = I and flag=1flag = 1
4: end if
5: while off (A)>tol(A) > tol do
6: for i=1i = 1 to n1n - 1 do
7: for j=i+1j = i + 1 to nn do
8: if aij0a_{ij} \neq 0 then
9: τ=(aiiajj)/(2aij)\tau = (a_{ii} - a_{jj}) / (2a_{ij})
10: t=sign(τ)/(τ+1+τ2)t = \mathrm{sign}(\tau) / (|\tau| + \sqrt{1 + \tau^2})
11: c=1/1+t2c = 1 / \sqrt{1 + t^2}
12: s=cts = c \cdot t
13: A=G(i,j,θ)TAG(i,j,θ)A = G(i, j, \theta)^{\mathrm{T}} AG(i, j, \theta)
14: if flag = 1 then
15: J=JG(i,j,θ)J = J \cdot G(i, j, \theta)
16: end if
17: end if
18: end for
19: end for
20: end while

循环 Jacobi 也具有 (渐进) 二次收敛性 [135, page 270].

Jacobi迭代法的优缺点

  • 优点:能够达到很高的计算精度 (特别是小特征值); 同时非常适合并行计算.

  • 缺点: 计算速度较慢; 矩阵稀疏性得不到充分的利用.