第4章 矩阵分析
本章包括
◆ 矩阵分析
◆ 线性方程组
◆ 矩阵分解
前面章节已经介绍过，MATLAB的基本运算单元是数组，因此本章内容将从矩阵分析的数值计算开始介绍。在MATLAB中，数值运算主要通过函数或者命令来实现，而由于在MATLAB中所有的数据都是以矩阵的形式出现的，其对应的数值运算包括两种类型：一种是针对整个矩阵的数值运算，也就是矩阵运算，例如求解矩阵行列式的函数det；另外一种是针对矩阵中的元素进行运算的函数，可以称为矩阵元素的运算，例如求解矩阵中每个元素的余弦函数cos等。
4.1 矩阵计算
矩阵分析是线性代数的重要内容，也是几乎所有MATLAB函数分析的基础。在MATLAB 7.0中，可以支持多种线性代数中定义的操作，正是其强大的矩阵运算能力才使得MATLAB成为优秀的数值计算软件。在本节中，将主要介绍关于矩阵分析的内容。
4.1.1 进行范数分析——使用norm函数
根据线性代数的知识，对于线性空间中的某个向量x={x1,x2,…,xn}，其对应的p级范数的定义为
，其中的参数p=1,2,…,n。同时，为了保证整个定义的完整性，定义范数数值
，
。
矩阵范数的定义是基于向量的范数的，具体的表达式为：
在实际应用中，比较常用的矩阵范数是1、2和∞阶范数，其对应的定义如下：
在上面的定义式
中，Smax{ATA}表示的是矩阵A的最大奇异值的平方，关于奇异值的定义将在后面章节中介绍。
说明
之所以在本节中介绍范数分析的内容，是因为矩阵的范数将直接影响通过MATLAB求解得到的数值解的精度。
在MATLAB中，求解向量和矩阵范数的命令如下：
◆ n = norm(A) 计算向量或者矩阵的2阶范数。
◆ n = norm(A,p) 计算向量或者矩阵的p阶范数。
说明
在命令n = norm(A,p)中，p可以选择任何大于1的实数，如果需要求解的是无穷阶范数，则可以将p设置为inf或者-inf。
例4.1 根据定义和norm来分别求解向量的范数。
step 1 进行范数运算。选择命令窗口菜单栏中的“FiIe”→“New”→“M-FiIe”命令，打开M文件编辑器，在其中输入以下程序代码：
%输入向量
x=[1:6];
y=x.^2;
%使用定义求解各阶范数
N2=sqrt(sum(y));
Ninf=max(abs(x));
Nvinf=min(abs(x));
%使用norm命令求解范数
n2=norm(x);
ninf=norm(x,inf);
nvinf=norm(x,-inf);
%输出求解的结果
disp('The method of definition:')
fprintf('The 2-norm is %6.4f\n',N2)
fprintf('The inf-norm is %6.4f\n',Ninf)
fprintf('The minusinf-norm is %6.4f\n',Nvinf)
fprintf('\n//-------------------------------//\n\n')
disp('The method of norm command:')
fprintf('The 2-norm is %6.4f\n',n2)
fprintf('The inf-norm is %6.4f\n',ninf)
fprintf('The minusinf-norm is %6.4f\n',nvinf)
在输入上面的代码后，将其保存为“normex.m”文件。
step 2 查看运算结果。在MATLAB的命令窗口中输入“normex”后，按“Enter”键，得到计算结果如下：
The method of definition:
The 2-norm is 9.5394
The inf-norm is 6.0000
The minusinf-norm is 1.0000
//-------------------------------//
The method of norm command:
The 2-norm is 9.5394
The inf-norm is 6.0000
The minusinf-norm is 1.0000
提示
从以上结果可以看出，根据范数定义得到的结果和norm命令得到的结果完全相同，通过上面的代码可以更好地理解范数的定义。
例4.2 根据定义和norm来分别求解HiIbert矩阵的范数。
step 1 进行范数运算。选择命令窗口菜单栏中的“FiIe”→“New”→“M-FiIe”命令，打开M文件编辑器，在其中输入以下程序代码：
%输入矩阵
A=hilb(5);
%使用定义求解各阶范数
N1=max(sum(abs(A)));
N2=norm(A);
Ninf=max(sum(abs(A')));
Nfro=sqrt(sum(diag(A'＊A)));
%使用norm命令求解范数
n1=norm(A,1);
n2=norm(A,2);
ninf=norm(A,inf);
nfro=norm(A,'fro');
%输出求解的结果
disp('The method of definition:')
fprintf('The 1-norm is %6.4f\n',N1)
fprintf('The 2-norm is %6.4f\n',N2)
fprintf('The inf-norm is %6.4f\n',Ninf)
fprintf('The Ff-norm is %6.4f\n',Nfro)
fprintf('\n//-------------------------------//\n\n')
disp('The method of norm command:')
fprintf('The 1-norm is %6.4f\n',n1)
fprintf('The 2-norm is %6.4f\n',n2)
fprintf('The inf-norm is %6.4f\n',ninf)
fprintf('The F-norm is %6.4f\n',nfro)
在输入上面的代码后，将其保存为“normex2.m”文件。
step 2 查看运算结果。在MATLAB的命令窗口中输入“normex2”，然后按“Enter”键，得到计算结果如下：
The method of definition:
The 1-norm is 2.2833
The 2-norm is 1.5671
The inf-norm is 2.2833
The Ff-norm is 1.5809
//-------------------------------//
The method of norm command:
The 1-norm is 2.2833
The 2-norm is 1.5671
The inf-norm is 2.2833
The F-norm is 1.5809
step 3 查看矩阵A的数值元素。在上面的步骤中，分别使用定义和norm命令来求解HiIbert矩阵A的范数，查看A的具体元素，如下所示：
%定义数值显示格式
>> format rat
>> A
A =
1 1/2 1/3 1/4 1/5
1/2 1/3 1/4 1/5 1/6
1/3 1/4 1/5 1/6 1/7
1/4 1/5 1/6 1/7 1/8
1/5 1/6 1/7 1/8 1/9
说明
在MATLAB中，Hilbert矩阵是著名的病态矩阵，主要用来分析矩阵的性能。用户可以使用hilb命令来创建该矩阵，其元素满足等式A(,i j)=1/(i+ j-1)。
4.1.2 进行范数分析——使用normest函数
当需要分析的矩阵比较大时，求解矩阵范数的时间就会比较长，因此当允许某个近似的范数满足某种条件时，可以使用normest函数来求解范数。在MATLAB的设计中，normest函数主要是用来处理稀疏矩阵的，但是该命令也可以接受正常矩阵的输入，一般用来处理维数比较大的矩阵。
normest函数的主要调用格式如下：
◆ nrm = normest(S) 估计矩阵S的2阶范数数值，默认的允许误差数值维为1e-6。
◆ nrm = normest(S,toI) 使用参数toI作为允许的相对误差。
例4.3 分别使用norm和normest命令来求解矩阵的范数。
step 1 在MATLAB的命令窗口中输入以下命令：
>> W = gallery('wilkinson',500);
t1=clock;
W_norm=norm(W);
t2=clock;
t_norm=etime(t2,t1);
t3=clock;
W_normest=normest(W);
t4=clock;
t_normest=etime(t4,t3);
提示
在上面的程序代码中，首先创建了Wilkinson的高维矩阵，然后分别使用norm和normest命令求解矩阵的范数，并统计了每个命令所使用的时间。
step 2 查看计算得到的矩阵范数。在命令窗口中输入对应的变量名称，得到的结果如下：
W_norm =
250.2462
t_norm =
0.7410
W_normest =
250.2368
t_normest =
0.3210
说明
从以上结果可以看出，两种方法得到的结果几乎相等，但在消耗的时间上，normest命令明显要少于norm命令。
step 3 修改矩阵的维度，重新求解范数。在MATLAB的命令窗口中输入以下命令：
%重新设置矩阵的维度
W = gallery('wilkinson',1000);
t1=clock;
W_norm=norm(W)
t2=clock;
t_norm=etime(t2,t1)
t3=clock;
W_normest=normest(W)
t4=clock;
t_normest=etime(t4,t3)
%显示结果
W_norm =
500.2462
t_norm =
8.4620
W_normest =
500.2116
t_normest =
4.4270
从以上结果可以看出，维度越大则两个命令求解所消耗时间的差别越大，建议在求解大型矩阵的范数时，使用normest命令。
提示
关于每个命令的计算时间，和运行命令的硬件环境以及是否首次运行软件有关，因此读者运行的时间结果可能和这里显示的不同，但是两者之间的大小关系不会变化。
4.1.3 条件数分析
在线性代数中，描述线性方程Ax = b的解对b中的误差或不确定性的敏感度的度量就是矩阵A的条件数，其对应的数学定义是：
k=‖A-1‖·‖A‖
根据基础的数学知识，矩阵的条件数总是大于等于1。其中，正交矩阵的条件数为1，奇异矩阵的条件数为∞，而病态矩阵的条件数则比较大。
依据条件数，方程解的相对误差可以由以下不等式来估算：
在MATLAB中，求取矩阵X的条件数的命令如下：
c = cond(X) 求矩阵X的条件数
例4.4 以MATLAB产生的Magic和HiIbert矩阵为例，使用矩阵的条件数来分析对应的线性方程解的精度。
step 1 进行数值求解。在MATLAB的命令窗口中输入以下命令：
>> M=magic(5);
>> b=ones(5,1);
%利用左除M求解近似解
>> x=M\b;
%准确的求解
>> xinv=inv(M)＊b;
%计算实际相对误差
>> ndb=norm(M＊x-b);
>> nb=norm(b);
>> ndx=norm(x-xinv);
>> nx=norm(x);
>> er=ndx/nx;
>> k=cond(M);
%计算最大可能的近似相对误差
>> erk1=k＊eps;
%计算最大可能的相对误差
>> erk2=k＊ndb/nb;
说明
在上面的程序代码中，首先产生Magic矩阵，然后使用近似解和准确解进行比较，得出计算误差。
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
>> k
k =
5.4618
>> er
er =
2.9403e-016
>> erk1
erk1 =
1.2128e-015
>> erk2
erk2 =
6.6426e-016
说明
从以上结果可以看出，矩阵M的条件数5.418，这种情况下引起的计算误差是很小的，其误差是完全可以接受的。
step 2 修改求解矩阵，重新计算求解的精度。在命令窗口中输入以下代码：
M=hilb(12);
b=ones(12,1);
x=M\b;
xinv=invhilb(12)＊b;
ndb=norm(M＊x-b);
nb=norm(b);
ndx=norm(x-xinv);
nx=norm(x);
er=ndx/nx;
k=cond(M);
erk1=k＊eps;
erk2=k＊ndb/nb;
step 4 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
>> k
k =
1.7945e+016
>> er
er =
0.0119
>> erk1
erk1 =
3.9846
>> erk2
erk2 =
5.3479e+007
说明
从以上结果可以看出，该矩阵的条件数为1.7945e+016，该矩阵在数学理论中就是高度病态的，这样会造成比较大的计算误差。
4.1.4 数值矩阵的行列式
在MATLAB中，求解矩阵行列式的命令比较简单，其调用格式如下：
d = det(X) 求解矩阵X的行列式，如果输入的参数X不是矩阵，而是一个常数，则该命令返回原来的常数。
例4.5 求解矩阵的行列式。
step 1 求解矩阵的行列式。在MATLAB的命令窗口中输入以下命令：
for i=1:3
A=randint(i+2);
a(i)=det(A);
disp('The Matrix is:');
disp(A);
disp('The determinant is');
disp(num2str(a(i)));
end
step 2 查看求解的结果。在输入上面的程序代码后，按“Enter”键，得到的结果如下：
The Matrix is:
1 0 0
0 0 1
0 0 0
The determinant is
0
The Matrix is:
0 1 0 0
1 1 0 1
0 0 1 0
1 0 1 0
The determinant is
1
The Matrix is:
0 0 1 1 1
1 1 0 0 1
0 1 1 1 0
1 0 1 1 1
0 1 1 0 1
The determinant is
-2
说明
在上面的程序代码中，首先使用randint命令产生随机矩阵，然后使用det命令来计算这些矩阵的行列式。
4.1.5 符号矩阵的行列式
需要说明的是，det命令除了可以计算数值矩阵的行列式之外，还可以计算符号矩阵的行列式，下面举例说明。
例4.6 求解符号矩阵的行列式。
step 1 在MATLAB的命令窗口中输入以下命令：
>> syms t;
>>A=[sin(t),cos(t);-cos(t),sin(t)];
>> B=det(A);
>> C=simple(B);
说明
在上面的程序代码中，使用det命令求解矩阵A的行列式表达式B，然后使用simple命令来化简表达式B。
step 2 查看运行的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
A =
[ sin(t), cos(t)]
[ -cos(t), sin(t)]
B =
sin(t)^2+cos(t)^2
C =
1
提示
从以上结果可以看出，使用det命令也可以求解符号矩阵的行列式。关于符号运算的其他命令，可以查看符号运算的章节。
4.1.6 矩阵的化零矩阵
对于非满秩的矩阵A，存在某矩阵Z，满足A·Z=0，同时矩阵Z是一个正交矩阵，也就是说Z ′·Z=I，则矩阵Z被称为矩阵A的化零矩阵。在MATLAB中，求解化零矩阵的命令为nuII，其具体的调用格式如下：
◆ Z = nuII(A) 返回矩阵A的化零矩阵，如果化零矩阵不存在，则返回空矩阵。
◆ Z = nuII(A,'r') 返回有理数形式的化零矩阵。
例4.7 求解非满秩矩阵A的化零矩阵。
step 1 在MATLAB的命令窗口中输入以下命令：
>> A = [1 2 3
4 5 6
7 8 9];
>> Z = null(A);
R=A＊Z;
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
>> Z
Z =
-0.4082
0.8165
-0.4082
>> R
R =
1.0e-015 ＊
0
-0.4441
0
step 3 求解有理数形式的化零矩阵。在MATLAB的命令窗口中输入以下命令：
>> ZR=null(A,'r');
>> RZ=A＊ZR;
step 4 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
>> ZR
ZR =
1
-2
1
>> RZ
RZ =
0
0
0
4.2 线性方程组
线性方程组是线性代数中的主要内容之一，也是理论发展最为完整的部分，在MATLAB中也包含了多种处理线性方程组的命令，限于篇幅，在本节中就不详细展开各种理论和对应的命令，主要处理三种方程：恰定方程组、超定方程组和欠定方程组，下面详细介绍。
4.2.1 非奇异线性方程组
在线性代数中，恰定方程组是指方程组的个数和未知量个数相等的方程组，在恰定方程组中矩阵A是方阵，矩阵B可能是向量也可能是方阵。对于恰定方程组，MATLAB提供了一个十分方便的命令，左除“\”。根据系数矩阵A的奇异属性，该命令可以得到不同的结果：
◆ 如果恰定方程是非奇异的，则左除命令给出了恰定方程组的精确解。
◆ 如果恰定方程组是奇异的，MATLAB会显示提示警告信息，同时给出解NaN。
下面分别使用例子来详细说明如何使用该命令求解方程组。
例4.8 求解非奇异矩阵的线性方程的解。
step 1 在MATLAB的命令窗口中输入以下命令：
>> A = pascal(4);
>> b = [1; 3; 4;6];
>> x=A\b;
>> Bsol=A＊x;
>> D=det(A);
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
A =
1 1 1 1
1 2 3 4
1 3 6 10
1 4 10 20
D =
1
x =
-4
10
-7
2
Bsol =
1
3
4
6
从以上结果可以看出，方程组的系数矩阵A的行列式为1，则系数矩阵A是非奇异的，因此MATLAB可以给出准确的数值解。同时，该数值解和系数矩阵相乘，得到的结果和原来的数值完全相同。
说明
在上面的代码中，首先使用pascal命令创建了数值矩阵，该命令是MATLAB的内置函数，该命令将产生一个对称的正定矩阵，矩阵的数值取自Pascal三角形。
4.2.2 奇异线性方程组
例4.9 求解奇异矩阵的线性方程的解。
step 1 在MATLAB的命令窗口中输入以下命令：
>> A = [ 1 3 7
-1 4 4
1 10 18 ];
>> b=[6;4;15];
>> x=A\b
Warning: Matrix is singular to working precision.
x =
NaN
Inf
-Inf
说明
从以上结果可以看出，MATLAB会显示提示信息，表示该矩阵是奇异矩阵，因此无法得到精确的数值解。
step 2 查看矩阵信息。用户可以使用多种命令来查看矩阵A是否是奇异的，如下：
>> det(A)
ans =
0
>> rank(A)
ans =
2
说明
从上面的结果可以看出，矩阵A的行列式为0，同时矩阵的秩为2，表示该矩阵A是严格奇异的，MATLAB将无法给出精确的数值解。
对于恰定方程组，如果数值解有解，则可以使用矩阵A的伪逆矩阵pinv(A)来得到方程的一个解，其对应的数值解为pinv(A)＊B，下面举例来说明。
例4.10 使用伪逆矩阵的方法求解奇异矩阵的线性方程的解。
step 1 在MATLAB的命令窗口中输入以下命令：
>> A = [ 1 3 7
-1 4 4
1 10 18 ];
>> b=[5;2;12];
>> x=pinv(A)＊b;
>> bsol=A＊x;
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
x =
0.3850
-0.1103
0.7066
bsol =
5.0000
2.0000
12.0000
说明
从以上结果可以看出，通过使用伪逆矩阵的方法，可以求解得到数值解，同时该数值解可以精确地满足结果。
step 3 修改需要求解的数值。在MATLAB的命令窗口中输入以下命令：
>> A = [ 1 3 7
-1 4 4
1 10 18 ];
>> b=[6;4;15];
>> x=pinv(A)＊b;
>> bsol=A＊x;
step 4 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
x =
0.0759
0.3126
0.6647
bsol =
5.6667
3.8333
15.1667
说明
从以上结果可以看出，通过该方法求解得到的结果并不完全满足原来的数值条件，并不具有上面步骤中的精度。
4.2.3 欠定线性方程组
在线性代数的理论中，欠定线性方程组是指方程组的未知量个数多于方程个数的问题，这类问题的解不是唯一解，MATLAB 7.0将首先寻求一个基本解，然后再寻求非零解。从解法的角度来看，MATLAB 7.0采用的是QR分解的方法来求解欠定线性方程组。下面举例来详细说明。
例4.11 求解欠定线性方程组。
step 1 在MATLAB的命令窗口中输入以下命令：
>> A = [ 1 2 3
4 5 6
7 8 9
10 11 12 ];
>>b = [1;3;5;7];
>> [Q,R] = qr(A);
>>y=Q'＊b;
>>xqr = R\y;
>> x=A\b;
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
Warning: Rank deficient, rank = 2, tol = 1.4594e-014.
x =
0.5000
0
0.1667
Warning: Rank deficient, rank = 2, tol = 1.4594e-014.
xqr =
0.5000
0
0.1667
Q =
-0.0776 -0.8331 0.5456 -0.0478
-0.3105 -0.4512 -0.6919 0.4704
-0.5433 -0.0694 -0.2531 -0.7975
-0.7762 0.3124 0.3994 0.3748
R =
-12.8841 -14.5916 -16.2992
0 -1.0413 -2.0826
0 0 -0.0000
0 0 0
说明
从以上结果可以看出，直接通过左除求解得到的数值解和使用QR分析得到的数值解是完全相同的，因此读者可以了解到MATLAB求解欠定方程组的原理。
4.2.4 超定线性方程组
超定线性方程组是指方程组的个数比未知数个数多的情况，对于这种情况，MATLAB提供了伪逆矩阵的方法来求解，关于该方法在前面已经介绍过，本节将利用一个具体的实例来说明如何求解超定线性方程组。
例4.12 求解超定线性方程组。
step 1 在MATLAB的命令窗口中输入以下命令：
>> A = magic(8);
>> A = A(:,1:6);
>> b = 260＊ones(8,1);
>> x = pinv(A)＊b;
>> bsol=A＊x;
>> xs=A\b;
>> bs=A＊xs;
>> nx=norm(x);
>> nxs=norm(xs);
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
A =
64 2 3 61 60 6
9 55 54 12 13 51
17 47 46 20 21 43
40 26 27 37 36 30
32 34 35 29 28 38
41 23 22 44 45 19
49 15 14 52 53 11
8 58 59 5 4 62
x =
1.1538
1.4615
1.3846
1.3846
1.4615
1.1538
Warning: Rank deficient, rank = 3, tol = 1.8829e-013.
xs =
4.0000
5.0000
0
0
0
-1.0000
bsol =
260.0000
260.0000
260.0000
260.0000
260.0000
260.0000
260.0000
260.0000
bs =
260.0000
260.0000
260.0000
260.0000
260.0000
260.0000
260.0000
260.0000
nx =
3.2817
nxs =
6.4807
说明
从以上结果可以看出，尽管左除和伪逆矩阵的方法求解得到的数值解都具有足够的精度，但是使用伪逆矩阵得到的数值解的范数要小。
4.3 矩阵分解
矩阵分解主要是指将一个矩阵分解为几个比较简单的矩阵连乘的形式，无论是在理论上还是在工程应用上，矩阵分解都是十分重要的。在MATLAB中，线性方程组的求解主要基于三种基本的矩阵分解，ChoIesky分解、LU分解和QR分解，对于这些分解MATLAB都提供了对应的函数。除了上面介绍的几种分解之外，在本节中还将介绍奇异值分解和舒尔求解两种比较常见的分解。
4.3.1 ChoIesky分解
ChoIesky分解是把一个对称正定矩阵A分解为一个上三角矩阵R和其转置矩阵的乘积，其对应的表达式为：A=R′R。从理论的角度来看，并不是所有的对称矩阵都可以进行ChoIesky分解，需要进行ChoIesky分解的矩阵必须是正定的。
在MATLAB中，进行ChoIesky分解的是choI命令：
◆ R = choI(X) 其中X是对称的正定矩阵，R是上三角矩阵，使得X=R′R。如果矩阵X是非正定矩阵，该命令会返回错误信息。
◆ [R,p]=choI(X) 该命令返回两个参数，并不返回错误信息。当X是正定矩阵时，返回的矩阵R是上三角矩阵，而且满足等式X=R′R，同时返回参数p=0；当X不是正定矩阵时，返回的参数p是正整数，R 是三角矩阵，且矩阵阶数是 p-1，并且满足等式X(1:p-1,1:p-1)=R＇R。
提示
对对称正定矩阵进行分解在矩阵理论中是十分重要的，用户可以首先对该对称正定矩阵进行Cholesky分解，然后经过处理得到线性方程的解，这些内容将在后面的步骤中通过实例来介绍。
例4.13 对对称正定矩阵进行ChoIesky分解。
step 1 在MATLAB的命令窗口中输入以下命令：
>> n = 5;
>>X = pascal(n)
>>R = chol(X)
>>C=transpose(R)＊R
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
X =
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
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
C =
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
说明
从以上结果中可以看出，R是上三角矩阵，同时满足等式C=RT R=X，表明上面的Cholesky分解过程成功。
step 3 修改矩阵信息。在MATLAB的命令窗口中输入以下命令：
>> X(n,n) = X(n,n)-1
>>[R1,p]=chol(X)
>>C1=transpose(R1)＊R1
>>C2=X(1:p-1,1:p-1)
step 4 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
X =
1 1 1 1 1
1 2 3 4 5
1 3 6 10 15
1 4 10 20 35
1 5 15 35 69
R1 =
1 1 1 1
0 1 2 3
0 0 1 3
0 0 0 1
p =
5
C1 =
1 1 1 1
1 2 3 4
1 3 6 10
1 4 10 20
C2 =
1 1 1 1
1 2 3 4
1 3 6 10
1 4 10 20
提示
从以上结果中可以看出，当将原来正定矩阵的最后一个元素减1后，矩阵将不是正定矩阵，并且满足条件X(1:p-1,1:p- =1) RT R。
4.3.2 使用ChoIesky分解求解方程组
例4.14 使用ChoIesky分解来求解线性方程组。
step 1 在MATLAB的命令窗口中输入以下命令：
>> A=pascal(4);
>>b=[1;4;6;13];
>>x=A\b
>>R=chol(A);
>>Rt=transpose(R);
>>xr=R\(Rt\b)R
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
x =
-9
23
-19
6
xr =
-9
23
-19
6
从以上结果可以看出，使用ChoIesky分解求解得到的线性方程组的数值解和使用左除得到的结果完全相同。其对应的数学原理如下：
对于线性方程组Ax=b，其中 A 是对称的正定矩阵，其A=RTR，则根据上面的定义，线性方程组可以转换为RTRx=b，该方程组的数值为x=R\(RT \b)。
提示
尽管两个方法得到的结果相同，但是其复杂度却不相同。假设A是n×n的方阵，则命令chol的计算复杂度是o( 3n)，而左除命令“\”的计算复杂度为o( 2n )。
4.3.3 不完全ChoIesky分解
对于稀疏矩阵，MATLAB提供了choIinc命令来做不完全的ChoIesky分解，该命令的另外一个重要功能是求解实数半正定矩阵的ChoIesky分解，其具体的调用格式如下：
◆ R = choIinc(X,droptoI) 其中参数X和R的含义和choI命令中的含义相同，其中droptoI表示的是不完全ChoIesky分解的丢失容限，当该参数为0时，则属于完全ChoIesky分解。
◆ R = choIinc(X,options) 其中参数options用来设置该命令的相关参数；具体来讲，options是一个结构体，包含了droptoI、michoI和rdiag三个参数。
◆ R = choIinc(X,'0') 完全ChoIesky分解。
◆ [R,p] = choIinc(X,'0') 和命令choI(X)相同。
◆ R = choIinc(X,'inf') 采用ChoIesky-Infinity方法进行分解，ChoIesky-Infinity方法是基于ChoIesky分解的，但是可以用来处理实半正定分解。
例4.15 使用choIinc命令对矩阵进行ChoIesky分解。
step 1 在MATLAB的命令窗口中输入以下命令：
>> H20 = sparse(hilb(20));
>> [R,p] = chol(H20);
>> Rinf = cholinc(H20,'inf');
>> Rfull=full(Rinf(14:end,14:end))
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
Rfull =
Inf 0 0 0 0 0 0
0 Inf 0 0 0 0 0
0 0 Inf 0 0 0 0
0 0 0 Inf 0 0 0
0 0 0 0 Inf 0 0
0 0 0 0 0 Inf 0
0 0 0 0 0 0 Inf
step 3 检验是否满足分解条件。在MATLAB的命令窗口中输入以下命令：
>> H=full(H20(14:end,14:end));
>>H20R=Rfull'＊Rfull;
step 4 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
H =
0.0370 0.0357 0.0345 0.0333 0.0323 0.0313 0.0303
0.0357 0.0345 0.0333 0.0323 0.0313 0.0303 0.0294
0.0345 0.0333 0.0323 0.0313 0.0303 0.0294 0.0286
0.0333 0.0323 0.0313 0.0303 0.0294 0.0286 0.0278
0.0323 0.0313 0.0303 0.0294 0.0286 0.0278 0.0270
0.0313 0.0303 0.0294 0.0286 0.0278 0.0270 0.0263
0.0303 0.0294 0.0286 0.0278 0.0270 0.0263 0.0256
H20R =
Inf NaN NaN NaN NaN NaN NaN
NaN Inf NaN NaN NaN NaN NaN
NaN NaN Inf NaN NaN NaN NaN
NaN NaN NaN Inf NaN NaN NaN
NaN NaN NaN NaN Inf NaN NaN
NaN NaN NaN NaN NaN Inf NaN
NaN NaN NaN NaN NaN NaN Inf
提示
从以上结果可以看出，尽管cholinc命令可以求解得到分解结果，但是该分解结果并不能保证开始的等式关系。
4.3.4 LU分解
LU分解又被称为是高斯消去法，它可以将任意一个方阵A分解为一个“心理”下三角矩阵L和一个上三角矩阵U的乘积，也就是A=LU。其中，“心理”下三角矩阵的定义为下三角矩阵和置换矩阵的乘积。
在MATLAB中，求解LU分解的命令为Iu，其主要调用格式如下：
◆ [L,U] = Iu(X) 其中X是任意方阵，L是“心理”下三角矩阵，U是上三角矩阵，这三个变量满足的条件式为X=LU。
◆ [L,U,P] = Iu(X) 其中X是任意方阵，L是“心理”下三角矩阵，U是上三角矩阵，P是置换矩阵，满足的条件式为PX=LU。
◆ Y = Iu(X) 其中X是任意方阵，把上三角矩阵和下三角矩阵合并在矩阵Y中给出，满足等式为Y=L+U-I，该命令将损失置换矩阵P的信息。
例4.16 使用Iu命令对矩阵进行LU分解。
step 1 在MATLAB的命令窗口中输入以下命令：
>> A = [ -1 8 -5
9 -1 2
2 -5 7 ];
>> [L1,U1]=lu(A);
>> Al=L1＊U1;
>> x=inv(A);
>> xl=inv(U1)＊inv(L1);
>> d=det(A);
>> dl=det(L1)＊det(U1);
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
L1 =
-0.1111 1.0000 0
1.0000 0 0
0.2222 -0.6056 1.0000
U1 =
9.0000 -1.0000 2.0000
0 7.8889 -4.7778
0 0 3.6620
Al =
-1 8 -5
9 -1 2
2 -5 7
x =
-0.0115 0.1192 -0.0423
0.2269 -0.0115 0.1654
0.1654 -0.0423 0.2731
xl =
-0.0115 0.1192 -0.0423
0.2269 -0.0115 0.1654
0.1654 -0.0423 0.2731
d =
-260
dl =
-260
从以上结果可以看出，方阵的LU分解满足以下等式条件：
A=LU、U-1L-1=A-1和det(A)=det(L)det(U)
step 3 使用三个输出变量的命令形式。在MATLAB的命令窗口中输入以下命令：
>> [L,U,P] = lu(A);
>> Lp=P＊L;
>> Ap=L＊U;
>> Pa=P＊A;
step 4 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
L =
1.0000 0 0
-0.1111 1.0000 0
0.2222 -0.6056 1.0000
U =
9.0000 -1.0000 2.0000
0 7.8889 -4.7778
0 0 3.6620
P =
0 1 0
1 0 0
0 0 1
Ap =
9 -1 2
-1 8 -5
2 -5 7
Pa =
9 -1 2
-1 8 -5
2 -5 7
Lp =
-0.1111 1.0000 0
1.0000 0 0
0.2222 -0.6056 1.0000
从以上结果可以看出，使用三个输出变量的命令满足以下等式关系：
PA =LU和PL ＇=L
在上面的等式中，＇L 表示的是使用两个输出变量求解的LU分解矩阵结果。
step 5 使用LU分解来求解线性方程组。在MATLAB的命令窗口中输入以下命令：
>> b=[2;3;5];
>>xb=A\b;
>>y1 = L\b;
>>xb1=U\yl;
>> y2 = L1\b;
>>xb2=U1\yl;
step 6 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
xb =
0.1231
1.2462
1.5692
xb1 =
0.1231
1.2462
1.5692
xb2 =
0.1231
1.2462
1.5692
说明
从以上结果可以看出，无论是两个输出变量还是三个输出变量的命令得到的结果，都和使用左除命令得到的结果相同。
4.3.5 不完全LU分解
对于稀疏矩阵，MATLAB提供了函数Iuinc进行不完全的LU分解，其调用格式如下：
◆ [L,U] = Iuinc(X,droptoI) 命令中各参数X和R的含义和Iu命令中的含义相同，其中droptoI表示的是不完全LU分解的丢失容限，当该参数为0时，则属于完全LU分解。
◆ [L,U] = Iuinc(X,options) 参数options设置了关于LU分解的各种参数。
◆ [L,U] = Iuinc(X,'0') 0级不完全LU分解。
◆ [L,U,P] = Iuinc(X,'0') 0级不完全LU分解。
例4.17 使用Iuinc命令对稀疏矩阵进行LU分解。
step 1 加载稀疏矩阵，并绘制稀疏矩阵图形。在MATLAB的命令窗口中输入以下命令：
%加载稀疏矩阵
>> load west0479;
>> S = west0479;
>> LU = lu(S);
%绘制稀疏矩阵的图形
>> subplot(1,2,1);
>> spy(S);
>>title('S')
%使用LU求解得到的结果
>> subplot(1,2,2);
>> spy(LU);
>>title('LU')
step 2 查看求解的结果。在输入上面的程序代码后，按“Enter”键，得到的图形如图4.1所示。
图4.1 系数矩阵和LU分解结果图形
提示
由于加载的系数矩阵维度比较大，因此如果直接使用数据查看很难看出或者分辨出矩阵的性质。在上面的实例中，使用了Spy命令来查看矩阵的属性。
step 3 使用Iuinc命令对稀疏矩阵进行LU分解。在命令窗口中输入以下命令：
%进行0级LU分解
>> [L,U,P] = luinc(S,'0');
>> D = (L＊U).＊spones(P＊S)-P＊S;
%打开新的图形窗口
%绘制使用lunic命令得到的结果
>>figure
>>subplot(221);
>>spy(L);
>>title('L:luinc(S,0)')
>> subplot(222);
>>spy(U);
>>title('U:luinc(S,0)')
>>subplot(223);
>>spy(P＊S);
>>title('P＊S')
>>subplot(224);
>>spy(L＊U);
>>title('L＊U')
step 4 查看求解的结果。在输入上面的程序代码后，按“Enter”键，得到的图形如图4.2所示。
图4.2 使用Iuinc命令得到的结果
step 5 使用不同的误差容忍度对稀疏矩阵进行LU分解。在命令窗口中输入以下命令：
>>[IL1,IU1,IP1] = luinc(S,1e-8);
>> [IL2,IU2,IP2] = luinc(S,1e-4);
Warning: Incomplete upper triangular factor has 7 zero diagonals.
It cannot be used as a preconditioner for an iterative method
>> [IL3,IU3,IP3] = luinc(S,1e-2);
Warning: Incomplete upper triangular factor has 22 zero diagonals.
It cannot be used as a preconditioner for an iterative method
>> [IL4,IU4,IP4] = luinc(S,1);
Warning: Incomplete upper triangular factor has 87 zero diagonals.
It cannot be used as a preconditioner for an iterative method
>> figure;
>> subplot(221);
>> spy(IL1＊IU1);
>> title('luinc(S,1e-8)')
>> subplot(222);
spy(IL2＊IU2);
title('luinc(S,1e-4)')
>> subplot(223);
spy(IL3＊IU3);
title('luinc(S,1e-2)')
>> subplot(224);
spy(IL4＊IU4);
title('luinc(S,1)')
step 6 查看求解的结果。在输入上面的程序代码后，按“Enter”键，得到的图形如图4.3所示。
图4.3 使用不同的误差容忍度
step 7 绘制“droptoI-nnz”图形。在本命令中，droptoI表示的是LU分解的丢失容限，nnz则表示系数矩阵中非零元素的个数，用户可以绘制两个变量之间的关系。在MATLAB的命令窗口中输入以下代码：
>> nz(1)=nnz(IL1);
tol=1.e-8;
nz(1)=nnz(IL1);
tol(1)=1.e-8;
nz(2)=nnz(IL2);
tol(2)=1.e-4;
nz(3)=nnz(IL3);
tol(3)=1.e-2;
nz(4)=nnz(IL4);
tol(4)=1;
>> semilogx(tol,nz,'g','LineWidth',1.5)
>> set(gca,'Ytick',[0:650:7800])
set(gca,'Ylim',[07800])
>> title('Drop tolerance vs nnz(luinc(S,droptol))')
xlabel('Drop tolerance')
ylabel('nnz')
>> grid
step 8 查看求解的结果。在输入上面的程序代码之后，按“Enter”键，得到的结果如图4.4所示。
图4.4 丢失容限和nnz的关系
step 9 绘制“droptoI-norm”图形。在本命令中，droptoI表示的是LU分解的丢失容限，norm则表示LU分解的相对误差，用户可以绘制两个变量之间的关系。在MATLAB的命令窗口中输入以下代码：
%计算相对误差
>> normlu(1)=norm(IL1＊IU2-IP1＊S,1)/norm(S,1);
>> normlu(1)=norm(IL1＊IU1-IP1＊S,1)/norm(S,1);
>> normlu(2)=norm(IL2＊IU2-IP2＊S,1)/norm(S,1);
>> normlu(3)=norm(IL3＊IU3-IP3＊S,1)/norm(S,1);
>> normlu(4)=norm(IL4＊IU4-IP4＊S,1)/norm(S,1);
>> loglog(tol,normlu,'g','LineWidth',1.5)
>> title('Drop tolerance vs norm norm(L＊U-P＊S,1)/norm(S,1))')
>>xlabel('Drop tolerance')
>>ylabel('norm')
>> grid
step 10 查看求解的结果。在输入上面的代码后，按“Enter”键，得到的图形如图4.5所示。
图4.5 丢失容限和相对误差的图形
4.3.6 QR分解
矩阵的正交分解又被称为QR分解，也就是将一个m×n的矩阵A分解为一个正交矩阵Q和一个上三角矩阵R的乘积，也就是说A=QR。
在MATLAB中，进行QR分解的命令为qr，其调用格式如下：
◆ [Q,R]=qr(A) 矩阵R和矩阵A的大小相同，Q是正交矩阵，满足等式A=QR，该调用方式适用于满矩阵和稀疏矩阵。
◆ [Q,R] = qr(A,0) 比较经济类型的QR分解。假设矩阵A是一个m×n的矩阵，其中m>n，则命令将只计算前n列的元素，返回的矩阵R是n×n矩阵；如果m≤n，该命令和上面的命令[Q,R] = qr(A)相等。该调用方式适用于满矩阵和稀疏矩阵。
◆ [Q,R,E] = qr(A) 该命令中，Q是正交矩阵，R是上三角矩阵，E是置换矩阵，满足条件关系式A·E =Q·R，该调用方式适用于满矩阵。
例4.18 使用qr命令对矩阵进行QR分解。
step 1 在MATLAB的命令窗口中输入以下命令：
>> A = magic(5);
[Q,R] = qr(A)
C=Q＊R
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
A =
17 24 1 8 15
23 5 7 14 16
4 6 13 20 22
10 12 19 21 3
11 18 25 2 9
Q =
-0.5234 0.5058 0.6735 -0.1215 -0.0441
-0.7081 -0.6966 -0.0177 0.0815 -0.0800
-0.1231 0.1367 -0.3558 -0.6307 -0.6646
-0.3079 0.1911 -0.4122 -0.4247 0.7200
-0.3387 0.4514 -0.4996 0.6328 -0.1774
R =
-32.4808 -26.6311 -21.3973 -23.7063 -25.8615
0 19.8943 12.3234 1.9439 4.0856
0 0 -24.3985 -11.6316 -3.7415
0 0 0 -20.0982 -9.9739
0 0 0 0 -16.0005
C =
17.0000 24.0000 1.0000 8.0000 15.0000
23.0000 5.0000 7.0000 14.0000 16.0000
4.0000 6.0000 13.0000 20.0000 22.0000
10.0000 12.0000 19.0000 21.0000 3.0000
11.0000 18.0000 25.0000 2.0000 9.0000
从以上结果可以看出，矩阵R是上三角矩阵，同时满足等式A=QR，在以下步骤中，将需要证明Q矩阵是正交矩阵。
step 3 证明矩阵Q的正交性。在MATLAB的命令窗口中输入以下命令：
>> detQ=det(Q);
>> for i=1:4
A=Q(:,i);
for j=(i+1):5
B=Q(:,j);
C=A'＊B;
disp(num2str(C))
end
end
step 4 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
detQ =
1.0000
C=
5.5511e-017
5.5511e-017
-2.7756e-017
6.9389e-018
2.7756e-017
0
1.3878e-017
0
-6.9389e-017
-2.2204e-016
从以上结果可以看出，矩阵Q的行列式是1，同时，矩阵中每列数据之间的点乘结果都近似为0，因此Q是正交矩阵。
提示
由于Q是正交矩阵，因此经过QR分解得到的矩阵R和原来的矩阵A是等秩的，可以通过分析R的秩来计算矩阵A的秩。
4.3.7 操作QR分解结果
在MATLAB中，除了提供qr命令之外，还提供了qrdeIete和qrinsert命令来处理矩阵运算的QR分解。其中，qrdeIete的功能是删除QR分解得到矩阵的行或者列；qrinsert的功能则是插入QR分解得到矩阵的行或者列。下面以qrdeIete命令为例，说明如何调用该命令。
◆ [Q1,R1] = qrdeIete(Q,R,j) 返回矩阵A1的QR分解结果，其中A1是矩阵A删除第j列得到的结果，而矩阵A=QR。
◆ [Q1,R1] = qrdeIete(Q,R,j,'coI') 计算结果和[Q1,R1] = qrdeIete(Q,R,j)相同。
◆ [Q1,R1] = qrdeIete(Q,R,j,'row') 返回矩阵A1的QR分解结果，其中A1是矩阵A删除第j行的数据得到的结果，而矩阵A=QR。
例4.19 对矩阵QR分解得到的矩阵进行删除运算。
step 1 在MATLAB的命令窗口中输入以下命令：
A = magic(5);
[Q,R] = qr(A);j = 3;
[Q1,R1] = qrdelete(Q,R,j,'row');
C=Q1＊R1;
A2 = A;
A2(j,:) = [];
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
Q1 =
0.5274 -0.5197 -0.6697 -0.0578
0.7135 0.6911 0.0158 0.1142
0.3102 -0.1982 0.4675 -0.8037
0.3413 -0.4616 0.5768 0.5811
R1 =
32.2335 26.0908 19.9482 21.4063 23.3297
0 -19.7045 -10.9891 0.4318 -1.4873
0 0 22.7444 5.8357 -3.1977
0 0 0 -14.5784 3.7796
C =
17.0000 24.0000 1.0000 8.0000 15.0000
23.0000 5.0000 7.0000 14.0000 16.0000
10.0000 12.0000 19.0000 21.0000 3.0000
11.0000 18.0000 25.0000 2.0000 9.0000
A2 =
17 24 1 8 15
23 5 7 14 16
10 12 19 21 3
11 18 25 2 9
在上面的结果中，满足等式C=Q1 · R1，其中C就是矩阵A删除对应数据行的结果，也就是上面结果中的A2矩阵。
step 3 证明Q1矩阵的正交性。在MATLAB的命令窗口中输入以下命令：
>> detQ1=det(Q1);
>> for i=1:3
A=Q1(:,i);
for j=(i+1):4
B=Q1(:,j);
C=A'＊B;
disp(num2str(C))
end
end
step 4 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
detQ1 =
1.0000
C=
2.7756e-017
1.1102e-016
0
-5.5511e-017
0
0
从以上结果可以看出，矩阵Q1的行列式是1，同时，矩阵中每列数据之间的点乘结果都近似为0，因此Q1是正交矩阵。
说明
在MATLAB中，qrinsert命令和qrdelete命令的调用方法几乎相同，在这里就不详细介绍该命令的用法了，下面将用一个简单的例子说明如何使用qrinsert命令。
例4.20 对矩阵QR分解得到的矩阵进行插入运算。
step 1 在MATLAB的命令窗口中输入以下命令：
>> A = magic(5);
[Q,R] = qr(A);
j = 3;
x = 1:5;
[Q1,R1] = qrinsert(Q,R,j,x,'row');
>> Aqr=Q1＊R1;
>> A2 = [A(1:j-1,:); x; A(j:end,:)];
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
Q1 =
0.5231 0.5039 -0.6750 0.1205 0.0411 0.0225
0.7078 -0.6966 0.0190 -0.0788 0.0833 -0.0150
0.0308 0.0592 0.0656 0.1169 0.1527 -0.9769
0.1231 0.1363 0.3542 0.6222 0.6398 0.2104
0.3077 0.1902 0.4100 0.4161 -0.7264 -0.0150
0.3385 0.4500 0.4961 -0.6366 0.1761 0.0225
R1 =
32.4962 26.6801 21.4795 23.8182 26.0031
0 19.9292 12.4403 2.1340 4.3271
0 0 24.4514 11.8132 3.9931
0 0 0 20.2382 10.3392
0 0 0 0 16.1948
0 0 0 0 0
Aqr =
17.0000 24.0000 1.0000 8.0000 15.0000
23.0000 5.0000 7.0000 14.0000 16.0000
1.0000 2.0000 3.0000 4.0000 5.0000
4.0000 6.0000 13.0000 20.0000 22.0000
10.0000 12.0000 19.0000 21.0000 3.0000
11.0000 18.0000 25.0000 2.0000 9.0000
A2 =
17 24 1 8 15
23 5 7 14 16
1 2 3 4 5
4 6 13 20 22
10 12 19 21 3
11 18 25 2 9
从以上结果中可以看出，满足等式Aqr=Q1 · R1，其中Aqr就是矩阵A删除对应数据行的结果，也就是上面结果中的A2矩阵。
step 3 证明Q1矩阵的正交性。在MATLAB的命令窗口中输入以下命令：
>> detQ1=det(Q1);
for i=1:5
A=Q1(:,i);
for j=(i+1):6
B=Q1(:,j);
C=A'＊B;
disp(num2str(C))
end
end
step 4 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
detQ1 =
1.0000
C=
-5.5511e-017
5.5511e-017
0
-6.9389e-018
-4.3368e-018
-5.5511e-017
0
-1.3878e-017
-1.7347e-018
0
0
-3.9899e-017
-2.6368e-016
-4.8572e-017
-3.9031e-017
4.3.8 奇异值分解
奇异值分解在矩阵分析中有着重要的地位，对于任意矩阵A∈Cm×n，存在酉矩阵（Unitray matrix），U=[u1,u2,…,un]，V =[v1,v2,…,v n ]，使得：
U T AV =diag(σ 1,σ2,…,σp)
其中参数σ1≥σ2≥…≥σp，p=min{m,n}。在上面的式子中，{σi,ui,vi}分别是矩阵 A的第i个奇异值、左奇异值和右奇异值，它们的组合就称为奇异值分解三对组。
在MATLAB中，计算奇异值分解的命令如下：
◆ [U,S,V] = svd(X) 奇异值分解。
◆ [U,S,V] = svd(X,0) 比较经济的奇异值分解。
◆ s = svds(A,k,0) 向量s中包含矩阵A分解得到的k个最小奇异值。
◆ [U,S,V] = svds(A,k,0) 给出矩阵A的k个最大奇异值分解对应的矩阵。
例4.21 对矩阵进行奇异值分解。
step 1 在MATLAB的命令窗口中输入以下命令：
>> X=[1 2
3 4
5 6
7 8];
>> [U,S,V] = svd(X)
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
U =
-0.1525 -0.8226 -0.3945 -0.3800
-0.3499 -0.4214 0.2428 0.8007
-0.5474 -0.0201 0.6979 -0.4614
-0.7448 0.3812 -0.5462 0.0407
S =
14.2691 0
0 0.6268
0 0
0 0
V =
-0.6414 0.7672
-0.7672 -0.6414
step 3 使用最经济的方法进行分解。在MATLAB的命令窗口中输入以下命令：
>> X=[1 2
3 4
5 6
7 8];
>> [U,S,V] = svd(X,0)
step 4 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
U =
-0.1525 -0.8226
-0.3499 -0.4214
-0.5474 -0.0201
-0.7448 0.3812
S =
14.2691 0
0 0.6268
V =
-0.6414 0.7672
-0.7672 -0.6414
例4.22 使用svd和svds命令对稀疏矩阵进行奇异值分解。
step 1 在MATLAB的命令窗口中输入以下命令：
>> load west0479
s = svd(full(west0479));
sl = svds(west0479,4);ss = svds(west0479,6,0);
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
sl =
1.0e+005 ＊
3.1895
3.1725
3.1695
3.1685
ss =
1.0e-004 ＊
0.5616
0.5169
0.4505
0.4020
0.0424
0.0098
说明
在上面的程序结果中，sl表示的是该稀疏矩阵的4个最大奇异值，ss则是表示该稀疏矩阵的6个最小奇异值。
step 3 绘制数据结果图形。在MATLAB的命令窗口中输入以下命令：
>> plot(sl,'ro')
>> hold on
>> plot(s,'gp')
>> set(gca,'Xtick',[0:1:5])
>> set(gca,'Xlim',[0 4.5])
>> xlabel('n')
>> ylabel('The singular value')
step 4 查看求解的结果。输入上面的程序代码后，按“Enter”键，得到的结果如图4.6所示。
图4.6 最大4个奇异值
说明
从图4.6中可以看出，使用svd和svds命令求解的结果在某个误差容限之内，结果十分相近。
4.4 特征值分析
在线性代数的理论中，对于n×n方阵A，其特征值λ和特征向量x满足以下等式：
Ax= λx
其中，λ是一个标量，x是一个向量。把矩阵A的n个特征值放置在矩阵的对角线上就可以组成一个矩阵D，也就是D=diag(λ1,λ2,…,λn)，然后将各特征值对应的特征向量按照对应次序排列，作为矩阵V的数据列。如果该矩阵V是可逆的，则关于特征值的问题可以描述为：
A·V =V ·D⇒A=V ·D·V-1
在MATLAB中，提供了多种关于矩阵特征值处理的函数，可以使用这些函数来分析矩阵的特征值的多种内容，下面分别详细介绍。
4.4.1 特征值和特征向量
在MATLAB中，求解矩阵特征值和特征向量的数值运算方法简单描述为：对矩阵进行一系列的House-hoIder变换，产生一个准上三角矩阵，然后使用QR法迭代进行对角化。对于一般读者来讲，可以不用了解这些计算原理。关于矩阵的特征值和特征向量的命令比较简单，具体的调用格式如下：
◆ d= eig(A) 仅计算矩阵A的特征值，并且以向量的形式输出。
◆ [V,D] = eig(A) 计算矩阵A的特征向量矩阵V和特征值对角阵D，满足等式AV =VD。
◆ [V,D] = eig(A,'nobaIance') 当矩阵A中有截断误差数量级相差不大时，该指令更加精确。
◆ [V,D] = eig(A,B) 计算矩阵A的广义特征向量矩阵V和广义特征值对角阵D，满足等式AV =BVD。
◆ d = eigs(A,k,sigma) 计算稀疏矩阵A的k个由sigma指定的特征向量和特征值，关于参数sigma的取值，请读者查看相应的帮助文件。
提示
当只需了解矩阵的特征值的时候，推荐使用第一条命令，这样可以节约系统的资源，同时可以有效地得出结果。
例4.23 对基础矩阵求解矩阵的特征值和特征向量。
step 1 对矩阵进行特征值分析。在MATLAB的命令窗口中输入以下命令：
>> A=pascal(5);
>> [V D]=eig(A)
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
V =
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
D =
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
step 3 检验分析得到的结果。在MATLAB的命令窗口中输入以下命令：
>> detV=det(V);
>> S=A＊V-V＊D;
step 4 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
detV =
1.0000
S =
1.0e-015 ＊
0.0620 -0.0154 0.0304 0.1110 0
0.1057 -0.0271 0.0694 0.0833 0
0.0649 -0.0717 0.0900 -0.0555 0.1110
0.0684 -0.0283 0.0416 0.1527 0.1110
0.0753 0.0933 -0.0191 0.0555 0.1110
提示
从以上结果中可以看出，V矩阵的行列式为1，是可逆矩阵，同时求解得到的矩阵结果满足等式AV =VD。
例4.24 计算当矩阵的元素和截断误差相当的时候，矩阵的特征值和特征向量。
step 1 对矩阵进行特征值分析。在MATLAB的命令窗口中输入以下命令：
>> B = [ 3 -4 -1.9 5＊eps
-2 3 2 -eps/5
-eps/3 eps/2 -1 0
-1.5 -.5 .1 1 ];
>> [VB,DB] = eig(B);
>> ER1=B＊VB - VB＊DB;
>> [VN,DN] = eig(B,'nobalance');
>> ER2=B＊VN - VN＊DN;
说明
在上面的程序代码中，ER1表示的是默认的误差，ER2则表示的是经过“nobalance”参数处理后的误差。
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
VB =
0.8016 -0.3920 0.0000 0.0420
-0.5668 -0.2772 0.0000 0.4408
-0.0000 -0.0000 0.0000 -0.8396
-0.1903 -0.8772 1.0000 -0.3148
VN =
1.0000 -0.4469 -0.0000 0.0500
-0.7071 -0.3160 -0.0000 0.5250
-0.0000 -0.0000 -0.0000 -1.0000
-0.2374 -1.0000 -1.0000 0.2187
DB =
5.8284 0 0 0
0 0.1716 0 0
0 0 1.0000 0
0 0 0 -1.0000
DN =
5.8284 0 0 0
0 0.1716 0 0
0 0 1.0000 0
0 0 0 -1.0000
ER1 =
0.0000 0.0000 0.0000 -0.0000
-0.0000 -0.0000 -0.0000 0.0000
0.0000 -0.0000 0.0000 -0.0000
-0.0000 -0.0000 -0.0000 -0.9970
ER2 =
1.0e-014 ＊
0.1776 0.0569 0.0290 0.0701
0.0888 -0.0430 -0.0173 -0.0333
0.0006 0.0021 0.0020 -0.0222
0 -0.0167 0 0
说明
从以上结果可以看出，如果矩阵中的元素和截断误差相当时，如果在命令中不选用“nobalance”参数，则求取的结果是有严重计算错误的。当矩阵A来自程序中的中间计算结果而且结果中包含了eps元素时，就会发生上面的情况。
4.4.2 稀疏矩阵的特征值和特征向量
例4.25 使用eigs命令求取稀疏矩阵的特征值和特征向量。
step 1 生成稀疏矩阵，并求取特征值。在MATLAB的命令窗口中输入以下命令：
>>A = delsq(numgrid('C',30));
>>d = eig(full(A));
>>[dum,ind] = sort(abs(d));
>>dlm = eigs(A);
>>dsm = eigs(A,6,'sm');
>>dsmt=sort(dsm);
>>subplot(2,1,1)
>>plot(dlm,'r+')
>>hold on
>>plot(d(ind(end:-1:end-5)),'rs')
>>hold off
>>legend('eigs(A)','eig(full(A))',3)
>>set(gca,'XLim',[0.5 6.5])
>>grid
>>title('Six largest magnitude eigenvalues')
>>subplot(2,1,2)
>>plot(dsmt,'r+')
>>hold on
>>plot(d(ind(1:6)),'rs')
>>hold off
>>legend('eigs(A,6,''sm'')','eig(full(A))',2)
>>grid
>>set(gca,'XLim',[0.5 6.5])
>>title('Six samllest magnitude eigenvalues')
step 2 查看求解的结果。在输入了上面的程序代码后，按“Enter”键，得到的图形如图4.7所示。
图4.7 计算的图形结果
说明
在本例中，矩阵A是一个对称的正定矩阵，维度是632，其对应的特征值均匀分布在(0 8)上，但是特征值4重复出现了18次。
4.4.3 特征值问题的条件数
在前面的章节中，曾经介绍过如果在MATLAB中求解代数方程的条件数，这个命令不能用来求解矩阵特征值对扰动的灵敏度。矩阵特征值条件数定义是对矩阵的每个特征值进行的，其具体的定义如下：
在上面的等式中，νli，νri分别是特征值λi所对应的左特征行向量和右特征列向量。其中，θ(·,·)表示的是两个向量的夹角。
在MATLAB中，计算特征值条件数的命令如下：
◆ c = condeig(A) 向量c中包含了矩阵A中关于各特征值的条件数。
◆ [V,D,s] = condeig(A) 该命令相等于[V,D] = eig(A)和c = condeig(A)的组合。
例4.26 使用命令分别求解方程组的条件数和特征值条件数。
step 1 在MATLAB的命令窗口中输入以下命令：
>>A=magic(10);
>>cequ=cond(A);
>>ceig=condeig(A);
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
cequ =
2.0660e+018
ceig =
1.0000
3.0261
4.9128
4.7390
1.0802
2.8862
1.4891
2.3247
2.3247
3.6460
说明
从以上结果来看，方程的条件数很大，但是矩阵特征值的条件数则比较小，这就表明了方程的条件数和对应矩阵特征值条件数是不等的。
step 3 重新计算新的矩阵，进行分析。在MATLAB的命令窗口中输入以下命令：
>> A=eye(5,5);
>> A(3,2)=1;
>>A(2,5)=1;
>> cequ=cond(A);
>>ceig=condeig(A);
step 4 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.465190e-032.
A =
1 0 0 0 0
0 1 0 0 1
0 1 1 0 0
0 0 0 1 0
0 0 0 0 1
cequ =
4.0489
ceig =
1.0e+031 ＊
0.0000
0.0000
2.0282
0.0000
2.0282
说明
从以上结果中可以看出，在上面的例子中方程组的条件数很小，而对应的特征值条件数则有两个分量，相当大。
例4.27 对亏损矩阵进行条件数分析。
step 1 在MATLAB的命令窗口中输入以下命令：
>> A=gallery(5);
>> [V,D,c_eig]=condeig(A)
>> condA=cond(A)
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
V =
0.0000 -0.0000+0.0000i -0.0000-0.0000i 0.0000+0.0000i 0.0000-0.0000i
0.0206 0.0206+0.0001i 0.0206-0.0001i 0.0207+0.0001i 0.0207-0.0001i
0.1398 -0.1397+0.0001i -0.1397-0.0001i -0.1397+0.0000i -0.1397-0.0000i
-0.9574 0.9574 0.9574 0.9574 0.9574
-0.2519 0.2519-0.0000i 0.2519+0.0000i 0.2519-0.0000i 0.2519+0.0000i
D =
-0.0408 0 0 0 0
0 -0.0119 + 0.0386i 0 0 0
0 0 -0.0119-0.0386i 0 0
0 0 0 0.0323 + 0.0230i 0
0 0 0 0 0.0323-0.0230i
c_eig =
1.0e+010 ＊
2.1293
2.0796
2.0796
2.0020
2.0020
condA =
2.0253e+018
说明
在上面的步骤中，矩阵A是代数重复度为5、几何重复度为1的亏损矩阵，这是一个十分特殊的矩阵，其方程组的条件数和特征值的条件数都很大。eig命令是不能适用于这类矩阵的，所求得的结果是不可信的。
4.4.4 特征值的复数问题
在理论中，即使是实数矩阵，其对应的特征值也有可能是复数，在实际应用中，经常需要将一对共轭复数特征值转换为一个实数块，为此MATLAB提供了以下命令：
◆ [VR,DR] = cdf2rdf(VC,DC) 将复数对角形转换成实数对角形。
◆ [VC,DC] = rsf2csf(VR,DR) 将实数对角形转换成复数对角形。
说明
以上命令的参数中，DC表示的是含有复数的特征值对角阵，VC表示对应的特征向量矩阵，DR表示的是含有实数的特征值对角阵，VR表示对应的特征向量矩阵。
例4.28 对矩阵的复数特征值进行分析。
step 1 在MATLAB的命令窗口中输入以下命令：
>> X=[ 1 2 3
0 4 5
0 -5 4];
>> [VC,DC] = eig(X);
>> [VR,DR] = cdf2rdf(VC,DC);
>> XR=VR＊DR/VR;
>> XC=VC＊DC/VC;
step 2 查看求解的结果。在命令窗口中输入计算的变量名称，得到的结果如下：
VC =
1.0000 -0.0191-0.4002i -0.0191 + 0.4002i
0 0-0.6479i 0 + 0.6479i
0 0.6479 0.6479
DC =
1.0000 0 0
0 4.0000 + 5.0000i 0
0 0 4.0000-5.0000i
VR =
1.0000 -0.0191 -0.4002
0 0 -0.6479
0 0.6479 0
DR =
1.0000 0 0
0 4.0000 5.0000
0 -5.0000 4.0000
XC =
1.0000 2.0000 + 0.0000i 3.0000
0 4.0000 5.0000
0 -5.0000 4.0000
XR =
1.0000 2.0000 3.0000
0 4.0000 5.0000
0 -5.0000 4.0000
4.5 小结
矩阵是MATLAB最重要的数值对象，几乎所有复杂的MATLAB操作都是基于矩阵的。因此，熟悉常见的矩阵分析是十分必要的。本章主要介绍了矩阵的常见运算、方程组的求解、矩阵分解和特征值的分析等。在讲解这些内容的时候，都分别讲解了这些分析方法的理论背景以及适用范围，并结合具体的例子进行讲解。