Linear Algebraic Equation Solvers

Systems of linear algebraic equations occur in most fields of science and engineering, for example,

structural analysis (civil engineering)
heat conduction (mechanical and chemical engineering)
analysis of power grids (electrical engineering)
production planning (economics)
regression analysis (statistics)

due to, for example, the discretisation of ordinary or partial differential equations or integral equations.

A linear equation is an equation that can be expressed as,
a₁x₁ + a₂x₂ + … + a_nx_n = b

where a₀, a₁ … and b are constants.

A set of linear equations is called a system of linear equations or a linear system.

A linear system of n linear equations in n variables,

a_1,1x₁ + a_1,2x₂ + … a_1,nx_n = b₁
a_2,1x₁ + a_2,2x₂ + … a_2,nx_n = b₂
…
a_n,1x₁ + a_n,2x₂ + … a_n,nx_n = b_n

can be expressed in the explicit matrix form as

a_1,1    a_1,2    …  a_1,n
a_2,1    a_2,2    …  a_2,n
…
…
a_n,1  a_n,2 …  a_n,n

x₁
x₂
…
…
x_n

b₁
b₂
…
…
b_n

or more conveniently still, by matrix notation as

Ax = b

where A is an n × n matrix containing the values of i,j s and x and b are n element vectors storing the values of xi s and bi s.

If there are fewer equations than unknowns, M < N, or M = N but the equations are degenerate, (effectively fewer equations than unknowns), there is either no solution, or more than one solution vector x. In the latter case the solution space consists of a particular solution xp added to any linear combination of typically N − M vectors. The task of finding the solution space of A involves Singular Value Decomposition of matrix A.

If there are more equations than unknowns, M > N, there is, in general, no solution vector x to the linear system, and the system is said to be overdetermined. Frequently the best compromise solution is sought, the one that comes closest to satisfying all the equations simultaneously. If closeness is defined based on least squares, (minimising the sum of the squares of the differences of the left and right hand sides of the equations in the linear system), then an overdetermined linear problem reduces to a usually solvable linear problem called the Linear least-squares problem.

These problems are not discussed further.

There are two methods to solve a system of linear equations, direct or iterative.

Direct Methods

Direct Methods solve systems of linear equations by rearranging matrix A so that values of b can be substituted for x. They are finite methods; that is they will find a solution in a finite, predetermined number of steps.

Gaussian Elimination

One of the oldest and most useful methods of solving a system of linear equations Ax = b, is Gaussian elimination. The method uses the systematic subtraction of multiples of equations from other equations to eliminate variables equation by equation until the system is upper triangular. This triangularisation process is referred to as decomposition.

In an upper triangular system, the matrix has only zeros below the main diagonal; the first equation contains all the variables, the second equation is missing x1, the third equation is missing x1 and x2 ... and the last (nth) equation only contains xn. The last equation gives xn. This is then back substituted into the next to last equation to find xn-1. This process of back‐substitution is repeated until the first equation is reached and x1, is solved using the calculated values of x2 ... xn.

Gaussian elimination and back‐substitution are easily programmed and a C++ implementation, based on the FORTRAN code in Chemical Engineering Volume 3 2nd Edition page 317 and Pascal programs Gauss and Solve, Numerical Methods for Engineers and Computer Scientists pages 77 and 83, is provided here, Gaussian Elimination Solver.

Gauss‐Jordan Elimination

Another direct method to solve the system of equations is Gauss‐Jordan elimination. This involves complete elimination, in which zeros are produced for all off‐diagonal elements, both above and below the diagonal.

For large systems the multiplication and division work for Gaussian triangularisation with back substitution tends towards a factor of two‐thirds of the work for Gauss‐Jordan complete reduction, [Numerical Methods for Engineers and Computer Scientists]. Therefore, I haven't implemented Gauss‐Jordan elimination as a routine. A C++ routine is provided in Numerical Recipes in C++ page 42.

LU Decomposition

Both Gaussian elimination and Gauss‐Jordan elimination have the disadvantage that all right‐hand sides must be known in advance. If they are not known LU decomposition can be used.

The basis of LU decomposition is to write matrix A as the product of two matrices

L ⋅ U = A(eqn 1)

where L is lower triangular (has elements only on the diagonal and below) and U is upper triangular (has elements only on the diagonal and above).

Equation 1 can be decomposed to solve the linear set
A ⋅ x = (L ⋅ U) ⋅ x = L ⋅ (U ⋅ x) = b (eqn 2)

by first solving for vector y where
L ⋅ y = b(eqn 3)

and then solving
U ⋅ x = y(eqn 4)

The solution of a triangular set of equations is quite trivial, as used in the back substitution in Gauss elimination and Gauss‐Jordan solvers. Equation 3 is solved by forward substitution and equation 4 is solved by back substitution as is done for Gauss elimination and Gauss‐Jordan solvers.

I have not implemented an LU decomposition algorithm, as I have already implemented Gaussian elimination with back substitution as it is the algorithm provided in most textbooks and so was the algorithm I came across when I was first trying to solve a system of linear equations. And going forward I expect to be using iterative methods to solve a set of linear equations as iterative solvers overcome the following problems with Direct Methods. A C++ routine is provided in Numerical Recipes in C++ page 49.

Problems with Direct Methods

The large number of calculations required for elimination and back substitution combined with the limited fixed precision of the numbers stored in a computer mean that it can be difficult for direct methods to obtain a solution with an accuracy comparable with the precision of the data. These problems are described in Numerical Problems with Direct Methods.

In many problems many of the coefficients off the principle diagonal will be zero, this means that matrix will use a lot more memory than necessary, and although since mid-2000 when P.C.s began having at least 1GB of RAM this probably is no longer a problem, all these zero values will still be checked by the algorithm algorithm, impacting performance.

For problems where the non-zero coefficients lie on and immediately above and below the main diagonal the matrix is termed a banded matrix. If only the main diagonal has non-zero values the matrix is called a diagonal matrix; if the diagonals immediately above and below the main diagonal have non-zero values the matrix is called a tridiagonal matrix. If the two diagonals immediately above and below the main diagonal have non-zero values the matrix is called pentadigaonal.

Specialised algorithms tailored to the structure of the matrix reduce the memory and calculations required to store and solve the system, for example the Pascal Procedures TRIDIAG and PENTDIAG in Numerical Methods for Engineers and Computer Scientists.

These problems are solved using Iterative Methods to solve the system of equations.

Iterative Methods

Iterative methods solve the system of linear equations, Ax = b, by generating a series of increasingly better approximations to the solution vector, x. Iterative methods are a better option for solving larger matrices for the following reasons,

The matrix is not modified (e.g. no decomposition), avoiding the problem of rounding errors. And the matrix can be used to solve multiple problem cases
Each iteration corrects for any rounding errors
Can solve matrices stored in Sparse Format

The theory of iterative methods for solving the equations is straightforward. Vector x is the exact solution of the linear system Ax = b. However we don't know x, we only know a slightly wrong solution x + δx where δx is error in our guess. When multiplied by A our slightly wrong value of x give a slightly wrong value of the desired result, b namely
A ⋅ (x + δx) = b + δb

By implementing an algorithm to calculate the value of δx that minimises δb we can find the values for x. As listed in the pros of iterative methods above, this means that the matrix is not manipulated and a general routine (common ones listed below) is used to calculate the values of x.

A number of algorithms have been developed to calculate the new values for x that will minimise δb. Common algorithms are,

the Jacobi method, Gauss‐Seidel method, and Succesive Over‐Relaxation (SOR) methods are Stationary methods, which are older, simpler to understand but converge too slowly to the solution to be effective
Nonstationary methods, which are harder to understand, but can be highly effective, include
the conjugate gradient method which has a higher rate of convergence and is a special case of the ...
bi‐conjugate gradient method. An implementation of the algorithm is provided. This is the most advanced algorithm provided in Numerical Recipes in C++
the IML++ (Iterative Methods Library) maintained by the National Institute of Standards and Technology, Oak Ridge National Laboratory, is a C++ templated library of modern iterative methods for solving both symmetric and non-symmetric linear systems of equations and has algorithms for Richardson Iteration, Chebyshev Iteration, Conjugate Gradiant (CG), Conjugate Gradiant Squared (CGS), BiConjugate Gradiant (BiCG), BiConjugate Gradiant Stablilised (BiCGSTAB), Generalised Minimum Residual (GMRES) and Quasi-Minimal Residual Without Lookahead (QMR)

IML++

IML++ uses the template facilities of C++ and the SparseLib++ Sparse Matrix Class Library by the National Institute of Standards and Technology, and at this time was too complex for me to implement for my needs.
PDF and PostScript versions of Templates: Building Blocks for the Solution of Iterative Systems can be found here. (The List of Symbols in this document was helpful)
GMRES combines the conjugate gradient method with a preconditioning technique, and Numerical Recipes in C++ implies is an improved solver compared to the bi-conjugate gradient method. A PDF of the original paper presented by Youcef Saad and Martin H. Schultz in 1986 can be found here

Sparse Matrices

In science and engineering many systems of equations, for example the discretisation of partial differential equations, the calculation of flows between equipment, etc. often result in the creation of spare systems of equations, that is for each equation most of the elements are zero. In a Full Format matrix, all these zeros are stored and manipulated during decomposition, which uses memory and requires unnecessary computes. This resulted in the development of Sparse Matrix Formats, where only the non zero values are stored along with their location, as described in the PDF of the SparseLib++ user's guide. (A Postscript version is also available on the SparseLib++ website)

While the ability to manufacture cheap memory modules and storage and 64bit CPUs has generally made memory size requirements less of an issue, but not manipulating zero elements still improves performance.