Jacobi方法

特征值与特征向量的计算

张瑞
中国科学技术大学数学科学学院

Jacobi方法解实对称阵

  • 幂法仅用来求解矩阵的按模最大的特征值。
  • 如果要求解矩阵的所有特征值，通常对矩阵做一系列的正交相似变换， 将矩阵变为上（下）三角阵或对角阵，然后得到矩阵的所有特征值。
  • Jacobi方法用来求解实对称阵的所有特征值。

注意到， 若$A$为$n$阶对称矩阵，则存在正交矩阵$Q$，使得

\[Q^T A Q=\begin{pmatrix} \lambda_1 \\ & \lambda_2 & \\ & & \ddots & \\ & & & \lambda_n \\ \end{pmatrix} \]
  • 但是，直接计算$Q$是非常困难的。
  • 因此，构造一系列特殊形式的正交矩阵$Q_1$, $Q_2$, $\cdots$, 对$A$做正交相似变换，使矩阵的非对角元的比重逐次减少，对角元的比重逐次增大，
  • 直到近似为一个对角阵。

Givens旋转变换

$n\times n$的Givens旋转矩阵为

\[G(i,j,\theta)=\begin{pmatrix} 1 & & & & & & \\ & \ddots & & & & & \\ & & \cos\theta & & \sin\theta & & \\ & & & \ddots & & & \\ & & -\sin\theta & & \cos\theta & & \\ & & & & & \ddots & \\ & & & & & & 1 \end{pmatrix} \]

其中，$\sin\theta$在$(i,j)$处，$-\sin\theta$在$(j,i)$位置。 则$G(i,j,\theta)$为正交阵。

用它对对称矩阵$A=(a_{ij})$作正交相似变换，得到

\[B=(b_{ij})=Q^T(p,q,\theta)AQ(p,q,\theta) \]

易知，矩阵$B$只是$p$行, $q$行与$p$列, $q$列与矩阵$A$的元素不同，且

(1)
\[\begin{cases} b_{ip}=b_{pi}=a_{pi}\cos\theta-a_{qi}\sin\theta, i\neq p,q \\ b_{iq}=b_{qi}=a_{pi}\sin\theta+a_{qi}\cos\theta, i\neq p,q \\ b_{pp}=a_{pp}\cos^2\theta+a_{qq}\sin^2\theta-a_{pq}\sin2\theta \\ b_{qq}=a_{pp}\sin^2\theta+a_{qq}\cos^2\theta+a_{pq}\sin2\theta \\ b_{pq}=b_{qp}=a_{qp}\cos2\theta+\frac{a_{pp}-a_{qq}}2\sin2\theta \end{cases} \]

旋转的目的是降低非对角元的比重，因而取合适的$\theta$使得$b_{pq}=b_{qp}=0$。 即有

\[a_{qp}\cos2\theta+\frac{a_{pp}-a_{qq}}2\sin2\theta=0 \]

因而$\theta$满足，

\[\cot2\theta=-\frac{a_{pp}-a_{qq}}{2a_{pq}}=s \]

由三角函数的万能公式，

\[\cot2\theta=\frac{\cos2\theta}{\sin2\theta}=\frac{\cos^2\theta-\sin^2\theta}{2\sin\theta\cos\theta}=\frac{1-\tan^2\theta}{2\tan\theta} \]

记$t=\tan\theta$，则

  • $t$为方程$t^2+2st-1=0$的根，可以取绝对值较小的那个。
  • $s=0$时，取$t=1$。

因而

\[\begin{cases} \cos\theta=\frac1{\sqrt{1+t^2}} \\ \sin\theta=\frac{t}{\sqrt{1+t^2}} \\ \end{cases} \]

(1)还可以进一步简化为

\[\begin{cases} b_{ip}=b_{pi}=a_{pi}\cos\theta-a_{qi}\sin\theta, i\neq p,q \\ b_{iq}=b_{qi}=a_{pi}\sin\theta+a_{qi}\cos\theta, i\neq p,q \\ b_{pp}=a_{pp} - a_{pq}\tan\theta \\ b_{qq}=a_{qq} + a_{pq}\tan\theta \\ b_{pq}=b_{qp}=0 \end{cases} \]

依据上式，即可得到矩阵$B$。

事实上，

\[\begin{aligned} b_{pp}=&a_{pp}\cos^2\theta+a_{qq}\sin^2\theta-a_{pq}\sin2\theta \\ =&a_{pp}+a_{pp}(\cos^2\theta-1)+a_{qq}\sin^2\theta-a_{pq}\sin2\theta \\ =&a_{pp}+\sin^2\theta(a_{qq}-a_{pp})-a_{pq}\sin2\theta \\ =&a_{pp}+\sin\theta\cos\theta(\tan\theta(a_{qq}-a_{pp})-2a_{pq}) \\ \end{aligned} \]

由

\[\frac{a_{qq}-a_{pp}}{2a_{pq}}=\cot2\theta=\frac{1-\tan^2\theta}{2\tan\theta} \]

得到

\[\tan\theta(a_{qq}-a_{pp})=a_{pq}(1-\tan^2\theta) \]

因而

\[\begin{aligned} b_{pp}=&a_{pp}+\sin\theta\cos\theta(-1-\tan^2\theta)a_{pq} \\ =&a_{pp}-\frac{\sin\theta}{\cos\theta}a_{pq} \\ =&a_{pp}- a_{pq}\tan\theta \end{aligned} \]

定理 1.
Givens正交变换后的矩阵$B$的对角元比重大于矩阵$A$的对角元比重，即

\[\begin{aligned} \sum_{i\neq j}b^2_{ij}=\sum_{i\neq j}a^2_{ij}-2a_{pq}^2 \\ \sum_{i}b^2_{ii}=\sum_{i}a^2_{ii}+2a_{pq}^2 \\ \end{aligned} \]

Jacoib方法

定义 1.
对称矩阵$A$的Jacobi方法是：

  1. 取$A^{(0)}=A$。
  2. 记$A^{(k)}=(a^{(k)}_{ij})$。 选取$p$, $q$满足$|a^{(k)}_{pq}|=\max_{i\neq j} |a^{(k)}_{ij}|$实施Given旋转变换
    \[A^{(k+1)}=(a^{(k+1)}_{ij})=Q^T(p,q,\theta)A^{(k)}Q(p,q,\theta) \]
    构造序列$\{A^{(k)}\}$。
  3. 直到非对角元比重足够小，即
    \[\sum_{i\neq j} \left(a^{(k+1)}_{ij}\right)^2<\epsilon \]

记$S(A^{(k)})=\displaystyle\sum_{i\neq j} (a^{(k)}_{ij})^2$为矩阵$A^{(k)}$ 非对角元的平方和。则有

\[S(A^{(k+1)})=S(A^{(k)})-2(a^{(k)}_{pq})^2 \]

由$|a^{(k)}_{pq}|=\displaystyle \max_{i\neq j}|a^{(k)}_{ij}|$，因而

\[(n^2-n)(a^{(k)}_{pq})^2\geq S(A^{(k)}) \]

进而有

\[S(A^{(k+1)})\leq\left(1-\frac2{n(n-1)}\right)S(A^{(k)}) \]

即，Jacobi方法收敛。

例 1. (例8.8) 用Jacobi方法计算矩阵 $A=\begin{pmatrix}3 & 1 & 2 \\1 & 3 & 4 \\2 & 4 & 6 \\\end{pmatrix}$ 的全部特征值。

解. 取非对角元最大位置$(3,2)$，旋转后有 

 array([[ 3.        , -0.31726406,  2.21344607],
       [-0.31726406,  0.22799813,  0.        ],
       [ 2.21344607,  0.        ,  8.77200187]])

再取位置$(3,1)$，有 

array([[ 2.24892176, -0.30043869,  0.        ],
       [-0.30043869,  0.22799813, -0.10194645],
       [ 0.        , -0.10194645,  9.52308011]])

可以看到Given旋转不能保持0元的位置不动。

因而，Jacoib方法一般用于低阶的“满矩阵”。

6步后 

array([[2.29261064e+00, 0.00000000e+00, 2.41928292e-11],
       [0.00000000e+00, 1.83189762e-01, 3.21856907e-07],
       [2.41928292e-11, 3.21856907e-07, 9.52419960e+00]])

目录

本节读完

例 2.

2.