在前面的博客中我提到了如何實現正定矩陣的Cholesky分解，並提供了源代碼，通過該代碼可以將一個正定矩陣分解為一個上三角矩陣和其轉置的乘積，在此基礎上，對上三角矩陣進行求逆是十分簡單的運算，在得到其逆矩陣之后，也就能求出原正定矩陣的逆矩陣了。
數學原理如下：
對於u的逆矩陣，可以使用下列函數進行計算：
[m,n]=size(u);
if m~=n
error('Error using in invutri.Matrix must besquare.');
end
for i=1:n
for j=1:i-1
if u(i,j)~=0
error('Error using in invutri.Matrix must beup-triangle.');
end
end
end
for i=1:n
if u(i,i)==0
error('Error using in invutri.Matrix issingular.')
end
end
mu=[u,eye(n)];
for i=n:-1:1
mu(i,:)=mu(i,:)/mu(i,i);
for j=i-1:-1:1
mu(j,:)=mu(j,:)-mu(i,:)*mu(j,i);
end
end
ut=mu(:,n+1:2*n);
end
這樣，結合前面的工作，我們就得到了兩個m函數文件：cholesky.m和invutri.m文件，分別完成正定矩陣的Cholesky分解和上三角矩陣的求逆。接下來，我們實現正定矩陣的求逆過程，只需要三行代碼即可完成工作：
u=cholesky(A);
tu=invutri(u);
B=tu*tu';
則B就是正定矩陣A的逆矩陣。
測試代碼如下：
>> clear
>> A=[1 2 3 4;2 5 9 10;3 9 22 20;4 10 20 37]
A =
1 2 3 4
2 5 9 10
3 9 22 20
4 10 20 37
>> u=cholesky(A)
u =
1 2 3 4
0 1 3 2
0 0 2 1
0 0 0 4
>> tu=invutri(u)
tu =
1.0000 -2.0000 1.5000 -0.3750
0 1.0000 -1.5000 -0.1250
0 0 0.5000 -0.1250
0 0 0 0.2500
>> B=tu*tu'
B =
7.3906 -4.2031 0.7969 -0.0938
-4.2031 3.2656 -0.7344 -0.0313
0.7969 -0.7344 0.2656 -0.0313
-0.0938 -0.0313 -0.0313 0.0625
>> B*A
ans =
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
>> A*B
ans =
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
表明B確實是A的逆矩陣
這幾段代碼都是使用的MATLAB中的基本語句，而沒有使用inv、chol等MATLAB中現成的函數，從數學原理上分析操作過程。但是也要小心的是，在使用時，cholesky函數並沒有檢驗矩陣是否是正定的，所以計算結果有問題的時候，最好看看矩陣是否是正定的。
您可能也喜歡：
- Cholesky分解
- 對稱陣的LDL分解
- 原來這個就是Given陣~
- 關於Cholesky分解的日志
- Cholesky分解
- 矩陣的QR分解
- Gauss-Seidel迭代法
- 左除，右除
- 矩陣的Doolittle分解
- 金字塔矩陣練習
- 矩陣的擬上三角化（Hessenberg矩陣）
- MATLAB LU函數