雅可比矩阵 Python numpy

雅可比方法

该方法是求解对称矩阵全部特征值和特征向量的一种方法，它基于以下结论：
①任何实对称矩阵A可以通过正交相似变换成对角型，即存在正交矩阵Q,使得
$Q^TAQ=diag(λ1,λ2,…,λn)$

其中λi(i=1,2,…,n)是A的特征值，Q中各列为相应的特征向量。

②在正交相似变换下,矩阵元素的平方和不变。 即设$A=(a_{ij})_{n×n}$,Q为正交矩阵,记$B=Q^TAQ=(b_{ij})_{n×n}$, 则
$∑_{i,j=1}^n a_{i,j}^2=∑_{i,j=1}^n b_{i,j}^2$

Jacobi方法的基本思想是通过一次正交变换,将A中的一对非零的非对角元素化成零并且使得非对角元素的平方和减小。反复进行上述过程,使变换后的矩阵的非对角元素的平方和趋于零，从而使该矩阵近似为对角矩阵，得到全部特征值和特征向量。

1、矩阵的旋转变换

设A为n阶实对称矩阵，考虑矩阵
$V_{ij}(ϕ) =\left[ \begin{matrix} 1 & 0 & 0 & 0 & 0 \\ 0 & cos(ϕ) & 0 & sin(ϕ) &0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & -sin(ϕ) &0& cos(ϕ) &0\\ 0 & 0 & 0 & 0 & 1 \end{matrix} \right]$
易见$V_{ij}(ϕ)$是正交矩阵，记$A^{(1)}=V_{ij}AV_{ij}^T = (a_{ij}^{(1)})$，可得
$\left\{ \begin{array}{c} a_{ii}^{(1)}=a_{ii} cos^2(ϕ)+a_{jj} sin^2(ϕ) + a_{ij}sin(2ϕ)\\ a_{jj}^{(1)}=a_{ii} sin^2(ϕ)+a_{jj} cos^2(ϕ) - a_{ij}sin(2ϕ)\\ a_{ik}^{(1)}=a_{ki}^{(1)} = a_{ik}cos(ϕ)+a_{jk} sin(ϕ) \\ a_{kj}^{(1)}=a_{jk}^{(1)}=-a_{ki} sin(ϕ)+a_{kj} cos(ϕ) \\ a_{kl}^{(1)}=a_{lk}^{(1)}=a_{kl}\\ a_{ij}^{(1)}=a_{ji}^{(1)} =(ajj - aii) sin(ϕ)cos(ϕ)+a_{ij}cos(2ϕ)\\ \end{array}\right.$
若$a_{ij}≠0$，取ϕ使得 $tan(2ϕ) = \frac{2a_{ij}}{(a_{ii}-a_{jj})}$，其中$|ϕ|<\frac{\pi}{4}$，则有
$a_{ij}^{(1)} = a_{ji}^{(1)} =0$

推导如下

注意到$B=V_{ij}A$的第i,j行元素 以及$BV_{ij}^T = A$的第i,j行元素 为
$\left\{ \begin{array}{c} b_{ik}=a_{ik} cos(ϕ)+a_{jk} sin(ϕ)\\ b_{jk}=-a_{ik} sin(ϕ)+a_{jk} cos(ϕ)\\ \\ a_{ki}^{(1)}=b_{ki} cos(ϕ)+b_{kj} sin(ϕ) \\ a_{kj}^{(1)}=-b_{ki} sin(ϕ)+b_{kj} cos(ϕ) \end{array}\right.$

分别令$k=i，j$，可求得$a_{ii}^{(1)},a_{ij}^{(1)},a_{jj}^{(1)},a_{ji}^{(1)}$，其中有
$\left\{ \begin{array}{c} a_{ii}^{(1)}=b_{ii} cos(ϕ)+b_{ij} sin(ϕ) \\ a_{jj}^{(1)}=-b_{ji} sin(ϕ)+b_{jj} cos(ϕ) \end{array}\right.$
将$b$带入$a^{(1)}$,
$\left\{ \begin{array}{c} a_{ii}^{(1)}=(a_{ii} cos(ϕ)+a_{ji} sin(ϕ))\cdot cos(ϕ)+(a_{ij} cos(ϕ)+a_{jj} sin(ϕ) )\cdot sin(ϕ) \\ a_{jj}^{(1)}=-(-a_{ii} sin(ϕ)+a_{ji} cos(ϕ))\cdot sin(ϕ)+(-a_{ij} sin(ϕ)+a_{jj} cos(ϕ))\cdot cos(ϕ) \end{array}\right.$
基于对称矩阵具有$a_{ji} = a_{ij}$，化简可得
$\left\{ \begin{array}{c} a_{ii}^{(1)}=a_{ii} cos^2(ϕ)+a_{jj} sin^2(ϕ) + a_{ij}sin(2ϕ)\\ a_{jj}^{(1)}=a_{ii} sin^2(ϕ)+a_{jj} cos^2(ϕ) - a_{ij}sin(2ϕ)\\ \end{array}\right. \tag{1}$
基于对称性，有$a_{ij}^{(1)}=a_{ji}^{(1)}$，可求得：
$a_{ji}^{(1)}=a_{ij}^{(1)}=-b_{ii} sin(ϕ)+b_{ij} cos(ϕ)$
.
将$b_{ik}$带入，可得
$a_{ji}^{(1)}=-(a_{ii} cos(ϕ)+a_{ji} sin(ϕ))\cdot sin(ϕ)+ (a_{ij} cos(ϕ)+a_{jj} sin(ϕ))\cdot cos(ϕ)$
.
基于对称矩阵具有$a_{ji} = a_{ij}$，化简可得
$a_{ij}^{(1)}=a_{ji}^{(1)} =(ajj - aii) sin(ϕ)cos(ϕ)+a_{ij}cos(2ϕ) \tag{2}$
.
当$a_{ij}≠0$，取ϕ使得 $tan(2ϕ) = \frac{2a_{ij}}{(a_{ii}-a_{jj})}$，其中$|2ϕ|<\frac{\pi}{2}$，则
$tan(2ϕ) = \frac{c \cdot sin(2ϕ)}{c \cdot cos(2ϕ)}$ ,可得
$\left\{ \begin{array}{c} 2a_{ij}=c \cdot sin(2ϕ)\\ a_{ii}-a_{jj}=c \cdot cos(2ϕ)\\ \end{array}\right.$
则有$a_{ij}^{(1)} = a_{ji}^{(1)} =c \cdot cos(2ϕ) \cdot (0.5\cdot sin(2ϕ) )- 0.5c \cdot sin(2ϕ)cos(2ϕ)$，化简可得
$a_{ij}^{(1)} = a_{ji}^{(1)} = 0 \tag{3}$

当$k≠i,j$时，仅有$v_{kk}=1$，故有
$\left\{ \begin{array}{c} b_{ki}=v_{ii} \cdot a_{ik} + 0 = a_{ik}\\ b_{kj}=v_{jj} \cdot a_{jk} + 0 = a_{jk}\\ \end{array}\right.$
将$b$带入$a^{(1)}$,且基于矩阵对称性，可得
$\left\{ \begin{array}{c} a_{ik}^{(1)}=a_{ki}^{(1)} = a_{ik}cos(ϕ)+a_{jk} sin(ϕ) \\ a_{kj}^{(1)}=a_{jk}^{(1)}=-a_{ki} sin(ϕ)+a_{kj} cos(ϕ) \\ \end{array}\right. \tag{4}$
且当 $l≠i,j$时，$V^T$中仅有$v_{kk}=1$，故有$a_{kl}^{(1)}= v_{ll} \cdot a_{kl} + 0$ , 可得
$a_{kl}^{(1)}=a_{lk}^{(1)}=a_{kl} \tag{5}$

对 $A^{(1)}$重复上述过程，得到$A^{(2)}$…，可以证明输入这种变换不一定使得矩阵中非对角元素中零元素的个数单调增加，但可以保证非对角元素的平方和递减。

以A和$A^{(1)}$为例，设$E(A)=∑_{k≠l}a_{kl}^2$，$E(A)^{(1)}=∑_{k≠l}（a_{kl}^{(1)}）^2$，即
.
$E(A)^{(1)}=∑_{k≠l 且 k,l ≠ i,j}（a_{kl}^{(1)}）^2 + 2∑_{k≠l 且 k,l = i,j}(（a_{ik}^{(1)}）^2 + （a_{jk}^{(1)}）^2)$
由 (5）、（3）二式可得
$E(A)^{(1)}= ∑_{k≠l 且 k,l ≠ i,j}a_{kl}^2 + 0$
而$E(A)$有
$E(A)=∑_{k≠l 且 k,l ≠ i,j}a_{kl}^2 + 2∑_{k≠l 且 k,l = i,j}(a_{ik}^2 + a_{jk}^2)\\ \ \\ = ∑_{k≠l 且 k,l ≠ i,j}a_{kl}^2 + 2a_{ij}^2$
由于（3）式条件为$a_{ij}≠0$，故在上述旋转变换下， $E(A)^{(1)} = E(A) - 2a_{ij}^2 ＜ E(A)$，非对角元素平方和严格单调递减，由(1)式可知，对角元素平方和单调增加，当然，由旋转特性（旋转后向量长度不变）也可以得到结论。

2、Jacobi方法

通过一系列旋转变换将$A$变成$A^{(k+1)}$，求得A的全部特征值和特征向量，该方法被称为Jacobi方法。计算过程如下

①令k=0，$A^{(k)}=A$
②求整数$i,j$，使得$\{a_{i,j}^{(k)}\}=max |a_{l,m}^{(k)}|$，其中$l≠m$
③计算旋转矩阵
$o = cot(2ϕ) = \frac{(a_{ii}^{(k)}-a_{jj}^{(k)})}{2a_{ij}^{(k)}}$
$t = tan(ϕ) = sign(o) \cdot (\sqrt{o^2+1} - |o|)$
$c = cos(ϕ) =\frac{1}{\sqrt{1+t^2}}$
$s = sin(ϕ) = t \cdot c$
$V^{(k)} = V_{i,j}(ϕ)$

④计算$A^{(k)}$

$a_{l,m}^{(k+1)}=a_{m,l}^{(k+1)} = a_{lm}^{(k)}$，其中$l,m ≠ i,j$
$a_{i,j}^{(k+1)}=a_{j,i}^{(k+1)} = 0$

$a_{i,i}^{(k+1)}=c^2 \cdot a_{i,i}^{(k)} + s^2 \cdot a_{j,j}^{(k)}+ 2s\cdot c\cdot a_{i,j}^{(k)}$
$a_{j,j}^{(k+1)}= s^2 \cdot a_{i,i}^{(k)}+c^2 \cdot a_{j,j}^{(k)} - 2s\cdot c\cdot a_{i,j}^{(k)}$

再计算$a_{i,l}^{(k+1)}、a_{j,l}^{(k+1)}$，其中$l ≠ i,j$

$a_{i,l}^{(k+1)}= a_{l,i}^{(k+1)} = c\cdot a_{i,l}^{(k)}+s\cdot a_{j,l}^{(k)}$
$a_{j,l}^{(k+1)}= a_{l,j}^{(k+1)} = -s\cdot a_{i,l}^{(k)}+c\cdot a_{j,l}^{(k)}$

⑤计算$E(A)^{(k+1)}$

$E(A)^{(k+11)}=∑_{l≠m}（a_{l,m}^{(k+1)}）^2$

⑥若$E(A)^{(k+11)} < ε$，则退出。可得特征值和特征向量

补充一点，在步骤③中，运用三角函数知识可由$cot(2ϕ)$求得$tan(ϕ)$
$cot(2ϕ)=\frac{1}{tan (2ϕ)} \\ \ \\ ∵tan(2ϕ)=\frac{2 tan(ϕ)}{1-tan^2(ϕ)}\\ \ \\ ∴cot(2ϕ)=\frac{1-tan^2(ϕ)}{2 tan(ϕ)}$
.
令$o=cot(2ϕ)$,$t=tan(ϕ)$, 则有$t^2 +2 o\cdot t - 1=0$
.
运用二次求根公式 $x=\frac{-b±\sqrt{b^2 - 4ac}}{2a}$，当$\Delta=b^2 - 4ac>0时$可得
.
$t=\frac{-2o ±\sqrt{4o^2 + 4}}{2}=-o ±\sqrt{o^2 + 1}$
.
由于$\sqrt{o^2 + 1}>|o|$，故$t=-o ± |o+\delta|$，在区间$[-\pi/4,\pi/4]$上,

假设$ϕ>0$，$t=-o - |o+\delta|<0$与其$tan(ϕ)$值矛盾，故取$t=-o+\sqrt{o^2 + 1}$

假设$ϕ<0$，$t=-o + |o+\delta|>0$与其$tan(ϕ)$值矛盾，故取$t=-o - \sqrt{o^2 + 1}$

综上，有：
$o>0$时，取$t=-o + \sqrt{o^2 + 1}$，等价于$t=-|o| + \sqrt{o^2 + 1}$
$o<0$时，取$t=-o - \sqrt{o^2 + 1}$，等价于$t=-(-|o|) - \sqrt{o^2 + 1}$
设sign为符号函数，则有
$t=sign(o)\cdot(-|o| + \sqrt{o^2 + 1})$
.
.
若ϕ很小，则o就很大，可能导致$o^2$在计算机上溢出，可取极限分析
$t=lim_{o->+∞} \sqrt{o^2 +1}-o$

采用$\sqrt{o^2+1}+o$进行分子有理化，得到
$\frac{(\sqrt{o^2 +1}-o)\cdot (\sqrt{o^2 +1}+o)}{\sqrt{o^2 +1}+o}=\frac{1}{\sqrt{o^2 +1}+o}$

故当ϕ很小时，有$tan(ϕ)=lim_{o->+∞} \sqrt{o^2 +1}-o = \frac{1}{2o}$
.
由$tan(ϕ)=\frac{sin(ϕ)}{cos(ϕ)}$，可得$sin(ϕ)=tan(ϕ)\cdot cos(ϕ)$
.
由$1+tan²ϕ=\frac{1}{cos²ϕ}$三角函数关系，可推得$cos(ϕ)=\frac{1}{\sqrt{tan^2(ϕ)+1}}$

更多参考雅可比矩阵（Jacobi）