线性代数Python计算:矩阵对角化
矩阵对角化。
线性变换
T
T
T的矩阵
A
∈
P
n
×
n
\boldsymbol{A}\in P^{n\times n}
A∈Pn×n的对角化,即寻求对角阵
Λ
\boldsymbol{\Lambda}
Λ,使得
A
\boldsymbol{A}
A~
Λ
\boldsymbol{\Lambda}
Λ,需分几步走:
(1)解方程
det
(
λ
I
−
A
)
=
0
\det(\lambda\boldsymbol{I}-\boldsymbol{A})=0
det(λI−A)=0,得根
λ
1
,
λ
2
,
⋯
,
λ
k
∈
P
\lambda_1,\lambda_2,\cdots,\lambda_k\in P
λ1,λ2,⋯,λk∈P
为
A
\boldsymbol{A}
A的特征值;
(2)对每一个特征值
λ
i
\lambda_i
λi,解齐次线性方程组
(
λ
i
I
−
A
)
x
=
o
(\lambda_i\boldsymbol{I}-\boldsymbol{A})\boldsymbol{x}=\boldsymbol{o}
(λiI−A)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}
(λI−A)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(λI−A)=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} (λ1I−A)=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}
Λ=
−1000−10005
。
写博不易,敬请支持:
如果阅读本文于您有所获,敬请点赞、评论、收藏,谢谢大家的支持!
代码诚可贵,原理价更高。若为AI学,读正版书好。
更多推荐
所有评论(0)