第3章 数组与矩阵——3.4 矩阵运算(2)

3.4  矩阵运算（2）

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(n³ )，而左除运算符“\”的计算复杂度只有O(n² )。

例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