主要还是调包:
from numpy.linalg import eig
特征值分解: A = P*B*PT 当然也可以写成 A = PT*B*P 其中B为对角元为A的特征值的对角矩阵。
首先A得正定,然后才能在实数域上分解,
>>> A = np.random.randint(-10,10,(4,4)) >>> A array([[ 6, 9, -10, -1], [ 5, 9, 5, -5], [ -8, 7, -4, 4], [ -1, -9, 0, 6]]) >>> C = np.dot(A.T, A) >>> C array([[126, 52, -3, -69], [ 52, 292, -73, -80], [ -3, -73, 141, -31], [-69, -80, -31, 78]]) >>> vals, vecs = eig(C) >>> vals array([357.33597086, 174.10172008, 8.84429957, 96.71800949]) >>> vecs array([[-0.28738314, -0.51589436, -0.38221983, -0.71075449], [-0.87487263, 0.12873861, -0.24968051, 0.39456798], [ 0.2572149 , -0.69304313, -0.33950158, 0.58161018], [ 0.29300052, 0.48679627, -0.82237845, -0.02955945]])
故使用时应先将特征值转换为矩阵:
>>> Lambda = np.diag(vals) >>> Lambda array([[357.33597086, 0. , 0. , 0. ], [ 0. , 174.10172008, 0. , 0. ], [ 0. , 0. , 8.84429957, 0. ], [ 0. , 0. , 0. , 96.71800949]]) >>> np.dot(np.dot(vecs, Lambda), vecs.T) # 与C=A.T*A相等 array([[126., 52., -3., -69.], [ 52., 292., -73., -80.], [ -3., -73., 141., -31.], [-69., -80., -31., 78.]]) >>> np.dot(np.dot(vecs.T, Lambda), vecs) array([[171.65817919, 45.58778569, 53.20435074, 13.37512137], [ 45.58778569, 125.15670964, 28.22684299, 134.91290105], [ 53.20435074, 28.22684299, 129.48789571, 80.5284382 ], [ 13.37512137, 134.91290105, 80.5284382 , 210.69721545]])
故验证了使用np中的eig分解为A=P*B*PT 而不是A=PT*B*P,其中P=vecs,
即 C = vecs * np.diag(vals) * vecs.T # 这里简写*为矩阵乘法
然后再来看使用np中的eig分解出来的vec中行向量是特征向量还是列向量是特征向量,只需验证:A*vecs[0] = vals[0]*vecs[0]
>>> np.dot(C, vecs[0]) array([-12.84806258, -80.82266859, 6.66283128, 17.51094927]) >>> vals[0]*vecs[0] array([-102.69233303, -184.34761071, -136.58089252, -253.97814676]) >>> np.dot(C, vecs[:,0]) array([-102.69233303, -312.62346098, 91.91213634, 104.69962583]) >>> vals[0]*vecs[:, 0] array([-102.69233303, -312.62346098, 91.91213634, 104.69962583])
后者两个是相等的,故使用np中的eig分解出的vecs的列向量是特征向量。
然后我们可以验证P是单位正交矩阵:
>>> np.dot(vecs.T, vecs) array([[ 1.00000000e+00, -7.13175042e-17, -2.45525952e-18, 2.75965773e-16], [-7.13175042e-17, 1.00000000e+00, 2.49530948e-17, -5.58839097e-16], [-2.45525952e-18, 2.49530948e-17, 1.00000000e+00, -7.85564967e-17], [ 2.75965773e-16, -5.58839097e-16, -7.85564967e-17, 1.00000000e+00]]) >>> np.dot(vecs, vecs.T) array([[ 1.00000000e+00, 2.97888865e-16, -2.68317972e-16, 1.69020590e-16], [ 2.97888865e-16, 1.00000000e+00, -4.40952204e-18, -6.24188690e-17], [-2.68317972e-16, -4.40952204e-18, 1.00000000e+00, -1.13726775e-17], [ 1.69020590e-16, -6.24188690e-17, -1.13726775e-17, 1.00000000e+00]]) # 可以看到除对角元外其他都是非常小的数
即PT*P = P*PT = E , PT=P-1。事实上,在求解P的过程中就使用了施密特正交化过程。
另一方面,我们从数学角度来看:
首先补充一些数学知识:
首先AB相似:P-1*A*P=B,AB合同:CT*A*C=B,
二次型:系数在K中的一个n元二次多项式。由其生成的矩阵称为二次型的矩阵,二次型的矩阵一定是对称矩阵!
正定矩阵:实二次型xT*A*x > 0, x为列向量。
性质:假设A为正定矩阵
1、正定矩阵特征值全大于0
2、行列式 |A| >0
3、A合同于单位阵E,即存在可逆方阵C, s.t. CT*E*C = A = CT*C, 显然可得A为对称正定
正交矩阵:A*AT=AT*A=E ,
性质:
1、A的各行/列是单位向量且两两正交
2、AT=A-1
3、|A|=1
4、(Ax,Ay)=(x,y)x,y∈R
酉矩阵:A*AH=AH*A=E 显然为正交矩阵在复数域上的推广。其中H为共轭转置。
性质:
1、A的各行/列是单位向量且两两正交
2、AH=A-1
3、|A|=1
(这里补充一个厄米特矩阵:AH = A)
正规矩阵:A*AH=AH*A (以上的矩阵均有这个性质,故正规矩阵最为广泛)
正规矩阵的充要条件是:存在酉矩阵U,使得A酉相似于对角矩阵B,即UH*A*U=U-1*A*U=B。
A = P*B*P ,其中B为对角元素为A的特征值的对角阵,P的列向量为特征值对应的特征向量(因为B每行乘以P每列)
原文地址:https://www.cnblogs.com/cymwill/p/9937850.html