第三十篇 雅可比主对角线化求对称矩阵的特征值
对于标准特征值方程



由特征值问题编程基础可知，对于任何非0解矩阵[P]，标准方程可以转化为具有相同特征值的方程



其中



这种转换技术的关键核心在于[A *]的特征值比原始[A]的特征值更容易找到。

然而，如果[A]是对称的，则变换后的矩阵不太可能保持对称。很容易证明下面变换



将保持[A∗]的对称性。为了使上面方程中给出的特征值[A∗]与[A]的特征值相同，必须把这两个性质综合起来



这种类型的矩阵称为“正交矩阵”，具有这种性质的矩阵称为“旋转矩阵”。



将此变换应用于下面的矩阵，我们有





其中很明显[A∗]对于任何α值都是对称的。这种情况，明显可以选择一个α值使得[A∗]成为一个对角矩阵，因为如果这样，对角线就是特征值。下面的情况，非对角线项将被消除



得出，tan α = 1和α = π/4，给出sin α = cos α = 1/√2。

得到的变换矩阵是



即[A]的特征值分别为3和1。

对于大于2 × 2的矩阵[A]，变换矩阵[P]必须通过在其他主对角线上放1，在所有非对角线上放0来“填充”，被消去的行和列选择上面的矩阵。例如，如果[A]是4 × 4，则变换矩阵可以选择6种形式中的一种，这取决于要消去初始矩阵中的哪些非对角项，例如:



上面第一个矩阵经过[P]T [A][P]变换后将原矩阵[A]中的a12和a21项消去，而第二个矩阵将消去a24和a42项。1和0的作用是让[A]的其他行和列保持不变。这意味着在一次转换中变为零的非对角线项在随后的转换中会恢复为非零值(尽管通常是很“小”的值)，因此正如期望的那样，该方法是迭代的。

这种类型的迭代的最早形式称为“雅可比对角化”，它通过消除每一次迭代剩余的“最大的”非对角项连续进行迭代。

对于任何对称矩阵[A]，得到广义方程为



得到[A∗]形式的非对角线项为



求α使这一项等于零



因此



因此，为了建立一个雅可比对角化的简单程序，必须在[a]中搜索“最大的”非对角线项，并找到它所在的行和列。“旋转矩阵”α按照之前的方法构建。矩阵[P]可以使用一个numpy库的transpose转化，然后矩阵乘积形成方程的[A *]。重复这个过程，直到[A∗]的主对角线在可接受的公差内收敛到[A]的特征值为止。

计算实例：

使用雅可比对角化去估算下面对称矩阵的特征值



下面的结果保留到小数点后四位，但实际计算的精确度更高。

第一次迭代

最大的非主对角项为a23 = a32 = −9.0,因此根据之前方程



第一次转换矩阵将包含下面的项



因此



通过转化矩阵得到



详细的数值为



最后



第二次迭代

最大的非主对角项为a12 = a21 = −7.7782，因此



第二次转化矩阵将包括下面的项



因此，



同上面一样，矩阵乘积将等于



可以看到，虽然位置(2,3)和(3,2)不再为零，但与初始矩阵中的值相比足够“小”。随着迭代的进行，旋转角度αk→0，变换矩阵[Pk]→[I]和变换矩阵[Ak]趋向于一个对角矩阵，特征值在对角线上。

对于这个例子，经过六次迭代，容差为1.0e-5

得到



因此[A]的特征值λ = 0.4659, 20.9681，−0.9340。特征向量将通过将每个特征值代入求线性方程的解。

程序如下

分为一个主程序和两个子程序，分别为判断收敛的子程序checkit，和高斯消元求特征向量的子程序eliminate。详情可参照之前文章的线性方程求解部分

主程序：
#雅可比主对角线化求对称矩阵的特征值
import numpy as np
import B
import math
n=3;tol=1.0e-5;limit=100
enew=np.zeros((n,1))
eold=np.zeros((n,1))
p=np.zeros((n,n))
a1=np.zeros((n,n))
a=np.array([[10,5,6],[5,20,4],[6,4,30]],dtype=np.float)
a2=a
pi=math.acos(-1)
x=np.zeros((n,1))
x=np.ones((3,1),dtype=np.float)
print('雅可比主对角线化求对称矩阵的特征值')
print('矩阵A')
print(a[:])
print('前几次迭代值')
iters=0;eold[:]=0
while(True):
    iters=iters+1
    big=0
    for i in range(1,n+1):
        for j in range(i+1,n+1):
            if abs(a[i-1,j-1]>big):
                big=abs(a[i-1,j-1]);hold=a[i-1,j-1];nr=i;nc=j
    if abs(big)<1.0e-20:
        break
    den=a[nr-1,nr-1]-a[nc-1,nc-1]
    if abs(den)<1.0e-20:
        alpha=pi/4.0
        if hold<0:
            alpha=-alpha
    else:
        alpha=math.atan(2.0*hold/den)/2.0
    ct=math.cos(alpha);st=math.sin(alpha);p[:]=0
    for i in range(1,n+1):
        p[i-1,i-1]=1.0
    p[nr-1,nr-1]=ct;p[nc-1,nc-1]=ct;p[nr-1,nc-1]=-st;p[nc-1,nr-1]=st
    a=np.dot(np.dot(np.transpose(p),a),p)
    if iters<5:
        for i in range(1,n+1):
            for j in range(1,n+1):
                print('{:13.4e}'.format(a[i-1,j-1]),end='')
            print(end='\n')
        print(end='\n')
    for i in range(1,n+1):
        enew[i-1,0]=a[i-1,i-1]
    if B.checkit(enew,eold,tol) or iters==limit:
        break
    eold[:,0]=enew[:,0]
print('迭代到收敛次数',iters)
print('最后的转化矩阵A')
for i in range(1,n+1):
    for j in range(1,n+1):
        print('{:13.4e}'.format(a[i-1,j-1]),end='')
    print(end='\n')
for i in range(1,n+1):
    a1[:]=a2[:]
    for j in range(1,n+1):
        a1[j-1,j-1]=a1[j-1,j-1]-a[i-1,i-1]
    x[:]=0;a1[i-1,i-1]=1.0e20;x[i-1]=1.0e20;x[:]=B.eliminate(a1,x)
    l2=np.linalg.norm(x)
    print('特征值','{:13.4e}'.format(a[i-1,i-1]))
    print('特征向量')
    for i in range(1,n+1):
        print('{:13.4e}'.format(x[i-1,0]/l2),end=' ')
    print()
    
    


checkit

def checkit(loads,oldlds,tol):
#检查多个未知数的收敛
  neq=loads.shape[0]
  big=0.0
  converged=True
  for i in range(1,neq+1):
    if abs(loads[i-1,0])>big:
      big=abs(loads[i-1,0])
  for i in range(1,neq+1):
    if abs(loads[i-1,0]-oldlds[i-1,0])/big>tol:
      converged=False
  checkit=converged
  return  checkit

eliminate

def eliminate(a,b):
  n=a.shape[0]
##确定主对角线最大值
  for i in range(1,n):
    big=abs(a[i-1,i-1]);ihold=i
    for j in range(i+1,n+1):
      if abs(a[j-1,i-1])>big:
        big=abs(a[j-1,i-1]); ihold=j
    if ihold!=i:
      for j in range(i,n+1):
        hold=a[i-1,j-1]; a[i-1,j-1]=a[ihold-1,j-1]; a[ihold-1,j-1]=hold
      hold=b[i-1,0]; b[i-1,0]=b[ihold-1,0]; b[ihold-1,0]=hold
##消元阶段
    for j in range(i+1,n+1):
      fac=a[j-1,i-1]/a[i-1,i-1]
      for l in range(i,n+1):
        a[j-1,l-1]=a[j-1,l-1]-a[i-1,l-1]*fac
      b[j-1]=b[j-1]-b[i-1]*fac
##从后迭代
  for i in range(n,0,-1):
    hold=0.0
    for l in range(i+1,n+1):
      hold=hold+a[i-1,l-1]*b[l-1]
    b[i-1]=(b[i-1]-hold)/a[i-1,i-1]
  return b

终端输出结果如下：

2阶实对称矩阵特性
定理：2阶实对称矩阵H的特征值是实数

H=[a,b;b,c]     

    a,b,c是实数，λ 是特征值

A=[a-λ,b;b,c-λ]     

    特征值求解方法为：(a- λ )(c- λ) - b2  = 0

    求解方程得到两个根为：
λ=（a+c）±（a+c）2-4（ac-b2）2   

                   （a+c）2-4ac-b2=a-c2+4b2≥0

    所以，在a、b、c为实数时，特征值也是实数。

       

2、特征向量

根据特征值和特征向量的定义：HX=λX，（H-λE）X = 0；因此方程若有解，则

det（H-λE）=0；

设

A=[a-λ,b;b,c-λ]     

则有-b/(a-λ) = (c-λ)/b,   线性齐次方程组AX=0有非零解，其中之一解向量 [1,-b/(a-λ)],归一化后得到标准解。

特征值和特征向量是指对于矩阵A有，Av=lv，v为特征向量，l为特征值。就是求解一个高次方程：det(A-lI)=0
代码如下：
unit Matrix;
interface
uses
  Math, Windows, SysUtils, Variants, Classes;
Type
TSingleExtendedArray =array of extended;
TDoubleExtendedArray=array of array of extended;
TDoubleLongintArray=array of array of longint;
procedure CalculateEigenVV(var EigenLambda: TSingleExtendedArray; var EigenVector: TDoubleExtendedArray; C: TDoubleLongintArray; N: longint; eps: extended);
implementation
procedure CalculateEigenVV(var EigenLambda: TSingleExtendedArray; var EigenVector: TDoubleExtendedArray; C: TDoubleLongintArray; N: longint; eps: extended);
var
  i, j, fi, fj, p, q, TL, k: longint;
  Change: boolean;
  x, y, cn, sn, Omega, fm, Acurrency: Extended;
  Ja: TDoubleExtendedArray;
begin
  setlength(Ja, N, N);
  setlength(EigenVector, N, N);
  setlength(EigenLambda, N);
  for i := 0 to N - 1 do begin
    EigenVector[i, i] := 1.0; EigenLambda[i] := 0;
    for j := 0 to N - 1 do
      if i <> j then EigenVector[i, j] := 0;
  end;
  Acurrency := 0;
  for i := 0 to N - 1 do
    for j := 0 to N - 1 do begin
      Ja[i, j] := C[i, j];
      Acurrency := Acurrency + Ja[i, j] * Ja[i, j];
    end;
  Acurrency := sqrt(2 * Acurrency);
  Change := true;
  repeat
    if not Change then Acurrency := Acurrency * 0.5;
    Change := false;
    for p := 0 to N - 1 do begin
      for q := p + 1 to N - 1 do
        if Abs(Ja[p, q]) > Acurrency then begin
          x := -Ja[p, q]; y := (Ja[q, q] - Ja[p, p]) / 2;
          if (x <> 0) or (y <> 0) then Omega := x / sqrt(x * x + y * y) else Omega := 1;
          if (y < 0.0) then Omega := -Omega;
          sn := 1.0 + sqrt(1.0 - Omega * Omega);
          sn := Omega / sqrt(2.0 * sn);
          cn := sqrt(1.0 - sn * sn);
          fm := Ja[p, p];
          Ja[p, p] := fm * cn * cn + Ja[q, q] * sn * sn + Ja[p, q] * Omega;
          Ja[q, q] := fm * sn * sn + Ja[q, q] * cn * cn - Ja[p, q] * Omega;
          Ja[p, q] := 0; Ja[q, p] := 0;
          for i := 0 to N - 1 do
            if (i <> p) and (i <> q) then begin
              fm := Ja[p, i];
              Ja[p, i] := fm * cn + Ja[q, i] * sn;
              Ja[q, i] := -fm * sn + Ja[q, i] * cn;
            end;
          for i := 0 to N - 1 do
            if (i <> p) and (i <> q) then begin
              fm := Ja[i, p];
              Ja[i, p] := fm * cn + Ja[i, q] * sn;
              Ja[i, q] := -fm * sn + Ja[i, q] * cn;
            end;
          for i := 0 to N - 1 do begin
            fm := EigenVector[p, i];
            EigenVector[p, i] := fm * cn + EigenVector[q, i] * sn;
            EigenVector[q, i] := -fm * sn + EigenVector[q, i] * cn;
          end;
          for i := 0 to N - 1 do
            EigenLambda[i] := Ja[i, i];
          Change := true; break;
        end;
      if Change then break;
    end;
  until Acurrency <= eps;
  setlength(Ja, 0);
end;
end.

第二十九篇 广义特征值Kx=λMx转化成标准对称形式
特征值问题的几种解法要求方程用“标准形式”表示。



其中A是对称的，如果遇到广义特征值方程



它可以转换为像上边的标准形式的方程，但，即使是[K]和[M]是对称的，乘积[K]−1[M]和[M]−1[K]通常也不是对称的。为了让其对称性，可以使用以下的方法推导。

由上式开始，利用Cholesky方法因式分解[M]可得，参见LDLT分解



因此



现在让



或



经过替换之后得到



最后



这是一个标准的特征值方程，左边矩阵[L]−1[K][L]−T也是对称的。

变换后的上面方程与广义方程具有相同的特征值λ，但特征向量是不同的。但一旦求得特征向量{z}，很容易地从求得特征向量{x}。

程序如下：

其中分为一个主程序和两个子程序，分别为乔利斯基分解的子程序ldlt，从前迭代的子程序subfor。
#Kx=λMx转化为一般对称形式
import numpy as np
import B
n=4;tol=1.0e-5;limit=100
d=np.zeros((n,1))
e=np.zeros((n,1))
s=np.zeros((n,n))
c=np.zeros((n,n))
k=np.array([[8,4,-24,0],[4,16,0,4],[-24,0,192,24],[0,4,24,8]],dtype=np.float)
m=np.array([[0.06667,-0.01667,-0.1,0],[-0.01667,0.1333,0,-0.01667],[-0.1,0,4.8,0.1],[0,-0.01667,0.1,0.06667]],dtype=np.float)
x1=np.zeros((n,1))
x=np.ones((4,1),dtype=np.float)
print('矩阵K')
print(k[:])
print('矩阵M')
print(m[:])
print('初始猜测值',x[:,0])
B.ldlt(m,d)
d[:,0]=d[:,0]**0.5
for j in range(1,n+1):
    for i in range(j,n+1):
        m[i-1,j-1]=m[i-1,j-1]/d[j-1,0]
for j in range(1,n+1):
    e[:,0]=k[:,j-1]
    B.subfor(m,e)
    c[:,j-1]=e[:,0]
print('矩阵C')
for i in range(1,n+1):
    for j in range(1,n+1):
        print('{:13.4e}'.format(c[i-1,j-1]),end='')
    print(end='\n')
for j in range(1,n+1):
    e[:,0]=c[j-1,:]
    B.subfor(m,e)
    s[:,j-1]=e[:,0]
print('最后得出的对称矩阵S')
for i in range(1,n+1):
    for j in range(1,n+1):
        print('{:13.4e}'.format(s[i-1,j-1]),end='')
    print(end='\n')



ldlt

def ldlt(a,d):
  n=a.shape[0]
  for k in range(1,n):
    d[0]=a[0,0]
    if abs(a[k-1,k-1])>1.0e-10:
      for i in range(k+1,n+1):
        x=a[i-1,k-1]/a[k-1,k-1]
        for j in range(k+1,n+1):
          a[i-1,j-1]=a[i-1,j-1]-a[k-1,j-1]*x
        d[i-1]=a[i-1,i-1]
    else:
      print('有0向量在第',k,'行')

subfor

def subfor(a,b):
#一个下三角的从前迭代法
  n=a.shape[0]
  for i in range(1,n+1):
    total=b[i-1]
    if i>1:
      for j in range(1,i):
        total=total-a[i-1,j-1]*b[j-1]
    b[i-1]=total/a[i-1,i-1]

终端输出结果如下



获得对称形式后就可以按照之前的方法求得想要获取的特征值，详情请看。

向量迭代法求特征值移位向量迭代求特征值移位取逆迭代求特征值