526互联

线性代数 · 矩阵 · Matlab | Cholesky 分解代码实现

发布时间 2023-11-01 14:28:58作者: MoonOut

（搬运外网的代码，非原创；原网址 ）

（其实是专业课作业，但感觉国内博客没有合适的代码实现，所以就搬运到自己博客了）

背景 - Cholesky 分解：

  • 若 A 为 n 阶实对称正定矩阵，则存在非奇异下三角矩阵 L，使得 A = LL^T。
    • 是特殊的 LU 分解（下三角 上三角分解）。
    • 若限定 L 的对角元素为正，则这种分解是唯一的。
    • 计算 L 的方法：待定系数，硬算（貌似是这样的）。

Cholesky 分解法的 matlab 函数代码：

% Cholesky 分解法
function [retval] = cholesky_decompose(A)
n = length(A);
L = zeros(n,n);
for i=1:n
   L(i,i) = sqrt(A(i,i)-L(i,:)*L(i,:)');
   for j=(i+1):n
      L(j,i) = (A(j,i)-L(i,:)*L(j,:)')/L(i,i);
   end
end
retval = L';
end

matlab 脚本代码：

X = rand(3,3)
A = X*X'

L1 = cholesky_decompose(A)
L2 = chol(A)

运行结果：

octave:5> source("my_script.m")
X =

   0.010982   0.924195   0.472230
   0.749800   0.704545   0.495176
   0.480342   0.828100   0.611450

A =

   1.0773   0.8932   1.0593
   0.8932   1.3038   1.2464
   1.0593   1.2464   1.2903

L1 =

   1.0379   0.8606   1.0207
        0   0.7505   0.4904
        0        0   0.0902

L2 =

   1.0379   0.8606   1.0207
        0   0.7505   0.4904
        0        0   0.0902

与 matlab 的 chol 函数结果一致。