3.4 矩阵运算
矩阵运算是线性代数中极其重要的部分。本节将介绍MATLAB中与矩阵运算相关的内容，包括矩阵分析、矩阵分解、特征值和特征向量等。
3.4.1 矩阵分析
MATLAB提供的矩阵分析函数如表3-2所示。
表3-2 矩阵分析函数
1．向量和矩阵的范数运算
对于线性空间中的一个向量x={x1,x2,…,xn}，如果存在一个函数r(x)满足以下3个条件：
（1）r(x)>0，且r(x)=0的充要条件为x=0。
（2）r(ax)=|a|r(x)，其中a为任意标量。
（3）对向量x和y，有r(x+y)≤r(x)+r(y)。
则称r(x)为向量x的范数，一般记为||x||。范数的形式多种多样，下面式子中定义的范数操作就满足以上3个条件：
式中定义的
称为p阶范数，其中最有用的是1、2和∞阶范数。
矩阵的范数是基于向量的范数定义的，其定义式如下：
与向量的范数一样，矩阵的范数最常用的也是1、2和∞阶范数，它们的定义式如下：
在上式中，
则为矩阵A的最大奇异值的平方。
在MATLAB中，求向量范数的函数的具体用法如下。
● N=norm(x,p)：对任意大于1的p值，返回向量x的p阶范数。
● N=norm(x)：返回向量的2阶范数，相当于N=norm(x,2)。
● N=norm(x,inf)：返回向量的∞阶范数，相当于N=max(abs(x))。
● N=norm(x,-inf)：返回向量的-∞阶范数，相当于N=min(abs(x))。
● 在MATLAB中，求矩阵范数的函数的具体用法如下。
● N=norm(A)：计算矩阵的2阶范数，也就是最大奇异值。
● N=norm(A,p)：根据参数p的值不同，求不同阶的范数值。当p=1时，计算矩阵A的1阶范数，相当于max(sum(abs(A)))。当p=2时，计算矩阵A的2阶范数，相当于norm(A)。当p=inf时，计算矩阵A的∞阶范数，相当于max(sum(abs(A')))。当p=pro时，计算矩阵A的F范数（Frobenius范数），相当于sqrt(sum(diag(A'*A)))。
例3-36：求向量x的2阶范数示例。
在命令行窗口中输入：
输出结果：
○注意
当矩阵维数比较大时，会导致计算矩阵范数的时间比较长，并且当一个近似的范数值满足要求时，可以考虑使用函数normest()来估计2阶范数值。函数normest()最初开发时是为了提供给稀疏矩阵使用的，同时它也能接收满矩阵的输入，一般在满矩阵维数比较大时使用。
函数normest()的用法如下。
● normest(S)：估计矩阵S的2阶范数值。
● normest(S,tol)：使用tol作为允许的相对误差。
例3-37：求矩阵的范数示例。
在命令行窗口中依次输入：
输出结果：
2．矩阵的秩
矩阵A中线性无关的列向量个数称为列秩，线性无关的行向量个数称为行秩。在MATLAB中用函数rank()来计算矩阵的秩。函数rank()的用法如下。
● rank(A)：用默认允许误差计算矩阵的秩。
● rank(A,tol)：给定允许误差计算矩阵的秩，tol=max(size(A))·eps(norm(A))。
例3-38：求矩阵的秩示例。
在命令行窗口中依次输入：
输出结果：
3．矩阵的行列式
矩阵A={aij}n×n的行列式定义如下：
其中，k1,k2,…,kn是将序列1,2,…,n交换k次所得的序列。在MATLAB中用函数det()来计算矩阵的行列式。
例3-39：计算矩阵的行列式示例。
在命令行窗口中依次输入：
输出结果：
4．矩阵的迹
矩阵的迹定义为矩阵对角元素之和。在MATLAB中用函数trace()来计算矩阵的迹。
例3-40：计算矩阵的迹示例。
在命令行窗口中依次输入：
输出结果：
5．矩阵的化零矩阵
MATLAB中提供了求化零矩阵的函数null()，其用法如下。
● Z=null(A)：返回矩阵A的一个化零矩阵，如果化零矩阵不存在，则返回空矩阵。
● Z=null(A,'r')：返回有理数形式的化零矩阵。
例3-41：求矩阵的化零矩阵示例。
在命令行窗口中依次输入：
输出结果：
6．矩阵的正交空间
矩阵A的正交空间Q具有Q'·Q=I的性质，并且Q的列矢量构成的线性空间与矩阵A的列矢量构成的线性空间相同，且正交空间Q与矩阵A具有相同的秩。MATLAB中提供了函数orth()来求正交空间Q。
例3-42：矩阵的正交空间求解示例。
在命令行窗口中依次输入：
输出结果：
7．矩阵的约化行阶梯形式
矩阵的约化行阶梯形式是高斯-约旦消去法解线性方程组的结果，其形式为：
MATLAB中提供了函数rref()来求矩阵的约化行阶梯形式，其用法如下。
● R=rref(A)：返回矩阵A的约化行阶梯形式R。
● [R,jb]=rref(A)：返回矩阵A的约化行阶梯形式R，并返回1×r的向量jb，r为矩阵A的秩；A(:,jb)是矩阵A的列矢量构成的线性空间；R(1:r,jb)是r×r的单位矩阵。
● [R,jb]=rref(A,tol)：以tol作为允许的相对误差计算矩阵A的秩。
例3-43：求矩阵A的约化行阶梯形式示例。
在命令行窗口中依次输入：
输出结果：
8．矩阵空间之间的夹角
矩阵空间之间的夹角代表两个矩阵线性相关的程度。如果夹角很小，它们之间的线性相关度就很高；反之，它们之间的线性相关度就很低。在MATLAB中用函数subspace()来实现求矩阵空间之间的夹角，其调用格式如下。
● theta=subspace(A,B)：返回矩阵A和矩阵B之间的夹角。
例3-44：求矩阵A和矩阵B之间的夹角示例。
在命令行窗口中依次输入：
输出结果：
3.4.2 矩阵分解
矩阵分解是把一个矩阵分解成几个“较简单”的矩阵连乘的形式。无论是在理论上还是在工程应用上，矩阵分解都是十分重要的。本节将介绍几种矩阵分解方法，相关函数如表3-3所示。
表3-3 矩阵分解函数
在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函数进行矩阵分解示例。
在命令行窗口中依次输入：
输出结果：
例3-46：计算稀疏矩阵的Cholesky因子，并使用置换输出创建具有较少非零元素的Cholesky因子。
在命令行窗口中依次输入：
对稀疏矩阵分解进行图形化显示，如图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分解示例。
在命令行窗口中依次输入：
输出结果：
此外，对于稀疏矩阵，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分解分析矩阵的秩示例。
在命令行窗口中依次输入：
输出结果：
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：舒尔分解示例。
在命令行窗口中依次输入：
输出结果：
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-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：矩阵特征值和特征向量的计算示例。
在命令行窗口中输入：
输出结果：
