3.4 矩阵运算（2）
3.4.2 矩阵分解
矩阵分解是把一个矩阵分解成几个“较简单”的矩阵连乘的形式。无论是在理论上还是在工程应用上，矩阵分解都是十分重要的。本节将介绍几种矩阵分解方法，相关函数如表3-3所示。
表3-3 矩阵分解函数
函数 |
功能描述 |
chol |
Cholesky分解 |
cholinc |
稀疏矩阵的不完全Cholesky分解 |
lu |
矩阵LU分解 |
luinc |
稀疏矩阵的不完全LU分解 |
qr |
正交三角分解 |
svd |
奇异值分解 |
gsvd |
一般奇异值分解 |
schur |
舒尔分解 |
在MATLAB中，线性方程组的求解主要基于4种基本的矩阵分解，即对称正定矩阵的Cholesky分解、一般方阵的高斯消去法分解、矩形矩阵的正交分解和舒尔分解。
1．对称正定矩阵的Cholesky分解
Cholesky分解在MATLAB中用函数chol()来实现，其常用的调用方式如下。
● R=chol(X)：其中X为对称正定矩阵，R是上三角矩阵，使得X=R'∙R。如果X是非正定的，则结果将返回出错信息。
● [R,p]=chol(X)：返回两个参数，并且不会返回出错信息。当X是正定矩阵时，返回的上三角矩阵R满足X=R'∙R，且p=0；当X是非正定矩阵时，返回值p是正整数，R是上三角矩阵，其阶数为p-1，且满足X(1:p-1,1:p-1)=R'∙R。
考虑线性方程组Ax=b，其中A可以做Cholesky分解，使得A=R'∙R，这样线性方程组就可以改写成R'∙R∙x=b。由于左除运算符“\”可以快速处理三角矩阵，因此得出：
x=R\(R'\b)
如果A是n×n的方阵，则chol(A)的计算复杂度是O(n3 )，而左除运算符“\”的计算复杂度只有O(n2 )。
例3-45：利用chol函数进行矩阵分解示例。
在命令行窗口中依次输入：
clear all; A = pascal(5) % 产生5阶帕斯卡矩阵 E = eig(A) % 计算矩阵A的特征值 R = chol(A) % 进行分解 B = R' * R % 计算R'·R
输出结果：
A = 1 1 1 1 1 1 2 3 4 5 1 3 6 10 15 1 4 10 20 35 1 5 15 35 70 E = 0.0108 0.1812 1.0000 5.5175 92.2904 R = 1 1 1 1 1 0 1 2 3 4 0 0 1 3 6 0 0 0 1 4 0 0 0 0 1 B = 1 1 1 1 1 1 2 3 4 5 1 3 6 10 15 1 4 10 20 35 1 5 15 35 70
例3-46：计算稀疏矩阵的Cholesky因子，并使用置换输出创建具有较少非零元素的Cholesky因子。
在命令行窗口中依次输入：
load west0479; % 基于westo479矩阵创建一个稀疏正定矩阵 A = west0479; s =A'*A; %用两种不同的方法计算矩阵的Cholesky 因子。首先指定两个输出，然后指定三个输出以支持行和列重新排序 [R,flag] = chol(S); [RP,flagP,P]=chol(S); if ~flag & &~flagP % 号对于每次计算,都检查flag = 0以确认计算成功 disp ( 'Factorizations successful. ') else disp ( 'Factorizations failed.') end % 比较chol(S）和经过重新排序的矩阵 chol(P'*S*P)中非零元素的个数 subplot (1,2,1) spy(R) title ( 'Nonzeros in chol(S)') subplot (1,2,2) spy(RP) title ( 'Nonzeros in chol(P''*S*P)')
对稀疏矩阵分解进行图形化显示，如图3-3所示。
图3-3 稀疏矩阵分解图形化显示
2．一般方阵的高斯消去法分解
高斯消去法分解又称LU分解，它可以将任意一个方阵A分解为一个下三角矩阵L和一个上三角矩阵U的乘积，即A=LU。LU分解在MATLAB中用函数lu()来实现，其调用方式如下。
● [L,U]=lu(X)：X为一个方阵，L为下三角矩阵，U为上三角矩阵，满足关系X=L∙U。
● [L,U,P]=lu(X)：X为一个方阵，L为下三角矩阵，U为上三角矩阵，P为置换矩阵，满足关系P∙X=L∙U。
● Y=lu(X)：X为一个方阵，把上三角矩阵和下三角矩阵合并在矩阵Y中给出，矩阵Y的对角元素为上三角矩阵的对角元素，即Y=L+U-I。置换矩阵P的信息丢失。
考虑线性方程组Ax=b，其中，对矩阵A可以做LU分解，使得A=L∙U，这样线性方程组就可以改写成L∙U∙x=b。由于左除运算符“\”可以快速处理三角矩阵，因此可以快速解出：
x=U\(L\b)
利用LU分解来计算行列式的值和矩阵的逆，其命令形式如下：
● det(A)=det(L)*det(U)。
● inv(A)=inv(U)*inv(L)。
例3-47：进行LU分解示例。
在命令行窗口中依次输入：
clear all A = [2 4 5;8 9 6;1 3 5] [L1,U1] = lu(A) [L2,U2,P] = lu(A) Y1 = lu(A) L1 * U1 % 验证 Y2 = L2 + U2 -eye(size(A)) % 验证
输出结果：
A = 2 4 5 8 9 6 1 3 5 L1 = 0.2500 0.9333 1.0000 1.0000 0 0 0.1250 1.0000 0 U1 = 8.0000 9.0000 6.0000 0 1.8750 4.2500 0 0 -0.4667 L2 = 1.0000 0 0 0.1250 1.0000 0 0.2500 0.9333 1.0000 U2 = 8.0000 9.0000 6.0000 0 1.8750 4.2500 0 0 -0.4667 P = 0 1 0 0 0 1 1 0 0 Y1 = 8.0000 9.0000 6.0000 0.1250 1.8750 4.2500 0.2500 0.9333 -0.4667 ans = 2 4 5 8 9 6 1 3 5 Y2 = 8.0000 9.0000 6.0000 0.1250 1.8750 4.2500 0.2500 0.9333 -0.4667
此外，对于稀疏矩阵，MATLAB提供了函数luinc()来做不完全LU分解，其调用格式如下。
● [L U]=luinc(X,DROPTOL)：其中X、L和U的含义与函数lu()中的变量相同，DROPTOL为不完全LU分解的丢失容限。当DROPTOL设为0时，退化为完全LU分解。
● [L U]=luinc(X,OPTS)：其中OPTS为结构体，它有4个属性，即DROPTOL、MICHOL、RDIAG和THRESH。DROPTOL为不完全LU分解的丢失容限；当MICHOL为1时，采用改进算法的不完全LU分解，否则不采用改进算法；当RDIAG为1时，R的对角元素中的零值替换成DROPTOL的平方根，当其为0时不做此替换；THRESH是绕对角线旋转因子，其取值范围是[0,1]，当THRESH为0时强制绕对角线旋转，THRESH的默认值是1。
● [L,U,P]=luinc(X,'0')：0级不完全LU分解。
● [L,U]=luinc(X,'0')：0级不完全LU分解。
● Y=luinc(X,'0')：0级完全LU分解。
3．矩形矩阵的正交分解
矩形矩阵的正交分解又称QR分解。QR分解把一个m×n的矩阵A分解为一个正交矩阵Q和一个上三角矩阵R的乘积，即A=Q∙R。在MATLAB中QR分解由函数qr()来实现，下面介绍QR分解的调用方式。
● [Q,R]=qr(A)：其中矩阵R为与矩阵A具有相同大小的上三角矩阵，Q为正交矩阵，它们满足A=Q∙R。该调用方式适用于满矩阵和稀疏矩阵。
● [Q,R]=qr(A,0)：为“经济”方式的QR分解。设矩阵A是一个m×n的矩阵，若m>n，则只计算矩阵Q的前n列元素，R为n×n的矩阵；若m≤n，则与[Q,R]=qr(A)效果一致。该调用方式适用于满矩阵和稀疏矩阵。
● [Q,R,E]=qr(A)：R是上三角矩阵，Q为正交矩阵，E为置换矩阵，它们满足A∙E=Q∙R。程序选择一个合适的矩阵E使得abs(diag(R))是降序排列的。该调用方式适用于满矩阵。
● [Q,R,E]=qr(A,0)：为“经济”方式的QR分解，其中E是一个置换矢量，它们满足A(:,E)=Q∙R。该调用方式适用于满矩阵。
● R=qr(A)：返回上三角矩阵R，这里R=chol(A'∙A)。该调用方式适用于稀疏矩阵。
● R=qr(A,0)：以“经济”方式返回上三角矩阵R。
● [C,R]=qr(A,B)：其中矩阵B必须与矩阵A具有相同的行数，矩阵R是上三角矩阵，C=Q' ∙B。
例3-48：通过QR分解分析矩阵的秩示例。
在命令行窗口中依次输入：
clear all A = [2 4 5;8 9 6;1 3 5]; [Q1,R1] = qr(A) B = [2 4 5;8 9 6;1 3 5;5 4 10]; B_rank = rank(B) disp(['矩阵B的秩 = ',num2str(B_rank)]) [Q2,R2] = qr(B) Q1 * R1
输出结果：
Q1 = -0.2408 0.6424 -0.7276 -0.9631 -0.2511 0.0970 -0.1204 0.7241 0.6791 R1 = -8.3066 -9.9920 -7.5843 0 2.4818 5.3257 0 0 0.3395 B_rank = 3 矩阵B的秩 = 3 Q2 = -0.2063 0.5983 -0.2285 0.7398 -0.8251 0.0774 0.5457 -0.1241 -0.1031 0.6299 -0.3956 -0.6604 -0.5157 -0.4892 -0.7025 0.0348 R2 = -9.6954 -10.6236 -11.6551 0 3.0230 1.7138 0 0 -6.8719 0 0 0 ans = 2.0000 4.0000 5.0000 8.0000 9.0000 6.0000 1.0000 3.0000 5.0000
4．舒尔分解
舒尔分解定义式为
A=U∙S∙U'
其中A必须是一个方阵，U是一个酉矩阵，S是一个块对角矩阵，由对角线上的1×1和2×2块组成。特征值可以由矩阵S的对角块给出，而矩阵U给出比特征向量更多的数值特征。此外，对缺陷矩阵也可以进行舒尔分解。MATLAB中用函数schur()来进行舒尔分解，其调用格式如下。
● [U,S]=schur(A)：返回酉矩阵U和块对角矩阵S。
● S=schur(A)：仅返回块对角矩阵S。
● schur(A,'real')：返回的实特征值放在对角线上，而把复特征值放在对角线上的2×2块中。
● schur(A,' complex')：返回的矩阵S是上三角矩阵，并且如果矩阵A有复特征值，则矩阵S是复矩阵。
另外，函数rsf2csf()可以把实数形式的舒尔矩阵转换成复数形式的舒尔矩阵。
例3-49：舒尔分解示例。
在命令行窗口中依次输入：
clear all A = pascal(5); [U,S] = schur(A) U * S * U' - A % 验证
输出结果：
U = 0.1680 -0.5706 -0.7660 0.2429 0.0175 -0.5517 0.5587 -0.3830 0.4808 0.0749 0.7025 0.2529 0.1642 0.6110 0.2055 -0.4071 -0.5179 0.4377 0.4130 0.4515 0.0900 0.1734 -0.2189 -0.4074 0.8649 S = 0.0108 0 0 0 0 0 0.1812 0 0 0 0 0 1.0000 0 0 0 0 0 5.5175 0 0 0 0 0 92.2904 ans = 1.0e-13 * -0.0033 0.0067 0.0111 0.0133 0.0044 0.0067 0.0133 0.0266 0.0266 0.0178 0.0111 0.0178 0.0355 0.0533 0.0355 0.0089 0.0266 0.0533 0.1421 0.0711 0.0044 0.0266 0.0355 0.1421 0.1421
3.4.3 特征值和特征向量
1．特征值和特征向量的定义
MATLAB中的命令计算特征值和特征向量十分方便，可以得到不同的子结果和分解，这在线性代数学习中十分有意义。本节中的命令只能对二维矩阵进行操作。
假设A是一个n×n的矩阵，A的特征值问题就是找到下面方程组的解：
A∙V=λ∙V
其中，λ为标量，V为矢量，若把矩阵A的n个特征值放在矩阵P的对角线上，相应的特征向量按照与特征值对应的顺序排列，作为矩阵V的列，则特征值问题可以改写为：
A∙V=V∙D
如果V是非奇异的，则该问题可以认为是一个特征值分解问题，此时关系式如下：
A=V∙D∙V-1
广义特征值问题是指方程A∙x=λ∙B∙x的非平凡解问题，其中A、B都是n×n的矩阵，λ为标量。满足此方程的λ为广义特征值，对应的向量x为广义特征向量。
如果X是一个列向量为a的特征向量的矩阵，并且它的秩为n，那么特征向量线性无关。如果不是这样，则称矩阵为缺陷阵。如果X'∙X=I，则特征向量正交，这对于对称矩阵是成立的。
2．特征值和特征向量的相关函数
现将MATLAB中矩阵特征值与特征向量的相关函数的具体调用格式及其功能列出。
● eig(A)：求包含矩阵A的特征值的向量。
● [X,D]=eig(A)：产生一个矩阵A的特征值在对角线上的对角矩阵D和矩阵X，它们的列是相应的特征向量，满足AX=XD。为了得到有更好条件特征值的矩阵，要进行相似变换。
● [T,B]=balance(A)：找到一个相似变换矩阵T和矩阵B，使得它们满足B=T-1∙A∙T。B是用命令balance求得的平衡矩阵。
● eig(A,'nobalance')：不经过平衡处理求得矩阵A的特征值和特征向量，也就是不进行平衡相似变换。
● eigs(A)：返回一个由矩阵A的部分特征值组成的向量，和eig命令一样，但是不返回全部的特征值。如果不带有参量，则计算出最大的特征值。当计算所有特征值时，如果矩阵A的秩不小于6，则计算出6个特征值。
● eigs(f,n)：求出矩阵A的部分特征值。在使用一个矩阵列的线性运算符时，字符串f中包含的是M文件的文件名，n指定问题的阶次。用这种方法来求特征值比开始就用运算符来求要快。
● eigs(A,B,k,sigma)：求矩阵A的部分特征值，矩阵B的大小和A相同；如果没有给出B=eye(size(A))，那么k就是要计算的特征值的个数；如果k没有给出，就用小于6的数或者A的秩。
变量sigma是一个实数或复数的移位参数，或者下列文本字符串中的一个，文本字符串指明的是特征值的属性：“lm”为最大的特征值，“sm”为最小的特征值，“lr”为最大的实数部分，“sr”为最小的实数部分，“be”为同时求得最大和最小的实数部分。
● condeig(A)：返回一个由矩阵A的特征值条件数组成的向量。
● [V,D,s]=condeig(A)：返回[V,D]=eig(A)和s=condeig(A)。
3．特征值和特征向量的计算
例3-50：矩阵特征值和特征向量的计算示例。
在命令行窗口中输入：
A = [0.8 0.2;0.2 0.8]; [Q,d] = eig(A) Q * Q'
输出结果：
Q = -0.7071 0.7071 0.7071 0.7071 d = 0.6000 0 0 1.0000 ans = 1.0000 0.0000 0.0000 1.0000
