Introduction
Theory
HOWTO
Error Analysis
Examples
Questions
Applications in Engineering
Matlab
Maple
BLAS
Introduction
You can solve a system of linear equation using Gaussian elimination and backward
substitution, however there are two drawbacks of using such a technique:
the Gaussian elimination component of the algorithm is slow, O(n^{3}),
relative to the backward substitution component, O(n^{2}), and
the Gaussian elimination is subject to numeric problems.
In this topic, we will see how we can avoid doing Gaussian elimination
repeatedly and thereby perform the most expensive component of the algorithm
only once, using the faster O(n^{2}) operations to solve numerous
different problems.
Background
Useful background for this topic includes:
 Gaussian elimination and backwardsubstitution.
References
Interactive Maplet
A PLU Decomposition Tutorial
Theory
Solving Systems of Linear Equations
Consider the following system of n equations and n unknowns.
m_{1,1} x_{1} + m_{1,2} x_{2} + ⋅⋅⋅ m_{1,n} x_{n} = b_{1}
m_{2,1} x_{1} + m_{2,2} x_{2} + ⋅⋅⋅ m_{2,n} x_{n} = b_{2}
⋅
⋅
⋅
m_{n,1} x_{1} + m_{n,2} x_{2} + ⋅⋅⋅ m_{n,n} x_{n} = b_{n}
In linear algebra, you saw how you could convert this into a matrixvector problem
M x = b
where M = (m_{i,j}) is the matrix of coefficients,
x = (x_{i})^{T} is a column vector of unknowns
and
b = (b_{i})^{T} is a column vector of known values.
To solve such a system, you would create the augmented matrix (Mb) and
proceed with applying Gaussian elimination and backward substitution.
If you do not recall forward or backward substitution, please see the corresponding
HOWTO pages (follow the links).
Problems with Gaussian Elimination
There are two problems with Gaussian elimination:
1. Speed
Even if you do an optimal implementation of Gaussian elimination and
backward substitution for solving a system of n linear equations, you
require the following number of operations:
Table 1. Operations due to Gaussian elimination.
 n^{2}/2 − n/2 divisions and negations
 n^{3}/3 − n/3 multiplications and additions
Table 2. Operations due to backward substitution.
 n^{2}/2 − n/2 multiplications and subtractions
 n divisions
Thus, most of the time is spent in Gaussian elimination: for example, for a 1000 × 1000
matrix M would require 333832500 operations to do the Gaussian elimination, but only
500500 operations to do the backward substitution. This is summarized by stating that
Gaussian elimination is O(n^{3}) while backward substitution (and forward
substitution, for that matter) are O(n^{2}).
2. Numeric Instability
Try solving the following system of equations using 3 decimal digits of
precision:
Performing Gaussian elimination, we add 21.4/0.353 = 60.6 times row 1 onto row 2 to get:
which is now in echelon form, so we may perform backward substitution:
If we now try to solve this using Matlab, we get a very different answer, namely
x = (0.01280, 0.12536)^{T}. Our solution for x_{2} may
be reasonable (a relative error of 0.0029), but the relative error of 0 for any number is 1.0, so this is a very
poor approximation.
Try it again, only now swap the two rows:
and apply Gaussian elimination:
Applying backward substitution yields x_{2} = 0.126 and x_{1} = 0.0133.
Thus, the relative error of both of these numbers are 0.0051 and 0.039, respectively: a significant
overall improvement.
What happened here was that by adding a large multiple of the first row
onto the second, because of the finite precision of our approximation, we
ended up with a matrix of the form
and when we apply backward substitution, we get x_{2} = b_{1}/m_{12} and therefore
x_{1} = 0.
Solving Problem 1
Suppose we could write the matrix M as a product of two matrices:
M = L U where L is a lowertriangular
matrix and U is an uppertriangular matrix. In this case, we
could then rewrite the matrixvector problem M x as
L U x = b
Thus, let us define y = U x. In this case, the equation may
be written as:
L y = b
In this case, L is a lowertriangular matrix and therefore we do not
have to perform Gaussian elimination to solve for y and instead we
only use forward substitution. Similarly, once
we have found y, we may now solve:
U x = y
for x. Again, U is uppertriangular, so we again do not have
to perform Gaussian elimination, we need only use backward substitution.
From before we may now calculate the number of required operations to
perform both forward and backward substitution:
Table 3. Operations due to forward and backward substitution.
 n^{2} − n multiplications and subtractions
 2n divisions
Thus, if we can find L and U, we could easily solve
M x = b_{i} for N different known vectors
b_{i} using only Nn^{2} + Nn operations.
Fortunately, we are already half way to writing M as a product L U:
when we do Gaussian elimination, we are left with an uppertriangular matrix U.
By looking at how we got to finding U, we will see that we can, simultaneously,
find the corresponding lowertriangular matrix L.
First, recall that adding a multiple of one row onto another row of a matrix
is termed an elementary operation and may be represented by an
elementary matrix. If you are adding c times row ν onto
row μ, the corresponding elementary matrix E = (e_{i j})
is an n × n
identity matrix with the entry e_{μ ν} = c.
Consider the following example of a Gaussian elimination:
Adding 0.5 times row 1 onto row 2 yields:
Adding 0.2 times row 1 onto row 3 yields:
Finally, adding 0.1 times row 2 onto row 3 yields
Now, the inverse of the elementary matrices defined above are simply the
identity matrix with the negative of the offdiagonal entry, so by multiplying
on the left by each of the three inverses, we get:
If you multiply out the three elementary matrices in the given order, you will find that
the product is simply:
Thus, we have found an LU decomposition of the matrix M. It should
be noted that there are many LU decompositions.
Thus, if we wished to solve M x = b where b = (19., 1.5, 28.6)^{T},
we would apply forward substitution to solve L y = b to get that y = (19, 8, 24)^{T}
and then solve U x = y using backward substitution to find that x = (2, 2, 3)^{T}.
Therefore, we only have to solve the LU decomposition once but we can solve
hundreds of systems of equations of the form M x = b_{i} at
a much cheaper cost.
Solving Problem 2
In the above example, we solved the first problem of time: solve once for L and U and
then use these to solve for many solution vectors using the much less expensive forward and backward
substitutions. However, we did not solve the second problem, that of the numerical instability.
If we examine the problem, we note that it occurs because we are adding a large multiple of
one row onto another. Recall that with Gaussian elimination, if you found a zero on the diagonal,
you could swap the row containing the diagonal with a row beneath it which did not contain
a zero in that column. We will do the same thing here, but instead, with each step, we will first
swap the row with largest element (in absolute value) on or below the diagonal with the row
containing the diagonal. The only thing left to do is to keep track of the row swaps.
This may be done by factoring the matrix M as a triplet P L U where P is
a permutation matrix. As the inverse of a permutation matrix is simply its transpose, if
we try to solve P L U x = b, it is equivalent to solving L U x = P^{T}b.
Without proof, the HOWTO for this section covers how to find the three matrices
P^{T}, L and U in a reasonably efficient and straightforward manner.
HOWTO
Problem
Given the n × n matrix M, find a decomposition
M = PLU where
P is a permutation matrix, L is a
lowertriangular matrix, and U is an uppertriangular matrix. Next, given
a column vector b, solve the system of linear equations defined
by M x = b for the vector of unknowns x.
Assumptions
We will assume that M is invertible. If it
is not, one of the steps will fail and we will not be able
to find such a decomposition.
PLU Decomposition
Start with three candidate matrices:
 U = M
 L = 0_{n, n}
 P^{T} = I_{n}
where 0_{m, n} is the
m × n zero matrix and I_{n}
is the n × n identity matrix.
Then, for i = 1, 2, ..., n  1 do:
Find the row j with the largest entry in absolute value on
or below the diagonal of the ith row and swap rows
i and j in all three matrices, P^{T}, L, and U.
If this maximum entry is zero, terminate this loop and indicate that
the matrix is singular, that is, noninvertible.
For j = i + 1, ..., n, calculate the scalar
value s = u_{j, i}/u_{i, i}.
Add s times row i onto row j in the matrix U
and set the entry l_{j, i} = s.

Having iterated from i = 1, 2, ..., n, finish by
adding the identity matrix onto L = L + I_{n}.
These are the P^{T}, L, and U matrices of
the PLU decomposition of M.
Solving with Forward and Backward Substitution
Calculate P^{T}⋅b and use forward substitution to solve
for y in L⋅y = P^{T}⋅b.
Then, having found y, use backward substitution to solve for
x in U⋅x = y.
Error Analysis
We will look at the error analysis for using PLU decomposition
when we look at the topic of illconditioned matrices.
Examples
Example 1
Find the PLU decomposition of the matrix:
Doing this for the given matrix, we first set up the three candidate matrices, a 4 × 4 identity
matrix, a 4 × 4 zero matrix, and the matrix M. These are the candidate matrices for
the three matrices P^{T}, L, and U, respectively.
Starting with column 1, we swap rows 1 and 2 to place 12 on the diagonal:
Add 0.5 times row 1 onto row 2, add 0.2 times row 1 onto row 3, and
set λ_{2,1} = 0.5 and
λ_{3,1} = 0.2. The entry u_{4,1} is
already zero.
Moving onto column 2, swap rows 2 and 3 to place the largest value on
the diagonal:
The only nonzero entry below the diagonal of U in column 2 is
u_{4, 2},
and we eliminate this by adding 0.1 times row 2 onto row 4 and set
λ_{4, 2} = 0.1:
Moving onto column 3, we swap rows 3 and 4 to bring the largest entry on
or below the diagonal of column 4 onto the diagonal:
To complete the procedure, we add 0.25 times row 3 onto row 4 and set
λ_{4, 3} = 0.25:
There is nothing to do in the 4th column, so we are finished. We
add the identity matrix back to get the following three matrices:
Thus we have found the three matrices P^{T}, L, and U, respectively.
If we take the transpose of the first matrix and calculate the product
PLU, we get:
which is, of course, equal to M.
Example 2
Use the PLU decomposition from Example 1 to solve Mx = b
for x when b = (0.75, 22, 3.6, 0.2)^{T}.
First, we replace M with PLU, that is, PLUx = b.
By multiplying both sides on the left by the inverse of P, that is,
P^{T},
we get LUx = P^{T}b.
In our example, we calculate P^{T}b:
Next, writing our equation as L(Ux) = P^{T}b, we let Ux = y, and we solve
Ly = P^{T}b. We may do
this by writing the augmented matrix:
and using forward substitution:
 y_{1} = 22.
 y_{2} = 3.6 − 0.2 y_{1} = 8
 y_{3} = 0.2 − 0 y_{1} − 0.1 y_{2} = 1
 y_{4} = 0.75 + 0.5 y_{1} − 0 y_{2} − 0.25 y_{3} = 10
Finally, having found y, we go back to solve for x by
creating the augmented matrix of U and y:
and use backward substitution to find x:
 x_{4} = 10/10 = 1
 x_{3} = (1 − 1 x_{4}) / 15 = 0
 x_{2} = (8 + 2 x_{3} − 2 x_{4}) / 10 = 10/10 = 1
 x_{1} = (22 − 2 x_{2} − 1 x_{3} − 0 x_{4}) / 12 = 24/12 = 2
Thus, the solution to Mx = b is the vector
and to confirm, calculating Mx, we get the vector b.
You may wish to compare the amount of work done in this example
as compared to Example 1.
Questions
Question 1
Find the PLU decomposition of the matrix
Answer:
L = [1 0 0 0; 0.2 1 0 0; 0.1 0 1 0; 0.2 0.5 0.1 1]
U = [10 3 2 3; 0 8 1 2; 0 0 12 3; 0 0 0 8]
Pt = [1 0 0 0; 0 0 1 0; 0 0 0 1; 0 1 0 0]
Question 2
Using the PLU decomposition of M from Question 1, use
backward and forward substitution
to solve Mx = b
where b = (14.0, 12.9, 5.8, 19.4)^{T}.
Answer: x = (1, 0, 1, 2)^{T}.
Question 3
Find the PLU decomposition of the matrix
Answer:
U = [6 2 1; 0 8 3; 0 0 5]
L = [1 0 0; 0.2 1 0; 0.1 0.3 1]
Pt = [0 0 1; 1 0 0; 0 1 0]
Question 4
Using the PLU decomposition of M from Question 3, use
backward and forward substitution
to solve Mx = b
where b = (16.4, 15.4, 12.0)^{T}.
Answer: x = (2, 1, 2)^{T}.
Applications to Engineering
Given a circuit made up of linear elements (resistors, inductors,
capacitors, and voltage and current sources), when doing nodal analysis
on such a circuit, you must solve a system of linear equations.
Given such a circuit, you may be interested as to the state of the
circuit given various input currents. Instead of solving the system
of linear equations to find the voltages for each new set of input currents,
you can find the PLU decomposition once and then simply use
forward and backward substitution to find the voltages for a given
set of input currents.
Given a system of linear equations which describe the zero input
response of a linear system described by a differential equation
(more on this in ECE 342), the initial conditions only affect the
solution vector while the matrix depends on the system. Thus, you may
solve for various initial conditions using only forward and backward
substitution.
Nodal Analysis (Simple)
Consider the circuit shown in Figure 1.
Figure 1. A simple circuit.
In ECE 100 you learn that we may apply nodal analysis by
introducing nodes which separate the circuit elements and
grounding one of the nodes. This is shown in Figure 2.
Figure 2. Introducing a node.
Considering the one reminding node in this circuit, we apply Kirkhoff's current
law (KCL), as is shown in Figure 3.
Figure 3. Applying KCL to the circuit.
This gives us the equation:
−10 mA + 1/(200 Ω) (v_{1} − 0) = 0
This may be rearranged into the linear equation:
1/(200 Ω) v_{1} = 10 mA
which we may solve for v_{1} = 2 V.
Nodal Analysis (More Interesting)
As a more interesting circuit, consider the circuit shown in
Figure 4.
Figure 4. A more interesting circuit.
Five nodes separate the circuit elements, one of which is
grounded, as is shown in Figure 5.
Figure 5. Introducing nodes.
Then, applying KCL to Node 1, as is shown in Figure 6, we get the
equation:
10 mA + 1/(120 Ω)(v_{1} − 0) + 1/(240 Ω)(v_{1} − v_{2}) = 0
Figure 6. Applying KCL to the circuit.
By removing the units and with appropriate scaling, if we apply this processing to the remaining nodes, we get the four equations:
0.01 + 1/120 (v_{1} − 0) + 1/240 (v_{1} − v_{2}) = 0
1/320 (v_{2} − 0) + 1/240 (v_{2} − v_{1}) + 1/180 (v_{2} − v_{3}) + 1/200 (v_{2} − v_{4}) = 0
1/160 (v_{3} − 0) + 1/180 (v_{3} − v_{2}) + 1/360 (v_{3} − v_{4}) = 0
0.01 + 1/200 (v_{4} − v_{2}) + 1/360 (v_{4} − v_{3}) = 0
By expanding these out, we have a system of four linear equations in four unknowns, which may therefore
be written as:
Solving this system gives us the voltage vector v = (0.619, 0.544, 0.556, 1.83)^{T}.
Using this information, we can, for example, determine the voltage flowing between two nodes, the
corresponding energy losses, etc.
Notice that the current only affects the righthand vector, and therefore, we could change the
currents, but continue to use the same conductance matrix. Thus, PLU decomposition would be
more efficient than applying Gaussian elimination.
Matlab
The follow Matlab code finds the PLU decomposition of the
matrix M:
M = [6 1 3.25 10.25; 12 2 1 0; 2.4 10.4 1.8 2; 0 1 14.8 1.2]
[L U Pt] = lu( M )
Here, the lower triangular matrix is assigned to the first output, L,
upper triangular matrix is assigned to the second, U, and the
transpose of the permutation matrix is assigned to the third, Pt.
To do this manually, we would do the following:
M = [6 1 3.25 10.25; 12 2 1 0; 2.4 10.4 1.8 2; 0 1 14.8 1.2]
n = length( M );
L = zeros( n, n );
U = M;
Pt = eye( n, n );
for i=1:(n  1)
% find the row with the maximum entry in the ith column on or below the diagonal
[tmp, j] = max( abs( U(i:n, i) ) ) ;
j = j + (i  1); % account for offset
if ( tmp == 0 )
error( 'the matrix is singular (non invertible)' );
end
% swap rows 'i' and 'j'
U ([i, j], :) = U ([j, i], :);
Pt([i, j], :) = Pt([j, i], :);
L ([i, j], :) = L ([j, i], :);
for j=(i + 1):n
s = U(j, i)/U(i, i); % calculate multiplier
U(j, :) = U(j, :) + s*U(i, :); % add the multiplier times the ith row onto the jth
L(j, i) = s; % store the negated multiplier
end
end
L = L + eye( n, n ); % add on the identity matrix
To solve the system M x = b where
b = [0.75 22 3.6 0.2]'
let Matlab use forward and backward substitution as follows:
y = L \ (Pt * b)
x = U \ y
Short Cut
If you simply wish to solve for x, just enter:
x = M \ b
Maple
The following commands in Maple finds the PLU decomposition
of the given matrix M:
M := <<6, 12, 2.4, 0><1, 2, 10.4, 1><3.25, 1, 1.8, 14.8><10.25, 0, 2, 1.2>>;
LinearAlgebra[LUDecomposition]( M );
b = <0.75, 22, 3.6, 0.2>
LinearAlgebra[LinearSolve]( M, b );
For more help on either of these routines or on the LinearAlgebra
package, enter:
?LinearAlgebra
?LinearAlgebra,LUDecomposition
?LinearAlgebra,LinearSolve
BLAS
The basic linear algebra subroutines (BLASs) for
performing a PLU decomposition together with solving the result are the two functions:
int gsl_linalg_LU_decomp( gsl_matrix *M, gsl_permutation *Pt, int );
int gsl_linalg_LU_solve( gsl_matrix *M, gsl_permutation *Pt, gsl_vector *b, gsl_vector *x );
An example of the use of these functions which finds the
coefficients of the polynomial interpolating the points
(1, sin(1)), (2, sin(2)), ..., (n, sin(n)) is provided
here in both a version which performs the operations but also
prints the matrices and vectors plu.print.c and
a version which simply performs the operations
plu.c.
To compile and execute this code, you must have the GNU Scientific Library
installed on your system, either in cygwin or in Linux. The example below looks at how you compile this code
in Cygwin:
$ gcc Wall I/usr/include plu.print.c lgsl lgslcblas lm
$ ./a.exe
You may view the output.
Reference: GNU Scientific Library,
13.1 LU Decomposition.