1/6 LU 分解

LU 分解可以写成A = LU，这里的L代表下三角矩阵，U代表上三角矩阵。对应的matlab代码如下：

function[L, U] =zlu(A)

% ZLU - LU decomposition for matrix A

% work as gauss elimination

[m, n] = size(A);

if m ~= n

error('Error, current time only support square matrix');

end

L = zeros(n);

U = zeros(n);

for k = 1:n-1

gauss_vector = A(:,k);

gauss_vector(k+1:end) = gauss_vector(k+1:end) ./ gauss_vector(k);

gauss_vector(1:k) = zeros(k,1);

L(:,k) = gauss_vector;

L(k,k) = 1;

for l=k+1:n

A(l,:) = A(l,:) - gauss_vector(l)*A(k,:);

end

end

U = A;

这段代码的目的非常简单，就是使用高斯消元法给出L，U。但是计算的稳定性非常不好，这点可以通过这段代码的分解结果和matlab自带lu的分解结果相比较得出。比较的方法非常简单：就是计算l*u与原始矩阵想减之后的Frobinus范数大小，使用如下的代码做出两个结果的比较：

n = 1000;

my_error = zeros(1, 1000);

sys_error = zeros(1, 1000);

for i = 1:n

test = randn(5);

[zl, zu] = zlu(test);

[l, u] = lu(test);

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

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

end

disp(mean(my_error));

disp(var(my_error));

disp(mean(sys_error));

disp(var(sys_error));

在这段代码里面，随机的生成一个5x5的符合高斯分布的矩阵，然后使用自己写的lu分解和matlab自带的lu分解分别给出L和U,再计算norm(L*U - test)，从这里就可以看出我们自己计算出来的结果精度和matlab自带的lu真实的差异了。这个差异就体现为这些值的均值和方差。结果如下：

mean of my lu : 13.313846

variance of my lu : 43622.114147

mean of matlab lu : 0.000000

variance of matlab lu : 0.000000

从这个结果可以看出，我们自己写的lu分解的结果在均值和方差上比matlab自带的差了很多。个人认为原因有两点：第一个方法的原因，matlab给出的结果是pivoted LU，第二个是因为实现的原因，matlab基于成熟的LAPACK，肯定会比自己写的更好了。

这一步使用PA = LU来完成LU分解。代码如下：

function [P, L, U] = zplu(A)

% pivoted LU decompositon P*A = L*U

[m, n] = size(A);

if m ~= n

error('zplu:test', 'current time only support square matrix');

end

P = eye(n);

L = zeros(n, n);

for k = 1:n-1

%find the largest element in k column of A from row k to n

[max_value, max_index] = max(A(k:end, k));

max_index = max_index + k - 1;

if max_index ~= k

A([k max_index], :) = A([max_index k], :);

P([k max_index], :) = P([max_index k], :);

L([k max_index], :) = L([max_index k], :);

end

if A(k,k) ~= 0

gauss_vector = A(:,k);

gauss_vector(k+1:end) = gauss_vector(k+1:end) ./ gauss_vector(k);

gauss_vector(1:k) = zeros(k,1);

L(:,k) = gauss_vector;

L(k, k) = 1;

for l=k+1:n

A(l,:) = A(l,:) - gauss_vector(l)*A(k,:);

end

end

end

U = triu(A);

下面是运行前面检测程序的输出：

mean of my lu : 7.803258

variance of my lu : 1450.332510

mean of matlab lu : 0.000000

variance of matlab lu : 0.000000

两个结果相比较可以看到，Matlab的lu一样的稳定，但是使用pivot来调整矩阵A的次序可以极大的提高LU分解的稳定度，这个可以从下降了非常多的方差可以看出。

pivot LU是从k列的k+1到n个元素种选择最大的一个，调换到第k个位置。从我个人的角度理解，除以最大的元素使得高斯变换矩阵中非对角元素全部小于1。由于计算机种存储浮点数的机制，绝对值越靠近0，其精度越高。所以使用pivot这种方法可以极大的提高LU分解的稳定程度。但是也需要指出，使用pivot并不一定能提高LU分解的精度，对于特定的矩阵，不使用pivot说不定可以获得更好的性能。

为了进一步提高提高LU分解的稳定性，可以使用full pivoted LU。分解公式：P*A*Q = L * U; 对应的Matlab代码如下：

function [P, Q, L, U] = zflu(A)

%full pivoted LU decomposition

%

% full pivoted LU decomposition

[m, n] = size(A);

if m ~= n

error('current only support square matrix')

end

P = eye(n);

Q = eye(n);

for k=1:n-1

%find the larget element in A(k:n,k:n)

[max_value, row_index] = max(A(k:n, k:n));

[max_value, col_index] = max(max_value);

real_row = k-1 + row_index(col_index);

real_col = k-1 + col_index;

%exchange the row and column of matrix A

if real_row ~= k

A([k real_row],:) = A([real_row k], :);

P([k real_row],:) = P([real_row k], :);

end

if real_col ~= k

A(:, [k real_col]) = A(:, [real_col k]);

Q(:, [k real_col]) = Q(:, [real_col k]);

end

if A(k, k) ~= 0

rows = k+1:n;

A(rows, k) = A(rows, k) ./ A(k, k);

A(rows, rows) = A(rows, rows) - A(rows, k)*A(k, rows);

end

end

L = tril(A);

for k=1:n

L(k, k) = 1;

end

U = triu(A);

跑完test之后的结果如下：

mean of my lu : 7.77222e-16

variance of my lu : 4.3478e-29

mean of matlab lu : 3.69764e-16

variance of matlab lu : 2.03659e-32

可以看到使用full pivoted LU 分解可以在很大程度上保证分解的稳定性，即便是使用我们自己写的代码。但是即便如此，仍然推荐使用LAPACK中的代码，因为那里面的代码是经过严格的测试和分析的，在各种一场情况下应该都有很好的表现。

这里所介绍的LU分解可以使用另外一种基于GaxPy形式的运行，将在下面介绍。