matlab矩阵的分解

1、第三章 矩阵的分解(一) 矩阵的特征值与特征向量(Eigenvalues and EigenVectors) 1. 矩阵的特征值与特征向量 解Ax=x 运算式中的及其所对应的非零的向量x , 我们称/ x 为矩阵A的 特征值与特征向量。 改写原式为 , (A-I) x = 0 , I 是单位矩阵, 我们令P() = det(A-I) = 0, 则 P()的展开式称为矩阵 A 的特征多项式, 解出矩阵 A 的特征多项式 , 就可得矩阵 A 的所有eigenvalues 。再将每一个 eigenvalue 代入原式中, 即可求出其相对应的 eigenvectors 。 例 1 : 解矩阵A = -

2、9 -3 -16 ; 13 7 16; 3 3 10 的特征值与特征向量。 【解 1】 先利用函数 poly() 求出矩阵 A 的特征多项式, 再用roots()函数 , 求出特征多项式所有的根。 A= -9 -3 -16; 13 7 16; 3 3 10 ; poly(A) %利用一个向量来储存此多项式的系数 roots(poly(A) ans = 1.0000 -8.0000 -44.0000 240.0000ans = 10.0000 4.0000 -6.0000 上面输出结果中, 第一个 ans 是 A 的特征多项式的系数, 即 第二个 ans 是 A 的eigenvalues : 1

3、0, 4, -6 接着针对某个特征值 , 我们找出其对应之特征向量 利用 rref() 函数, 求出 (A-I) 的 row reduced echelon form或是 利用 null() 函数, 求出 (A-I) null space 的基底向量 A = -9 -3 -16; 13 7 16; 3 3 10 ; rref(A - 10*eye(size(A) null(A - 10*eye(size(A) ans = 1 0 1 0 1 -1 0 0 0ans = 0.5774 -0.5774 -0.5774 上面输出结果中, 第一个 ans 是的 reduced row echelon

4、form 即 令 , 得 为 10所对应的eigenvectors 第二个 ans 是 null space 的基底向量, 这个基底向量 的长度为1. 上述的解x, 当取 t=-1再除以norm(x), 即可得这个 基底向量。 依此方法, 将其他的 eigenvectors 求出。 【解 2】 使用 matlab 函数 eig(), eig() 有两种不同的输出形式 : eig(A) 只传回eigenvalues, 而 V,D = eig(A) 则传回D是由 eigenvalues形成的diagonal matrix; V是由eigenvalues所对应的 eigenvectors形成的mat

5、rix , 满足A*V=V*D , 当 V是nonsingular 时, 则表示矩阵A可对角化。 A = 1 1 -1; 2 0 1; 1 1 0 ; eig(A) % 传回eigenvalues所形成的行向量 V, D = eig(A) % D是由eigenvalues形成的diagonal matrix % V是由eigenvalues所对应的eigenvectors形成的matrix ans = 1.4142 -1.4142 1.0000V = -0.1913 -0.4619 0.0000 -0.7325 0.8446 0.7071 -0.6533 -0.2706 0.7071D = 1

6、.4142 0 0 0 -1.4142 0 0 0 1.0000 在上面的输出结果中, V的第一行行向量, 为eigenvalue 1.4142 所对应的eigenvector ; 第二行行向量, 为eigenvalue -1.4142所 对应的eigenvector ; 第三行行向量, 即为eigenvalue 1.000所 对应的eigenvector。 rank(V) %检查V的行向量是否线性独立 inv(V)*A*V %验证A是否对角化 ans = 3ans = 1.4142 -0.0000 -0.0000 0.0000 -1.4142 -0.0000 -0.0000 -0.0000

7、1.0000 从 rank(V) =3, V的行向量是线性独立, 得知A可对角化 而且V-1AV也几乎近似于diagonal matrix D . format long %观察数值的准确度 inv(V)*A*V ans = 1.41421356237309 -0.00000000000000 -0.00000000000000 0.00000000000000 -1.41421356237310 -0.00000000000000 -0.00000000000000 -0.00000000000000 1.00000000000000 2. 矩阵的对角化 Theorem : An n-by-

8、n matrix A is diagonalizable if and only if it has n linearly independent eigenvectors . 也就是说 ; 如果 A 有 n 个线性独立的 eigenvectors , 那么矩阵A 就是可以对角化的。此外 , 我们也可证明当特征值不同时, 其 所对应的特征向量是线性独立的。从上面的例子, 我们可观察 到这个结论。 例 : 求出矩阵 A=-28 13 -42; -10 13 -6; 19 -7 30 的eigenvalues & eigenvectors, 判断矩阵 A 是否可以对角化 ? 【解】 for

9、mat short A=-28 13 -42; -10 13 -6; 19 -7 30 V,D=eig(A) A = -28 13 -42 -10 13 -6 19 -7 30V = -0.7570 -0.6667 -0.6667 -0.5179 -0.6667 -0.6667 0.3984 0.3333 0.3333D = 3.0000 0 0 0 6.0000 0 0 0 6.0000 从上面D的结果eigenvalues are 3, 6, and 6 ; 6 的重根数为2, 而其 对应的eigenvectors, V的第二行行向量与第三行行向量, 看起来似 乎是一样的, 并不线性独立,

10、 故矩阵 A不可以对角化。 若我们观察V的秩数(rank), 会发现一个问题: rank(V) ans = 3 从ans =3 , 表示V的行向量是线性独立, 与上述看到的结果有异。 下面我们使用format long的格式注意其数值的变化 . format long V D format short %改回原来的格式 V = -0.75697811924511 -0.66666669551571 -0.66666663781762 -0.51793239737824 -0.66666662470442 -0.66666670862891 0.39840953644480 0.33333335

11、955974 0.33333330710693D = 3.00000000000014 0 0 0 5.99999936269826 0 0 0 6.00000063730158 因为计算误差的结果Matlab把V的第二, 第三行向量, 视为不同 而且线性独立 . 因此我们得到rank(V) =3 . 让我们真正来计算6对应的eigenvectors. V1=rref(A- 6*eye(3) %rref of matrix A-6I v1=null(A- 6*eye(3) V1 = 1 0 2 0 1 2 0 0 0v1 = 0.6667 0.6667 -0.3333 从上面的结果EigenV

12、alue 6的几何重根数(Geometric Multiplicity) = 1, 表示其特征向量与零向量形成的向量空间(eigen space)之维度只有1, 而代数重根数(Algebraic Multiplicity) = 2. 因此, EigenValue 6必然有 一个 2*2的Jordan block, 下面我们将说明矩阵的Jordan form . 3. 矩阵的Jordan Form 8, p. 358 所谓Jordan block J()是一个上三角方阵, 其主对角上的元素是, 主对角上面 的元素是 1, 其余则全为0 . 即 1 0 0 J() = 0 1 0 . . . 0

13、0 0 Theorem : Each p*p matrix A is similar to a matrix J in Jordan form , J=Q-1AQ . with J1 0 0 Q-1AQ = 0 J2 0 0 0 Jk where each Jr is an nr*nr Jordan block and k= k1+ k2+ ks equals the sum of the geometric multiplicities of the distinct eigenvalues1, 2, s of A. The same eigenvalue may occur in diff

14、erent Jordan blocks Jr , but the total number of blocks with that eigenvalue equals its geometric multiplicity kj , j = 1, 2 , , s . 也就是说 ; 针对某个eigenvalue可能出现在不同的Jordan Block, 但 其Block的总个数, 等于它的几何重根数 . 而在所有 Block主对角上, 此eigenvalue出现的总数, 就是它的 代数重根数 . 考虑nr*nr Jordan block Jr , 其作用在矩阵Q 的行向量, 若令有影响到的行向量

15、为 Vr1, Vr2, , Vrnr ; 从 AQ=QJ 我们可得到下列关系式 : A Vr1 = r Vr1 and A Vrj = rVrj + Vr,j-1 for j=2, , nr 其中Vr1是r的eigenvector, 而Vrj 则称为generalized eigenvector . 从上个例子中, 我们的第一个eigenvalue1=3 所对应的Jordan block, J(3)=3 而eigenvector为V11=V(:,1)=-0.7570 -0.5179 0.3984 , 第二个2=6 所对应的 Jordan block , J(6)=6 1 ; 0 6 , 因为它

16、的几何重根数=1; 而代数重根数 =2 . 若我们要找出矩阵Q 的行向量, 则必须找到一个generalized eigenvector, for 2=6 . A1=A-6*eye(3); %解(A-6I) x = b, b=2 2 -1' b=2 2 -1' %b is an eigenvector of A with %eigenvalue=6 matgv=rref(A1 b) matgv = 1.0000 0 2.0000 0.1111 0 1.0000 2.0000 0.4444 0 0 0 0 从matgv的结果得 (A-6I)x=b=2 2 -1'的解, g

17、eneralized eigenvectors, 为 t*b+0.1111 0.4444 0' 若取t=0, 则得到向量gv=0.1111 0.4444 0' 令 Q=V(:,1) V(:,2) gv ; 我们计算Q-1AQ的结果 gv=matgv(:,4) ; Q=V(:,1) b gv ; Jd=inv(Q)*A*Q Jd = 3.0000 -0.0000 0 -0.0000 6.0000 1.0000 0.0000 0 6.0000 上述的结果说明, 虽然矩阵A不可对角化, 但是可相似(similar to) 于一个矩阵Jd, in Jordan Form . 例 : 找

18、一个矩阵A=2 2 -1 ; -1 -1 1; -1 -2 2 的 Jordan Form. 【解】 先找出A的eigenvalues以及其对应的几何重根数与代数重根数 . 再找出矩阵Q 的行向量, 即eigenvectors / generalized eigenvectors. A=2 2 -1 ; -1 -1 1; -1 -2 2 ; eig(A) ans = 1.0000 + 0.0000i 1.0000 - 0.0000i 1.0000 为了解上面ans复数的部分, 是否为误差值 ? 我们观察其特征多项式 poly(A) %检视上面ans复数的部分是否为误差值 ans = 1.000

19、0 -3.0000 3.0000 -1.0000 从ans的结果, 我们得到矩阵 A的特征多项式P()=( -1 )3 , 故其 特征值=1 有代数重根数 3; rref(A-eye(3) ans = 1 2 -1 0 0 0 0 0 0 我们得到矩阵 A的特征向量x, 可写成 x= -2+ = 2 -1 0 ' + 1 0 1' 2 -1 0 ' , 1 0 1' 线性独立, 所以特征值 =1有几何重根数 2 ; 因此有两个Jordan blocks, 一个是1*1而另一个是 2*2 . 令矩阵Q的第一个行向量q1= 1 0 1' , 第二行向量 q2

20、 =2 -1 0 ' +1 0 1' , 则必须找到一个generalized eigenvector当 做第三行向量q3 , 所以解 (A-I)x = q2 =2+ - ' , 其解为q3 因为 1 2 -1 (2+) rref(A-I q2)= 0 0 0 + 0 0 0 (2+2) 若要有解, 则势必= -; 取=1 , 得 = -1 因此 q2 =1 -1 -1' , A1= A-eye(3) ; q2 =1 -1 -1' ; rref(A1 q2) ans = 1 2 -1 1 0 0 0 0 0 0 0 0 我们得q3 = 1 0 0 

21、9; and Q = q1 q2 q3 , 检验Q-1AQ 是否是 Jordan Form ? q1 = 1 0 1' ; q2 = 1 -1 -1' ; q3 = 1 0 0' ; Q =q1 q2 q3 ; inv(Q)*A*Q ans = 1 0 0 0 1 1 0 0 1 从上面的例子我们可见, 若要找出矩阵Q使得 Q-1AQ=J, 并不容易; 尤 其当矩阵的维度增大时, 其特征值的代数重根数与几何重根数, 也很可 能增大 ; 因此要判断Jordan Blocks的形式倍加困难 . 例如下列矩阵A 2 -1 0 0 0 2 0 0 A= 0 0 2 1 0 0

22、0 2 同学可自行检验, A 的特征值只有2, 其代数重根数是4而几何重根数 为2, 所以有两个 Jordan Blocks, 但是无法确定Block的大小是1*1与 3*3, 还是两个2*2 . 下面的定理5, p. 125可供我们判断 . Theorem: Let A be a n*n matrix andbe a eigenvalue of A. If k is the number of k*k blocks J(), andk = dim ker(A-I)k , then the following equations are valid :(a) 1 = 2 1 - 2 ,(b)

23、k = - k-1 + 2k - k+1 , 1<k<n , (c) n = -n-1 + n . 例 : 找一个矩阵A=2 0 0 0 ; -1 2 0 0; 0 0 2 1; 0 0 0 2 的 Jordan Form. A =2 0 0 0; -1 2 0 0; 0 0 2 1; 0 0 0 2; eigvl= 2 ; n= 4; A1=A-eigvl*eye(n); %A1=A-I deta=zeros(1,n); mu=zeros(1,n); for k=1:1:n A2=A1k ; %A2=(A-I)k z=null(A2,'r'); %basis fo

24、r the null space of A1k deta(k)=size(z,2); %dim ker(A-I)k end ; mu(1)= 2*deta(1)-deta(2) ; for k=2:n-1 mu(k)=-deta(k-1)+2*deta(k)-deta(k+1); end ; mu(n)=-deta(n-1)+deta(n); mu %显示不同大小之Jordan block的个数 mu = 0 2 0 0 从mu(1)=mu(3)=mu(4)=0, mu(2)=2 , 我们得知只有两个2*2的 Jordan Blocks, 因此其Jordan Form 为 1 1 0 0 J

25、= 0 1 0 0 0 0 1 1 0 0 0 1 4. 矩阵的QR decomposition 所谓矩阵的QR decomposition, 即是将矩阵A分解成A=Q*R, 其中Q 是垂直 矩阵 (orthonormal matrix ), R 是上三角矩阵 (upper triangular matrix ) . 我们利用 三种方法 , 说明如下 : (1) Gram-Schmidt process consider A= V1 V2 . . . Vn , each Vi is a column vector 采用下列公式 : U1=V1 , Uk =Vk - 1kU1 - -k-1,k

26、Uk-1 , 2kn ji = (Uj, Vi) / (Uj ,Uj ) if Uj 0 (0 向量 ) = 0 if Uj = 0 得 Q0 = U1 U2 Un . Q0 has orthogonal columns, some of which may be equal to 0 . 满足A=Q0*R0 , R0 = 1 12 13 1n 0 1 23 2n 0 0 1 3n . . . 0 0 . . . 1 n-1,n 0 0 . . . 0 1 R0 is unit-upper-triangular and nonsingular. 我们若从Q0中删除0 的行向量, 而且从R0 中

27、删除对应 的列向量, 然后 Q0 中每个非0的行向量再除以其2-norm, 而R0 中对应 的列向量乘上 此值, 则可得到Q , 单位垂直矩阵(orthonormal matrix)与R, 上三角矩阵 (upper-triangular) . 例 : 分解A = 1 2 0 -1 1 -1 3 2 1 -1 3 2 -1 1 -3 1 , 使得A=Q*R A = 1 2 0 -1; 1 -1 3 2; 1 -1 3 2; -1 1 -3 1 ; m, n = size(A) ; rk=rank(A); R=zeros(rk, n) ; Q=zeros(m,rk) ; i=0 ; %给一些初值

28、for k = 1:1:n , v = A(1:m, k) ; if v = zeros(m,1) i = i+1 ; R(i,k)=norm(v) ; %填入上三角矩阵主对角上的元素 Q(1:m, i) = v / norm(v) ; %填入垂直矩阵第i行向量 for j = k+1:1:n, R(i, j) = Q(1:m, i)'* A(1:m, j); %填入上三角矩阵第(i,j)的元素 A(1:m, j) = A(1:m, j)- R(i, j)*Q(1:m, i); %计算Q第j行向量 end %内回圈for结束 end %if- block end %外回圈for结束 Q

29、, R %列印Q, R Q*R Q = 0.5000 0.8660 -0.0000 0.5000 -0.2887 0.4082 0.5000 -0.2887 0.4082 -0.5000 0.2887 0.8165R = 2.0000 -0.5000 4.5000 1.0000 0 2.5981 -2.5981 -1.7321 0 0 0 2.4495ans = 1 2 0 -1 1 -1 3 2 1 -1 3 2 -1 1 -3 1 从ans=Q*R 的结果得知, 我们分解成功. Matlab 的m-file与 function(函数) 由于我们解题的指令趋于复杂, 因此亦可类似于写程式的方式, 将所需的指令 存成一个档案, 即Matlab所谓的m-file; 然后, 再执行这个档案. 若发生错误, 则 修改档案之后, 再来执行. 例如, 我们把刚才利用Gram-Schmidt process分解 矩阵A的方法, 存成一个m-file其file name是 qrdecomp.m