LU分解

function [L,U,p] = lutxloops(A)
%LU Triangular factorization
%   [L,U,p] = lup(A) produces a unit lower triangular matrix L,
%   an upper triangular matrix U and a permutation vector p,
%   so that L*U = A(p,:)

[n,n] = size(A);
p = (1:n)‘;

for k = 1:n-1

   % Find index of largest element below diagonal in k-th column
   m = k;
   for i = k+1:n
      if abs(A(i,k)) > abs(A(m,k))
         m = i;
      end
   end

   % Skip elimination if column is zero
   if (A(m,k) ~= 0)
   
      % Swap pivot row
      if (m ~= k)
         for j = 1:n;
            A([k m],j) = A([m k],j);
         end
         p([k m]) = p([m k]);
      end

      % Compute multipliers
      for i = k+1:n;
         A(i,k) = A(i,k)/A(k,k);
      end

      % Update the remainder of the matrix
      for j = k+1:n;
         for i = k+1:n;
            A(i,j) = A(i,j) - A(i,k)*A(k,j); 
         end
      end
   end
end
    
% Separate result
L = tril(A,-1) + eye(n,n);
U = triu(A);

之前分析过gauss elimination，其核心是LU分解。本段代码保存为lutxloops.m以上代码不再带有线性方程的右端项而已，本质上无区别。

L = tril(X,k) returns the elements on and below the kth diagonal of X. k = 0 is the main diagonal, k > 0 is above the main diagonal, and k < 0 is below the main diagonal.

U = triu(X) returns the upper triangular part of X.

以下代码是lutx.m,可以看到，少用了许多for循环。接下来就是测试这两个功能相同的程序，其各自的运行效率了。

function [L,U,p] = lutx(A)
%LUTX  Triangular factorization, textbook version
%   [L,U,p] = lutx(A) produces a unit lower triangular matrix L,
%   an upper triangular matrix U, and a permutation vector p,
%   so that L*U = A(p,:)

%   Copyright 2014 Cleve Moler
%   Copyright 2014 The MathWorks, Inc.

[n,n] = size(A);
p = (1:n)‘;

for k = 1:n-1

   % Find index of largest element below diagonal in k-th column
   [r,m] = max(abs(A(k:n,k)));
   m = m+k-1;

   % Skip elimination if column is zero
   if (A(m,k) ~= 0)
   
      % Swap pivot row
      if (m ~= k)
         A([k m],:) = A([m k],:);
         p([k m]) = p([m k]);
      end

      % Compute multipliers
      i = k+1:n;
      A(i,k) = A(i,k)/A(k,k);

      % Update the remainder of the matrix
      j = k+1:n;
      A(i,j) = A(i,j) - A(i,k)*A(k,j); 
   end
end

% Separate result
L = tril(A,-1) + eye(n,n);
U = triu(A);

以下是测试结果：

>> n = 525, A = randn(n,n); tic, lutxloops(A); toc

n =

525

Elapsed time is 2.314373 seconds.

>> n = 525, A = randn(n,n); tic, lutx(A); toc

n =

525

Elapsed time is 0.584623 seconds.

可见，运行效率差距是明显的。

>> n = 1920, A = randn(n,n); tic, lu(A); toc

n =

1920

Elapsed time is 0.191691 seconds.

内置的lu函数的效率最高。