接着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