4.1 因式分解
本節(jié)介紹線性代數(shù)的一些基本操作，包括行列式、逆和秩，LU分解和QR分解，以及范數(shù)等。其中LU分解和QR分解都是使用對角線上方或者下方的元素均為0的三角矩陣來進行計算。使用三角矩陣表示的線性方程組，可以通過向前或者向后置換很容易地得出結(jié)果。
4.1.1 行列式、逆和秩
在MATLAB中，用戶可以通過以下命令來計算矩陣A的行列式、逆和矩陣的秩。
（1）det(A)：求方陣A的行列式。
（2）rank(A)：求矩陣A的秩，即A中線性無關(guān)的行數(shù)和列數(shù)。
（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函數(shù)可以用來計算矩陣的行列式。
>> A1=[1 2;3 4] % 創(chuàng)建矩陣A1
A1 =
1 2
3 4
>> A2=[1 2;3,6] % 創(chuàng)建矩陣A2
A2 =
1 2
3 6
>> A3=[1 2;3 4;5 6] % 創(chuàng)建矩陣A3
A3 =
1 2
3 4
5 6
>> det1=det(A1) % 求方陣A1的行列式
det1 =
-2
>> det2=det(A2) % 求方陣A2的行列式
det2 =
0
>> det3=det(A3) % 注意非方陣的行列式?jīng)]有意義
??? Error using ==> det
Matrix must be square.
>> det_1=det(A1') % 實數(shù)矩陣的行列式和它的轉(zhuǎn)置的行列式相同
det_1 =
-2
>> det_2=det(A2')
det_2 =
0
>> det_3=det(A3')
??? Error using ==> det
Matrix must be square.
本例使用det函數(shù)計算A3的行列式時返回了錯誤信息，提醒用戶A3必須是是方陣才可以調(diào)用det函數(shù)。
【例4-2】 矩陣的逆計算示例。
本例在上例創(chuàng)建的矩陣基礎上進行演示。
>> 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
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函數(shù)計算矩陣的偽逆示例。
>> 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
本例調(diào)用pinv函數(shù)計算了矩陣A1、A2、A3的偽逆。因為矩陣A2行列式為0，矩陣A3不是方陣，所以不能求矩陣A2和A3的逆，但是可以求這兩個矩陣的偽逆。
【例4-4】 使用rank函數(shù)求解矩陣的秩示例。
>> rank1=rank(A1)
rank1 =
2
>> rank2=rank(A2)
rank2 =
1
>> rank3=rank(A3)
rank3 =
2
>> rank_1=rank(A1')
rank_1 =
2
>> rank_2=rank(A2')
rank_2 =
1
>> rank_3=rank(A3')
rank_3 =
2
從本例中可以看出矩陣的秩和它的轉(zhuǎn)置的秩相同。
通過上面的這4個例子，可以總結(jié)出以下規(guī)律。
（1）只有方陣的行列式才有意義。
（2）只有方陣的逆才有意義，但如果方陣的行列式為0，該方陣則不存在逆矩陣。
（3）如果方陣的逆矩陣存在，它的偽逆和逆相同。
（4）如果方陣的逆矩陣存在，它的逆矩陣的行列式det(A-1)等于1/det(A)。
（5）矩陣的秩和它的轉(zhuǎn)置的秩相同。
（6）實數(shù)矩陣的行列式和它的轉(zhuǎn)置矩陣的行列式相同。
4.1.2 Cholesky因式分解
分解法是將原方陣分解成一個上三角形矩陣（或是以不同次序排列的上三角陣）和一個下三角形矩陣，這樣的分解法又稱為三角分解法，它主要用于簡化大矩陣行列式值的計算過程、求逆矩陣和求解聯(lián)立方程組。需要注意的是：這種分解法所得到的上下三角陣并不是唯一的，可以找到多個不同的上下三角陣對，每對三角陣相乘都會得到原矩陣。
對線性系統(tǒng)的求解，MATLAB依據(jù)系數(shù)矩陣A的不同，而相應地使用不同的方法。如有可能，MATLAB會先分析矩陣的結(jié)構(gòu)。例如A是對稱且正定的，則使用Cholesky分解。
如果沒有找到可以替代的方法，則采用高斯消元法和部分主元法，主要是對矩陣進行LU因式分解或LU分解。這種方法就是令A＝LU，其中A是方陣，U是一個上三角矩陣，L是一個帶有單位對角線的下三角矩陣。
Cholesky因式分解是把一個對稱正定矩陣A分解為一個上三角矩陣R和其轉(zhuǎn)置矩陣的乘積，其對應的表達式為：A=R'R。從理論上說，并不是所有的對稱矩陣都可以進行Cholesky因式分解，只有正定矩陣才可以。
Pascal矩陣就是一個有趣的例子。下面以Pascal矩陣為例，說明Cholesky因式分解的使用方法。
【例4-5】 Cholesky因式分解示例。
>> A = pascal(6) % 創(chuàng)建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的元素是二項式系數(shù)，每一個元素都是上方和左方兩個元素的和。在MATLAB中，進行Cholesky因式分解的是chol函數(shù)。矩陣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的元素同樣也是二項式系數(shù)。
Cholesky因式分解允許線性方程組A x = b被R’Rx=b代替。在MATLAB環(huán)境中，這個線性方程組可以通過以下命令來求解：
>> x=R\(R'\b)
4.1.3 LU因式分解
LU分解法主要用于簡化大矩陣行列式值的計算過程、求逆矩陣和求解聯(lián)立方程組。需要注意的是：這種分解法得到的上下三角陣對并不是唯一的，可以找到多個不同的上下三角陣對，每對三角陣相乘都會得到原矩陣。
在MATLAB中，求矩陣A的LU分解的調(diào)用函數(shù)是lu，調(diào)用格式如下：
[L,U]=lu(A)
另外，矩陣A的LU分解為線性系統(tǒng)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 % 驗證結(jié)果
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 % 驗證結(jié)果
ans =
9 8 7
6 5 4
3 2 1
4.1.4 QR因式分解
如果A是正交矩陣，那么它滿足A’A=1。二維坐標旋轉(zhuǎn)變換矩陣就是一個簡單的正交矩陣：
矩陣的正交分解又稱做QR分解，是將矩陣分解成一個單位正交矩陣和上三角形矩陣。假設A是m×n的矩陣，那么A就可以分解成：
A=QR
其中Q是一個正交矩陣，R是一個維數(shù)和A相同的上三角矩陣，因此Ax=B可以表示為QRx=B或者等同于Rx=QB。這個方程組的系數(shù)矩陣是上三角的，因此容易求解。
在MATLAB中，用戶可以調(diào)用函數(shù)qr來求QR因式分解，這個命令可用于分解m×n的矩陣，假設A是m×n的矩陣。qr函數(shù)常用調(diào)用格式有以下幾種。
（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中行數(shù)小于列數(shù)，則給出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分解。
用戶只需要調(diào)用qr函數(shù)就可以實現(xiàn)對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 % 驗證結(jié)果
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 =
7
9
5
>> [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 范數(shù)
向量的范數(shù)是一個標量，用來衡量向量的長度。需注意不要把向量范數(shù)和向量中元素的個數(shù)相混淆。在MATLAB中，可以用命令norm得到不同的范數(shù)。
norm形式的表達式還有norm(x,-inf)，但它不是求向量的范數(shù)，而是求向量x的絕對值的最小值，即min(abs(x))。請讀者注意區(qū)分。
【例4-10】 向量范數(shù)的求解示例。
>> x=[2 4 5]
x =
2 4 5
>> norm1=norm(x) % 歐幾里德范數(shù)
norm1 =
6.7082
>> norm2=norm(x,1) % 1-范數(shù)
norm2 =
11
>> norm3=norm(x,inf) % ￥-范數(shù)
norm3 =
5
>> norm4=norm(x,4) % p-范數(shù)
norm4 =
5.4727
>> norm5=norm(x,-inf) % 向量絕對值的最小值
norm5 =
2
