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);

if k < n

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

    L(k, k) = 1;

U(k, k) = v(k);

将这段代码运行1500次，每次使用随机生成的高斯矩阵并和Matlab的lu分解结果对比，统计数据如下：

mean

of my lu     : 1.7676e-14

variance

of my lu : 3.07075e-25

of matlab lu     : 3.70416e-16

of matlab lu : 2.14989e-32

和前一种方法不同，目前的这个方法显式的计算出L和U，而在前面的方法是计算出L，然后就可以拿到U。

使用部分选主元的方法进行LU分解，对应如下代码：

function [P, L,

U] = zgaxpyplu(A)

%compute

LU decomposition based on Gaxpy operation with pivoting

%aimed

at improve the stability of LU decomposition

error(‘current

P =

eye(n);

L =

U =

zeros(n);

    v

= zeros(n, 1);

    v(k:end) =

A(k:end, 1);

    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);

%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);

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]);

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);

运行这段代码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中对应的方法是最好的选择。