LU分解(2)
接着上次LU分解的讲解，这次给出使用不同的计算LU分解的方法，这种方法称为基于GaxPy的计算方法。这里需要了解lapapck中的一些函数。lapack中有一个函数名为gaxpy，所对应的矩阵计算公式是：x = Gx + y; 对应的Matlab代码如下：
function[L, U] =zgaxpylu(A)
%calculate LU decomposition based on Gaxpy operation
%the same way as zlu.m but differnt approach
[m, n] = size(A);
if m ~= n
error('current support square matrix only')
end
L = eye(n);
U = eye(n);
for k = 1:n
v = zeros(n, 1);
if k == 1
v(k:end) = A(k:end, k);
else
U(1:k-1,k) = L(1:k-1,1:k-1)\A(1:k-1,k);
v(k:end) = A(k:end, k) - L(k:end, 1:k-1)*U(1:k-1,k);
end
if k < n
L(k+1:end,k) = v(k+1:end)/v(k);
L(k, k) = 1;
end
U(k, k) = v(k);
end
将这段代码运行1500次，每次使用随机生成的高斯矩阵并和Matlab的lu分解结果对比，统计数据如下：
mean of my lu : 1.7676e-14
variance of my lu : 3.07075e-25
mean of matlab lu : 3.70416e-16
variance of matlab lu : 2.14989e-32
function [P, L, U] = zgaxpyplu(A)
%compute LU decomposition based on Gaxpy operation with pivoting
%aimed at improve the stability of LU decomposition
[m, n] = size(A);
if m ~= n
error('current support square matrix only')
end
P = eye(n);
L = eye(n);
U = zeros(n);
for k = 1:n
v = zeros(n, 1);
if k == 1
v(k:end) = A(k:end, 1);
else
U(1:k-1, k) = L(1:k-1, 1:k-1) \ A(1:k-1, k);
v(k:n) = A(k:n, k) - L(k:n, 1:k-1)*U(1:k-1, k);
end
%find the largest element in v(k:end)
[max_value, max_index] = max(v(k:end));
max_index = max_index + k - 1;
if max_index ~= k
%exchange the max_index-th row and k-th row
A([k max_index], k+1:n) = A([max_index k], k+1:n);
%exchange the max_index-th row and k-th row in L
if k > 1
L([k max_index], 1:k-1) = L([max_index k], 1:k-1);
end
P([k max_index], :) = P([max_index k], :);
v([k max_index]) = v([max_index k]);
end
if (abs(v(k)) > 1e-6) && (k < n)
v(k+1:end) = v(k+1:end) / v(k);
L(k+1:end, k) = v(k+1:end);
end
U(k, k) = v(k);
end
运行这段代码1500次，每次都使用随机生成的矩阵和Matlab的LU分解做对比，结果如下：
mean of my lu : 2.02547e-15
variance of my lu : 9.06349e-29
mean of matlab lu : 3.65765e-16
variance of matlab lu : 1.98304e-32
将每个算法和同类型的算法相比，使用gaxpy的算法都比之前的方法要好。但是现在有一个问题是使用pivot的算法居然比不使用pivot的算法没有足够显著的改进，当然这也证明了使用基于gaxpy运算可以提高LU分解的稳定性。
总结：LU分解是用于求解线性方程组最基本的矩阵分解方法，对于LU分解可以使用高斯消元法，也可以使用基于gaxpy的运算方法。但是从目前的结果来看，使用gaxpy运算可以极大的提高LU分解的稳定性。但是作为成熟的方法，使用lapack中对应的方法是最好的选择。
