zoukankan      html  css  js  c++  java
  • LU分解(1)

    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 ~=

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

  • 相关阅读:
    单例模式
    HashSet、LinkedHashSet、SortedSet、TreeSet
    ArrayList、LinkedList、CopyOnWriteArrayList
    HashMap、Hashtable、LinkedHashMap
    andrew ng machine learning week8 非监督学习
    andrew ng machine learning week7 支持向量机
    andrew ng machine learning week6 机器学习算法理论
    andrew ng machine learning week5 神经网络
    andrew ng machine learning week4 神经网络
    vue组件监听属性变化watch方法报[Vue warn]: Method "watch" has type "object" in the component definition. Did you reference the function correctly?
  • 原文地址:https://www.cnblogs.com/lacozhang/p/3693203.html
Copyright © 2011-2022 走看看