MATLAB因式分解

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

A1 =

1 2

3 4

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

A2 =

1 2

3 6

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

A3 =

1 2

3 4

5 6

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

det1 =

-2

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

det2 =

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

??? Error using ==> det

Matrix must be square.

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

det_1 =

-2

>> det_2=det(A2')

det_2 =

>> det_3=det(A3')

??? Error using ==> det

Matrix must be square.

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

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

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

>> inv1=inv(A1)

inv1 =

-2.0000 1.0000

1.5000 -0.5000

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

Warning: Matrix is singular to working precision.

inv2 =

Inf Inf

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

??? Error using ==> inv

Matrix must be square.

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

detinv1 =

-0.5000

>> 1/det(A1)

ans =

-0.5000

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

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

pinv1 =

-2.0000 1.0000

1.5000 -0.5000

>> pinv2=pinv(A2)

pinv2 =

0.0200 0.0600

0.0400 0.1200

>> pinv3=pinv(A3)

pinv3 =

-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)

rank1 =

>> rank2=rank(A2)

rank2 =

>> rank3=rank(A3)

rank3 =

>> rank_1=rank(A1')

rank_1 =

>> rank_2=rank(A2')

rank_2 =

>> rank_3=rank(A3')

rank_3 =

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

通过上面的这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矩阵

A =

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 5 15 35 70 126

1 6 21 56 126 252

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

>> R = chol(A)

R =

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提供了以下表达式来快速求解：

x=U\(L\b)

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

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

A =

5 2 0

2 6 2

5 6 7

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

L =

1.0000 0 0

0.4000 1.0000 0

1.0000 0.7692 1.0000

U =

5.0000 2.0000 0

0 5.2000 2.0000

0 0 5.4615

>> L*U % 验证结果

ans =

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.586033e-017.

x =

-27 -26 -17

42 41 24

-16 -16 -8

>> A*x % 验证结果

ans =

9 8 7

6 5 4

3 2 1

4.1.4 QR因式分解

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

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

A=QR

其中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)

A =

8 1 6

3 5 7

4 9 2

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

Q =

-0.8480 0.5223 0.0901

-0.3180 -0.3655 -0.8748

-0.4240 -0.7705 0.4760

R =

-9.4340 -6.2540 -8.1620

0 -8.2394 -0.9655

0 0 -4.6314

>> Q*R % 验证结果

ans =

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]

A =

1 2 2

3 2 2

1 1 2

>> B=[7;9;5]

B =

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

Q =

-0.3015 0.9239 -0.2357

-0.9045 -0.3553 -0.2357

-0.3015 0.1421 0.9428

R =

-3.3166 -2.7136 -3.0151

0 1.2792 1.4213

0 0 0.9428

>> X=R\Q'*B

X =

1.0000

2.0000

1.0000

>> A\B

ans =

1.0000

2.0000

1.0000

4.1.5 范数

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

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

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

>> x=[2 4 5]

x =

2 4 5

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

norm1 =

6.7082

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

norm2 =

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

norm3 =

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

norm4 =

5.4727

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

norm5 =

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。