第三十一篇 豪斯霍德转换成三对角形式

三对角形式

三对角矩阵就是对角线、邻近对角线的上下次对角线上有元素，其他位置均为0的矩阵。

雅可比主对角化求特征值要求转化矩阵应该具有下面的性质

豪斯霍尔德技术要求选择p为

其中{w}是一个归一化的列向量，使其欧几里德范数等于单位，因此

例如，让

按照上面的乘积

然后可以得到

它具有下面的性质

为消去[A]在对角线外第一行的项，向量{w}应被取为

因此，假设[A]为3 × 3，对第一行的变换矩阵为

当乘积[P][A][P]被计算时，结果矩阵的第一行包含以下三项

让

上面方程的项可以写成

通过向量的点乘为点位1，得到

因此通过方程(1)后面两项的平方，再代换上面方程得到

然后可以写到

再取方程（2）的平方，然后把上面两个方程替换进去得到，

然后把向量写成用向量表示的形式

其中

于是，转化矩阵为

对于r的确定，根据之前的等式，应该选择r的符号与a12是相反的。
所以对于通常的第i行，{v}采取下面形式的向量

程序如下

#对称矩阵到三对角矩阵的豪斯霍德推导
import numpy as np
import B
n=4
v=np.zeros((n,1))
p=np.zeros((n,n))
a1=np.zeros((n,n))
a=np.array([[1,-3,-2,1],[-3,10,-3,6],[-2,-3,3,-2],[1,6,-2,1]],dtype=np.float)
print('对称矩阵到三对角矩阵的豪斯霍德推导')
print('系数矩阵A')
print(a[:])
for k in range(1,n-1):
    r=0
    for l in range(k,n):
        r=r+a[k-1,l]*a[k-1,l]
    r=r**0.5
    if r*a[k-1,k]>0:
        r=-r
    h=-1.0/(r*r-r*a[k-1,k])
    v[:]=0
    v[k,0]=a[k-1,k]-r
    for l in range(k+2,n+1):
        v[l-1,0]=a[k-1,l-1]
    p=np.dot(v,np.transpose(v))*h
    for l in range(1,n+1):
        p[l-1,l-1]=p[l-1,l-1]+1.0
    a1=np.dot(p,a)
    a=np.dot(a1,p)
print('转化后的主对角线')
for i in range(1,n+1):
    print('{:13.4e}'.format(a[i-1,i-1]),end='')
print()
print('转化后的非主对角线值')
for i in range(2,n+1):
    print('{:13.4e}'.format(a[i-2,i-1]),end='')

    
  

输出结果如下