matlab中的lu分解,LU分解 - 奋斗吧，骚年 - OSCHINA - 中文开源技术交流社区

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,:)

[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函数的效率最高。这是因为内部已经将其编译成了可执行的机器码。由于matlab是解释执行的编程语言，且其自身的特点是对于矩阵和向量的运算能力较强。因此在Matlab 中应尽量用向量、矩阵运算，不要用 For 循环。

MATLAB (matrix laboratory) is a multi-paradigm numerical computing environment and fourth-generation programming language. A proprietary programming language developed by MathWorks, MATLAB allows matrix manipulations, plotting of functions and data, implementation of algorithms, creation of user interfaces, and interfacing with programs written in other languages, including C, C++, Java, Fortran and Python.