第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的矩陣，其中>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操作都是基于矩陣的。因此，熟悉常見的矩陣分析是十分必要的。本章主要介紹了矩陣的常見運算、方程組的求解、矩陣分解和特征值的分析等。在講解這些內容的時候，都分別講解了這些分析方法的理論背景以及適用范圍，并結合具體的例子進行講解。