cholesky分解

? ? 接著LU分解繼續往下，就會發展出很多相關但是并不完全一樣的矩陣分解，最后對于對稱正定矩陣，我們則可以給出非常有用的cholesky分解。這些分解的來源就在于矩陣本身存在的特殊的

結構。對于矩陣A，如果沒有任何的特殊結構，那么可以給出A=L*U分解，其中L是下三角矩陣且對角線全部為1，U是上三角矩陣但是對角線的值任意，將U正規化成對角線為1的矩陣，產生分解A = L*D*U， D為對角矩陣。如果A為對稱矩陣，那么會產生A=L*D*L分解。如果A為正定對稱矩陣，那么就會產生A=G*G，可以這么理解G=L*sqrt(D)。

A=L*D*U分解對應的Matlab代碼如下：

function[L, D, U] =zldu(A)

%LDU decomposition of square matrix A. The first step for Cholesky

%decomposition

?

[m, n] = size(A);

if m ~= n

? ? error('support square matrix only')

end

?

L = eye(n);

U = eye(n);

d = zeros(n,1);

?

for k=1:n

?? ?

? ? v = zeros(n, 1);

? ? if k == 1

? ? ? ? v(k:end) = A(k:end, k);

? ? else

? ? ? ? m = L(1:k-1, 1:k-1) \ A(1:k-1, k);

? ? ? ? for j = 1:k-1

? ? ? ? ? ? U(j, k) = m(j) / d(j);

? ? ? ? end

?? ? ? ?

? ? ? ? v(k:end) = A(k:end, k) - L(k:end, 1:k-1)*m(:);

? ? end

?? ?

? ? d(k) = v(k);

?? ?

? ? if k < n

? ? ? ? L(k+1:end, k) = v(k+1:end)/v(k);

? ? end

?? ?

end

?

D = diag(d);

分解的穩定性和精度結果如下：

mean of my lu ? ? : 9.0307e-15

variance of my lu : 4.17441e-27

mean of matlab lu ? ? : 3.70519e-16

variance of matlab lu : 2.07393e-32

這里的計算是基于Gaxpy，所以穩定性和精確度相當之好。

?

A=L*D*L分解對應代碼如下，這里要求A必須為對稱矩陣：

function[D, L] =zldl(A)

%A = L*D*L' another version of LU decomposition for matrix A

?

[m, n] = size(A);

?

if m ~= n

? ? error('support square matrix only')

end

?

L = eye(n);

d = zeros(n,1);

?

for k=1:n

? ? v = zeros(n,1);

?? ?

? ? for j=1:k-1

? ? ? ? v(j) = L(k, j)*d(j);

? ? end

?? ?

? ? v(k) = A(k,k) - L(k, 1:k-1)*v(1:k-1);

?? ?

? ? d(k) = v(k);

?? ?

? ? L(k+1:end, k) = (A(k+1:end,k) - A(k+1:end, 1:k-1)*v(1:k-1)) / v(k);

end

?

D = diag(d);

對應分解的精確度和穩定度如下：

mean of my lu : 35.264
variance of my lu : 29011.2
mean of matlab lu : 5.88824e-16
variance of matlab lu : 8.40037e-32

使用如下的代碼做測試：

n = 1500;

my_error = zeros(1, n);

sys_error = zeros(1, n);

?

for i = 1:n

? ? test = gensys(5);

? ? [zd, zl] = zldl(test);

? ? [l, d] = ldl(test);

?

? ? my_error(i) = norm(zl*zd*(zl') - test, 'fro');

? ? sys_error(i) = norm(l*d*(l') - test, 'fro');

end

?

fprintf('mean of my lu ? ? : %g\n', mean(my_error));

fprintf('variance of my lu : %g\n', var(my_error));

?

fprintf('mean of matlab lu ? ? : %g\n', mean(sys_error));

fprintf('variance of matlab lu : %g\n', var(sys_error));

對于運算的精度如此之低的原因并不清楚

?

A=G*G’; cholesky分解對應的代碼如下：

function[G] =zgaxpychol(A)

%cholesky decomposition for symmetric positive definite matrix

%the only requirement is matrix A: symmetric positive definite

?

[m, n] = size(A);

?

if m ~= n

? ? error('support square matrix only')

end

?

G = eye(n);

?

for k=1:n

?? ?

? ? v = A(:,k);

?? ?

? ? if k > 1

? ? ? ? v(:) = v(:) - G(:,1:k-1)*G(k,1:k-1)';

? ? end

?? ?

? ? G(k:end, k) = v(k:end) / sqrt(v(k));

end

對應的測試結果如下

mean of my lu : 1.10711e-15
variance of my lu : 3.04741e-31
mean of matlab lu : 5.5205e-16
variance of matlab lu : 9.64928e-32

自己代碼的精確度和穩定性可以媲美Matlab的代碼，產生這種結果的原因應該是positive sysmetric definite matrix的原因，這段代碼基于gaxpy的結果，下面給出另外一種基于外積的運算結果。

function[G] =zopchol(A)

%cholesky decomposition based on rank-1 matrix update

?

[m, n] = size(A);

if m ~= n

? ? error('support square matrix only')

end

?

G = zeros(n);

?

for k=1:n

?? ?

? ? G(k,k) = sqrt(A(k,k));

? ? G(k+1:end, k) = A(k+1:end, k) / G(k,k);

?? ?

? ? %update matrix A

? ? for j = (k+1):n

? ? ? ? A(k+1:end,j) = A(k+1:end,j) - G(j,k)*G(k+1:end,k);

? ? end

end

?

對應的測試結果如下：

mean of my lu : 9.33114e-16
variance of my lu : 1.71179e-31
mean of matlab lu : 9.92241e-16
variance of matlab lu : 1.60667e-31

對應的測試程序如下，這里使用系統自帶的chol函數完成cholesky分解。

n = 1500;

my_error = zeros(1, n);

sys_error = zeros(1, n);

?

for i = 1:n

? ? test = genpd(5);

? ? [zg] = zopchol(test);

? ? l = chol(test, 'lower');

?

? ? my_error(i) = norm(zg*(zg') - test, 'fro');

? ? sys_error(i) = norm(l*(l') - test, 'fro');

end

?

fprintf('mean of my lu ? ? : %g\n', mean(my_error));

fprintf('variance of my lu : %g\n', var(my_error));

?

fprintf('mean of matlab lu ? ? : %g\n', mean(sys_error));

fprintf('variance of matlab lu : %g\n', var(sys_error));

將兩個結果想比較，可以發現兩個版本的cholesky分解的精確度和穩定度差不多。

Cholesky分解的核心在于矩陣對稱正定的結構，基于LU分解的再次擴展。