线性代数Python计算:矩阵对角化

在这里插入图片描述
线性变换 T T T的矩阵 A ∈ P n × n \boldsymbol{A}\in P^{n\times n} APn×n的对角化,即寻求对角阵 Λ \boldsymbol{\Lambda} Λ,使得 A \boldsymbol{A} A~ Λ \boldsymbol{\Lambda} Λ,需分几步走:
(1)解方程 det ⁡ ( λ I − A ) = 0 \det(\lambda\boldsymbol{I}-\boldsymbol{A})=0 det(λIA)=0,得根
λ 1 , λ 2 , ⋯   , λ k ∈ P \lambda_1,\lambda_2,\cdots,\lambda_k\in P λ1,λ2,,λkP
A \boldsymbol{A} A的特征值;
(2)对每一个特征值 λ i \lambda_i λi,解齐次线性方程组 ( λ i I − A ) x = o (\lambda_i\boldsymbol{I}-\boldsymbol{A})\boldsymbol{x}=\boldsymbol{o} (λiIA)x=o,得基础解系 α i 1 , α i 2 , ⋯   , α i n i \boldsymbol{\alpha}_{i1},\boldsymbol{\alpha}_{i2},\cdots,\boldsymbol{\alpha}_{in_i} αi1,αi2,,αini i = 1 , 2 , ⋯   , k i=1,2,\cdots,k i=1,2,,k
(3)若 n 1 + n 2 + ⋯ + n k = n n_1+n_2+\cdots+n_k=n n1+n2++nk=n,则
Λ = diag ( λ 1 , ⋯   , λ 1 ⏟ n 1 , ⋯   , λ k , ⋯   , λ k ⏟ n k ) \boldsymbol{\Lambda}=\text{diag}(\underbrace{\lambda_1,\cdots,\lambda_1}_{n_1},\cdots,\underbrace{\lambda_k,\cdots,\lambda_k}_{n_k}) Λ=diag(n1 λ1,,λ1,,nk λk,,λk)
A \boldsymbol{A} A~ Λ \boldsymbol{\Lambda} Λ。即 T T T在基 α 11 , ⋯   , α 1 n 1 , ⋯   , α k 1 , ⋯   , α k n 1 \boldsymbol{\alpha}_{11},\cdots,\boldsymbol{\alpha}_{1n_1},\cdots,\boldsymbol{\alpha}_{k1},\cdots,\boldsymbol{\alpha}_{kn_1} α11,,α1n1,,αk1,,αkn1下的矩阵为 Λ \boldsymbol{\Lambda} Λ
numpy.linalg提供了函数eigvals用来计算方阵的特征值,其调用接口为
eigvals(A) \text{eigvals(A)} eigvals(A)
参数A表示方阵 A \boldsymbol{A} A。返回值为 A \boldsymbol{A} A n n n个根(包括重根)。需要提起注意的是,此处返回的根有可能是复数。对函数eigvals算出 A \boldsymbol{A} A的每个ℝ中的特征值 λ \lambda λ,调用博文《线性方程组的通解》中定义的mySolve函数,计算 ( λ I − A ) x = o (\lambda\boldsymbol{I}-\boldsymbol{A})\boldsymbol{x}=\boldsymbol{o} (λIA)x=o的基础解系,记为属于 λ \lambda λ的线性无关的特征向量。把属于各不同特征值的特征向量按序排列,若这些向量的总数恰为 n n n,则构成向量空间的一个基底,且各特征值(包括重复值)构成的对角阵就是 T T T在这个基下的矩阵。
例1 用Python计算矩阵 A = ( 1 2 2 2 1 2 2 2 1 ) ∈ \boldsymbol{A}=\begin{pmatrix}1&2&2\\2&1&2\\2&2&1\end{pmatrix}\in A= 122212221 3 × 3 ^{3\times3} 3×3的对角化。

import numpy as np                      			#导入numpy
from fractions import Fraction as F                 #导入Fraction
np.set_printoptions(formatter=                      #设置输出数据格式
                    {'all':lambda x:str(F(x).limit_denominator())})
A=np.array([[1,2,2],                    			#设置矩阵
            [2,1,2],
            [2,2,1]],dtype='float')
n,_=A.shape                             			#A的阶数
w=np.linalg.eigvals(A)                  			#A的特征值
Lambda=np.diag(np.sort(w))              			#对角阵
v=np.unique(np.round(w,decimals=14))    			#特征值唯一化
I=np.eye(n)                             			#单位阵
o=np.zeros((n,1))                       			#零向量
lam=v[0]                                			#第1个特征值
P=(mySolve(lam*I-A,o))[:,1:]            			#属于第1个特征值的特征向量
for lam in v[1:]:                       			#其余个特征值
    X=mySolve(lam*I-A,o)                			#属于特征值的特征向量
    P=np.hstack((P,X[:,1:n]))
print('对角阵:')
print(Lambda)
print('基:')
print(P)

程序的第5~7行设置矩阵 A \boldsymbol{A} A的数据为A。第8行读取 A \boldsymbol{A} A的阶数为n。第9行调用numpy.linalg的eigvals函数计算 A \boldsymbol{A} A n n n个特征值存于w。第10行调用numpy的diag函数用w按升序排列的 n n n个值构造对角阵Lambda。第11行调用numpy的unique函数将w中的 n n n个特征值删掉重复,仅保留不同值。即 det ⁡ ( λ I − A ) = 0 \det(\lambda\boldsymbol{I}-\boldsymbol{A})=0 det(λIA)=0的所有重根仅算一个。得到的各不相同的特征值按升序存于v。注意,w中存储的 A \boldsymbol{A} A的特征值都是浮点型的,为能确定重根,需限制精度。np.round(w,decimals=14)调用numpy的round函数将w中的值限制有效位数为14。第12行设置 n n n阶单位阵I,第13行设置 n n n-维零向量o。第14行读取第1个特征值 λ 1 \lambda_1 λ1为lam,第15行调用mySolve函数(见博文《线性方程组的通解》)解齐次线性方程组 ( λ 1 I − A ) = o (\lambda_1\boldsymbol{I}-\boldsymbol{A})=\boldsymbol{o} (λ1IA)=o,将所得基础解系(存于返回值的第2列起的各列)置于P中。第16~18行的for循环将其余各特征值所属特征向量依次置于P后,最终得到对角化后的基底。运行程序,输出

对角阵:
[[-1  0 0]
 [ 0 -1 0]
 [ 0  0 5]]
基:
[[-1 -1 1]
 [ 1  0 1]
 [ 0  1 1]]

A = ( 1 2 2 2 1 2 2 2 1 ) ∈ \boldsymbol{A}=\begin{pmatrix}1&2&2\\2&1&2\\2&2&1\end{pmatrix}\in A= 122212221 3 × 3 ^{3\times3} 3×3在基 α 1 = ( − 1 1 0 ) , α 2 = ( − 1 0 1 ) \boldsymbol{\alpha}_1=\begin{pmatrix}-1\\1\\0\end{pmatrix},\boldsymbol{\alpha}_2=\begin{pmatrix}-1\\0\\1\end{pmatrix} α1= 110 ,α2= 101 α 3 = ( 1 1 1 ) \boldsymbol{\alpha}_3=\begin{pmatrix}1\\1\\1\end{pmatrix} α3= 111 下对角化为 Λ = ( − 1 0 0 0 − 1 0 0 0 5 ) \boldsymbol{\Lambda}=\begin{pmatrix}-1&0&0\\0&-1&0\\0&0&5\end{pmatrix} Λ= 100010005
写博不易,敬请支持:
如果阅读本文于您有所获,敬请点赞、评论、收藏,谢谢大家的支持!
代码诚可贵,原理价更高。若为AI学,读正版书好

评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值