4.1  因式分解

本节介绍线性代数的一些基本操作，包括行列式、逆和秩，LU分解和QR分解，以及范数等。其中LU分解和QR分解都是使用对角线上方或者下方的元素均为0的三角矩阵来进行计算。使用三角矩阵表示的线性方程组，可以通过向前或者向后置换很容易地得出结果。

4.1.1  行列式、逆和秩

在MATLAB中，用户可以通过以下命令来计算矩阵A的行列式、逆和矩阵的秩。

(1)det(A)：求方阵A的行列式。

(2)rank(A)：求矩阵A的秩，即A中线性无关的行数和列数。

(3)inv(A)：求方阵A的逆矩阵。如果A是奇异矩阵或者近似奇异矩阵，则会给出一个错误信息。

(4)pinv(A)：求矩阵A的伪逆。如果A是m×n的矩阵，则伪逆的尺寸为n×m。对于非奇矩阵A来说，有pinv(A)=inv(A)。

(5)trance(A)：求矩阵A的迹，也就是对角线元素之和。

下面用具体示例对矩阵行列式、逆和秩作简要的说明。

【例4-1】  矩阵的行列式计算示例。

det函数可以用来计算矩阵的行列式。

>> A1=[1 2;3 4]          %  创建矩阵A1

1     2

3     4

>> A2=[1 2;3,6]           %  创建矩阵A2

1     2

3     6

>> A3=[1 2;3 4;5 6]      %  创建矩阵A3

1     2

3     4

5     6

>> det1=det(A1)           %  求方阵A1的行列式

>> det2=det(A2)           % 求方阵A2的行列式

>> det3=det(A3)           %  注意非方阵的行列式没有意义

??? Error using ==> det

Matrix must be square.

>> det_1=det(A1')        %  实数矩阵的行列式和它的转置的行列式相同

>> det_2=det(A2')

>> det_3=det(A3')

??? Error using ==> det

Matrix must be square.

本例使用det函数计算A3的行列式时返回了错误信息，提醒用户A3必须是是方阵才可以调用det函数。

【例4-2】  矩阵的逆计算示例。

本例在上例创建的矩阵基础上进行演示。

>> inv1=inv(A1)

-2.0000    1.0000

1.5000   -0.5000

>> inv2=inv(A2)              %  A2行列式为0，不存在逆矩阵

Warning: Matrix is singular to working precision.

Inf   Inf

Inf   Inf

>> inv3=inv(A3)               %  非方阵不存在逆矩阵

??? Error using ==> inv

Matrix must be square.

>> detinv1=det(inv(A1))      %  A1的逆矩阵的行列式就等于1/det(A1)

>> 1/det(A1)

【例4-3】  使用pinv函数计算矩阵的伪逆示例。

>> pinv1=pinv(A1)   % A1的逆矩阵和它的伪逆是一样的

-2.0000    1.0000

1.5000   -0.5000

>> pinv2=pinv(A2)

0.0200    0.0600

0.0400    0.1200

>> pinv3=pinv(A3)

-1.3333   -0.3333    0.6667

1.0833    0.3333   -0.4167

本例调用pinv函数计算了矩阵A1、A2、A3的伪逆。因为矩阵A2行列式为0，矩阵A3不是方阵，所以不能求矩阵A2和A3的逆，但是可以求这两个矩阵的伪逆。

【例4-4】  使用rank函数求解矩阵的秩示例。

>> rank1=rank(A1)

>> rank2=rank(A2)

>> rank3=rank(A3)

>> rank_1=rank(A1')

>> rank_2=rank(A2')

>> rank_3=rank(A3')

从本例中可以看出矩阵的秩和它的转置的秩相同。

通过上面的这4个例子，可以总结出以下规律。

(1)只有方阵的行列式才有意义。

(2)只有方阵的逆才有意义，但如果方阵的行列式为0，该方阵则不存在逆矩阵。

(3)如果方阵的逆矩阵存在，它的伪逆和逆相同。

(4)如果方阵的逆矩阵存在，它的逆矩阵的行列式det(A-1)等于1/det(A)。

(5)矩阵的秩和它的转置的秩相同。

(6)实数矩阵的行列式和它的转置矩阵的行列式相同。

4.1.2  Cholesky因式分解

分解法是将原方阵分解成一个上三角形矩阵(或是以不同次序排列的上三角阵)和一个下三角形矩阵，这样的分解法又称为三角分解法，它主要用于简化大矩阵行列式值的计算过程、求逆矩阵和求解联立方程组。需要注意的是：这种分解法所得到的上下三角阵并不是唯一的，可以找到多个不同的上下三角阵对，每对三角阵相乘都会得到原矩阵。

对线性系统的求解，MATLAB依据系数矩阵A的不同，而相应地使用不同的方法。如有可能，MATLAB会先分析矩阵的结构。例如A是对称且正定的，则使用Cholesky分解。

如果没有找到可以替代的方法，则采用高斯消元法和部分主元法，主要是对矩阵进行LU因式分解或LU分解。这种方法就是令A＝LU，其中A是方阵，U是一个上三角矩阵，L是一个带有单位对角线的下三角矩阵。

Cholesky因式分解是把一个对称正定矩阵A分解为一个上三角矩阵R和其转置矩阵的乘积，其对应的表达式为：A=R'R。从理论上说，并不是所有的对称矩阵都可以进行Cholesky因式分解，只有正定矩阵才可以。

Pascal矩阵就是一个有趣的例子。下面以Pascal矩阵为例，说明Cholesky因式分解的使用方法。

【例4-5】  Cholesky因式分解示例。

>> A = pascal(6)             % 创建Pascal矩阵

1     1    1     1     1     1

1     2    3     4     5    6

1     3    6    10    15   21

1     4   10    20    35   56

1     6   21    56   126  252

矩阵A的元素是二项式系数，每一个元素都是上方和左方两个元素的和。在MATLAB中，进行Cholesky因式分解的是chol函数。矩阵A的Cholesky因式分解可以通过以下命令得到：

>> R = chol(A)

1     1    1     1     1    1

0     1    2     3     4    5

0     0    1     3     6   10

0     0    0     1     4   10

0     0    0     0     1    5

0     0    0     0     0    1

得到的矩阵R的元素同样也是二项式系数。

Cholesky因式分解允许线性方程组A x = b被R’Rx=b代替。在MATLAB环境中，这个线性方程组可以通过以下命令来求解：

>> x=R\(R'\b)

4.1.3  LU因式分解

LU分解法主要用于简化大矩阵行列式值的计算过程、求逆矩阵和求解联立方程组。需要注意的是：这种分解法得到的上下三角阵对并不是唯一的，可以找到多个不同的上下三角阵对，每对三角阵相乘都会得到原矩阵。

在MATLAB中，求矩阵A的LU分解的调用函数是lu，调用格式如下：

[L,U]=lu(A)

另外，矩阵A的LU分解为线性系统A*x=b提供了以下表达式来快速求解：

【例4-6】  矩阵A的LU分解示例。

>> A=[5 2 0;2 6 2;5 6 7]

5     2    0

2     6    2

5     6    7

>> [L,U]=lu(A)               %  分解所得L是带有单位对角线的下三角矩阵，U是上三角矩阵

1.0000         0         0

0.4000   1.0000         0

1.0000    0.7692    1.0000

5.0000    2.0000         0

0    5.2000    2.0000

0         0    5.4615

>> L*U                          %  验证结果

5     2    0

2     6    2

5     6    7

【例4-7】  矩阵A的LU分解实例。

>> A=[1 2 3;4 5 6;7 8 9];

>> [L,U]=lu(A);

>> B=[9 8 7;6 5 4; 3 2 1];

>> x=U\(L\B)

Warning: Matrix is close to singular or badlyscaled.

Results may be inaccurate. RCOND =1.e-017.

-27   -26  -17

42    41   24

-16   -16   -8

>> A*x           % 验证结果

9     8    7

6     5    4

3     2    1

4.1.4  QR因式分解

如果A是正交矩阵，那么它满足A’A=1。二维坐标旋转变换矩阵就是一个简单的正交矩阵：

矩阵的正交分解又称做QR分解，是将矩阵分解成一个单位正交矩阵和上三角形矩阵。假设A是m×n的矩阵，那么A就可以分解成：

其中Q是一个正交矩阵，R是一个维数和A相同的上三角矩阵，因此Ax=B可以表示为QRx=B或者等同于Rx=QB。这个方程组的系数矩阵是上三角的，因此容易求解。

在MATLAB中，用户可以调用函数qr来求QR因式分解，这个命令可用于分解m×n的矩阵，假设A是m×n的矩阵。qr函数常用调用格式有以下几种。

(1)[Q,R]=qr(A)：求得m×m阶矩阵Q和m×n阶上三角矩阵R。Q和R满足A=QR。

(2)[Q,R,P]= qr(A)：求得矩阵Q，上三角矩阵R和置换矩阵P。R的对角线元素按大小降序排列，且满足AP=QR。

(3)[Q,R]= qr(A,0)：求矩阵A的QR因式分解。如果在m×n的矩阵A中行数小于列数，则给出Q的前n列。

(4)[Q1,R1]=gradelete(Q,R,j)：求去掉矩阵A中第j列之后形成的矩阵的QR因式分解，矩阵Q和R是A的QR因子。

(5)[Q1,R1]=qrinset(Q,R,b,j)：求在矩阵A的第j列前插入列向量b后形成的矩阵的QR因式分解，矩阵Q和R是A的QR因子。如果j=n+1，那么插入的一列放在最后。

【例4-8】  QR分解示例。已知魔方矩阵

对其进行QR分解。

用户只需要调用qr函数就可以实现对A进行QR分解。具体过程如下：

>> A=magic(3)

8     1    6

3     5    7

4     9    2

>> [Q,R]=qr(A)               %  QR分解

-0.8480    0.5223   0.0901

-0.3180   -0.3655   -0.8748

-0.4240   -0.7705    0.4760

-9.4340   -6.2540   -8.1620

0   -8.2394   -0.9655

0         0   -4.6314

>> Q*R                     %  验证结果

8.0000    1.0000    6.0000

3.0000    5.0000    7.0000

4.0000    9.0000    2.0000

【例4-9】  利用QR分解求线性方程组的解。

用户可以通过求A的QR分解并计算R\Q'*B来求解X。具体过程如下：

>> A=[1 2 2;3 2 2;1 1 2]

1     2     2

3     2    2

1     1    2

>> B=[7;9;5]

>> [Q,R]=qr(A)

-0.3015    0.9239  -0.2357

-0.9045   -0.3553  -0.2357

-0.3015    0.1421   0.9428

-3.3166   -2.7136  -3.0151

0    1.2792   1.4213

0         0   0.9428

>>  X=R\Q'*B

>> A\B

4.1.5  范数

向量的范数是一个标量，用来衡量向量的长度。需注意不要把向量范数和向量中元素的个数相混淆。在MATLAB中，可以用命令norm得到不同的范数。

norm形式的表达式还有norm(x,-inf)，但它不是求向量的范数，而是求向量x的绝对值的最小值，即min(abs(x))。请读者注意区分。

【例4-10】  向量范数的求解示例。

>> x=[2 4 5]

2     4    5

>> norm1=norm(x)                 %  欧几里德范数

>> norm2=norm(x,1)                %  1-范数

>> norm3=norm(x,inf)             % ¥-范数

>> norm4=norm(x,4)               %  p-范数

>> norm5=norm(x,-inf)            % 向量绝对值的最小值