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形式的运行，将在下面介绍。