Matlab
implementation of LU Decomposition and Linear Solver
By Nasser Abbasi
Example
runs and output and Matlab verification
This is a MATLAB implementation for LU decomposition, forward substitution, backward substitution, and linear system solver.
The functions written are:
(1) nma_LU.m does LU decomposition with partial pivoting with threshold support.
(2) nma_ForwardSub.m solves Ly=b for y
(3) nma_BackSub.m solves Ux=y for x
(4) nma_LinearSolve.m driver to solve Ax=b for x by calling (1)->(2)->(3)
I've added partial pivoting (P matrix) support to the LU decomposition function. In addition, the LU function accepts an additional argument which allows the user more control on row exchange.
Matlab own lu() function will do row exchange once it encounters a pivot larger than the current pivot. This is a good thing to always try to do. But sometimes if the difference between the pivots is small, a user might not want this feature. Hence I added a 'threshold' second parameter to the nma_LU.m function to indicate how large a difference should exist for a row exchange to occur. Ofcourse a row exchange will always occur if the current pivot is zero and a non-zero pivot exist to do the exchange.
To get the same exact behavior as Matlab lu() simply make this parameter zero.
Below are examples calling the nma_LU, ForwardSub.m, BackSub.m and LinearSolve.m.
In each example below, I verified the output with Matlab (for the LU and Linear solver cases).
>> A=[1 1 2;2 -1 1;1 2
0]
A =
1
1 2
2
-1 1
1
2 0
>> [L,U,P]=nma_LU(A,0)
L =
1.00000000000000 0 0
0.50000000000000 1.00000000000000 0
0.50000000000000 0.60000000000000 1.00000000000000
U =
2.00000000000000 -1.00000000000000 1.00000000000000
0 2.50000000000000 -0.50000000000000
0 0 1.80000000000000
P =
0
1 0
0
0 1
1
0 0
>> [L,U,P]=lu(A)
L =
1.00000000000000 0 0
0.50000000000000 1.00000000000000 0
0.50000000000000 0.60000000000000 1.00000000000000
U =
2.00000000000000 -1.00000000000000 1.00000000000000
0 2.50000000000000 -0.50000000000000
0 0 1.80000000000000
P =
0
1 0
0
0 1
1
0 0
>>
>> A=rand(4);
>> [L,U,P]=nma_LU(A,0)
L =
1.00000000000000 0 0 0
0.01212703756687 1.00000000000000 0 0
0.07119243718995 0.20742768803520 1.00000000000000 0
0.43394327408595 0.19377225100868 0.40879105345917 1.00000000000000
U =
0.81316649730376 0.19872174266149 0.01527392702904 0.46599434167542
0 0.60138257315521 0.74660044907755 0.41299833684006
0 0 0.11623493184110 0.32625386921423
0 0 0 0.51620218784594
P =
0
0 1 0
0
0 0 1
1
0 0 0
0
1 0 0
>> [L,U,P]=lu(A)
L =
1.00000000000000 0 0 0
0.01212703756687 1.00000000000000 0 0
0.07119243718995 0.20742768803520 1.00000000000000 0
0.43394327408595 0.19377225100868 0.40879105345917 1.00000000000000
U =
0.81316649730376 0.19872174266149 0.01527392702904 0.46599434167542
0 0.60138257315521 0.74660044907755 0.41299833684006
0 0 0.11623493184110 0.32625386921423
0 0 0 0.51620218784594
P =
0
0 1 0
0
0 0 1
1
0 0 0
0
1 0 0
>>
>> A=[1 1 2;2 -1 1;1 2
0]
A =
1
1 2
2
-1 1
1
2 0
>> b=[1 2 1];
>> nma_LinearSolve(A,b)
ans =
1
0
0
>> A\b(:)
ans =
1
0
0
>>
>> A=rand(6);
>> b=rand(6,1);
>> nma_LinearSolve(A,b)
ans =
0.59090034220622
-0.56523444269280
0.95687095978224
-0.97248777153372
1.00007995741472
0.24035777097022
>> A\b(:)
ans =
0.59090034220622
-0.56523444269280
0.95687095978223
-0.97248777153372
1.00007995741472
0.24035777097022
>>
>> [L,U,P]=nma_LU(A,0);
>> nma_ForwardSub(L,b)
ans =
0.83849604493808
0.36727512318587
0.12405626870025
-0.14539724685973
0.17813906538571
-0.19809655526705
>>
>> nma_BackSub(U,ans)
ans =
0.29867870305809
-0.84855613142087
0.48347828223154
-1.68311779577975
1.49928530116874
1.53825192677360
>>
function
[L,U,P]=nma_LU(A,threshold)
%function [L,U,P]=nma_LU(A,threshold)
%
%does LU decomposition with permutation matrix
for
%pivoting reorder. Supports threasold
parameter.
%Similar to matlab lu function, but with
little more control to the
%user as to when row exhanges should be made.
%
%INPUT:
%
A: an nxn square matrix
%
threshold:
%
numerical positive value. row exchanges will be made
%
only if the abs difference between the largest pivot and
%
the current pivot is larger than this threshold.
%
Hence setting threashold to be 0 will cause a row exchange
% anytime when there is a larger pivot.
This is the DEFAULT
%
behaviour similar to Matlab lu().
%OUTPUT:
% P:
nxn permutation matrix such that PA=LU
% L:
unity lower triangular matrix
% U:
unity upper triangular
% By Nasser Abbasi
% HW 4, MATH 501, CSUF
%
if nargin ~=2
error '2 parameters are expected'
end
if ~isnumeric(A)
error 'input matrix A must be numeric'
end
if
~isnumeric(threshold)
error 'threshold must be numeric'
end
if threshold<0
error 'threshold must be positive'
end
[nRow,nCol]=size(A);
if nRow ~= nCol
error 'Matrix must be square'
end
%******************************
%* Internal functions to flip
%* rows for row pivoting
%******************************
function
flipRows()
[c,I]=max(abs(A(n:end,n)));
I=I+(n-1);
tmp=A(n,:);
A(n,:)=A(I,:);
A(I,:)=tmp;
%make sure we also flip the L matrix rows to
keep in sync
tmp=L(n,:);
L(n,:)=L(I,:);
L(I,:)=tmp;
%now make the elementary matrix for this move
E(n,:)=0;
E(n,I)=1;
E(I,:)=0;
E(I,n)=1;
end
P=diag(ones(nRow,1));
U=zeros(nRow);
L=zeros(nRow);
for n=1:nRow-1
currentPivot=A(n,n);
E=diag(ones(nRow,1));
maxPivot=max(A(n+1:end,n));
if
abs(currentPivot)<eps %zero, do row exchage always
if abs(maxPivot)<eps %
not possible to exchange
error 'unable to
complete LU decomposition, bad A'
else
flipRows();
end
else %not a zero pivot, but still can exchange,
check threshold
if
abs(currentPivot)<abs(maxPivot)
if
abs(currentPivot-maxPivot)>=threshold
flipRows();
end
end
end
P=P*E; %update the perumtation matrix
for i=n+1:nRow
L(i,n)=A(i,n)/A(n,n);
A(i,n)=0;
for j=n+1:nRow
A(i,j)=A(i,j)-L(i,n)*A(n,j);
end
end
end
L=L+diag(ones(nRow,1));
P=P';
U=A; %
that is all
end
function
y=nma_ForwardSub(L,w)
%function y=nma_ForwardSub(L,w)
%
%Forward substitution that solves the lower
triangular system
%Ly=w for y
%
%
%INPUT:
%
L: an nxn L, lower tirangular
square matrix
%
w: vector of size n
%
%OUTPUT:
% y:
the solution to L*y=w
% By Nasser Abbasi
% HW 4, MATH 501, CSUF
%
if nargin ~= 2
error 'Only two inputs are required'
end
if
~(isnumeric(L)&isnumeric(w))
error 'input must be numeric'
end
[nRow,nCol]=size(w);
if nRow>1 &
nCol>1
error 'w must be a vector not a matrix'
end
[nRow,nCol]=size(L);
if nRow ~= nCol
error 'Matrix L must be square'
end
if length(w) ~=
nRow
error 'w length does not match L matrix dimension'
end
y=zeros(nRow,1);
y(1)=w(1)/L(1,1);
w=w(:);
for n=2:nRow
y(n)=(
w(n) - L(n,1:n-1)*y(1:n-1) ) / L(n,n);
end
end
function
x=nma_BackSub(U,v)
%function x=nma_BackSub(U,v)
%
%Backward substitution that solves the lower
triangular system
%Ux=v for x
%
%
%INPUT:
%
U: an nxn U, upper tirangular
square matrix
%
v: vector of size n
%
%OUTPUT:
% x:
the solution to U*x=v
% By Nasser Abbasi
% HW 4, MATH 501, CSUF
%
if nargin ~= 2
error 'Only two inputs are required'
end
if
~(isnumeric(U)&isnumeric(v))
error 'input must be numeric'
end
[nRow,nCol]=size(v);
if nRow>1 &
nCol>1
error 'v must be a vector not a matrix'
end
[nRow,nCol]=size(U);
if nRow ~= nCol
error 'Matrix U must be square'
end
if length(v) ~=
nRow
error 'v length does not match U matrix dimension'
end
x=zeros(nRow,1);
x(nRow)=v(nRow)/U(nRow,nRow);
v=v(:);
for n=(nRow-1):-1:1
x(n)=(v(n)-(U(n,n+1:end)*x(n+1:end))) / U(n,n);
end
end
function
x=nma_LinearSolve(A,b)
%function x=nma_LinearSolve(A,b)
%
% Solves Ax=b
%
%
%INPUT:
% A: an nxn square matrix
%
b: vector of size n
%
%OUTPUT:
% x:
the solution to A*x=b
% By Nasser Abbasi
% HW 4, MATH 501, CSUF
%
if nargin ~= 2
error 'Only two inputs are required'
end
if
~(isnumeric(A)&isnumeric(b))
error 'input must be numeric'
end
[nRow,nCol]=size(b);
if nRow>1 &
nCol>1
error 'b must be a vector'
end
[nRow,nCol]=size(A);
if nRow ~= nCol
error 'Matrix A must be square'
end
if length(b) ~=
nRow
error 'b
length does not match A matrix dimension'
end
b=b(:);
[L,U,P] =
nma_LU(A,0);
y = nma_ForwardSub(L,P*b);
x = nma_BackSub(U,y);
end