Solving Linear Equations

  

 

Input description: An matrix A, and an vector b, together representing m linear equations on n variables.

Problem description: What is the vector x such that ?

Discussion: The need to solve linear systems arises in an estimated 75% of all scientific computing problems [DB74].   For example, applying Kirchhoff's laws to analyze electric circuits generates a system of equations, the solution of which gives currents through each branch of the circuit. Analysis of the forces acting on a mechanical truss generates a similar set of equations. Even finding the point of intersection between two or more lines reduces to solving a (small) linear system.       

Not all systems of equations have solutions; consider the equations 2x+3y = 5 and 2x+3y = 6. Some systems of equations have multiple solutions; consider the equations 2x+3y=5 and 4x+6y=10.   Such degenerate systems of equations are called singular, and they can be recognized by testing whether the determinant of the coefficient matrix is zero.    

Solving linear systems is a problem of such scientific and commercial importance that excellent codes are readily available. There is likely no good reason to implement your own solver, even though the basic algorithm (Gaussian elimination) is one we learned in high school.     This is especially true if you are working with large systems.

Gaussian elimination is based on the fact that the solution to a system of linear equations is invariant under scaling (multiplying both sides by a constant; i.e. if x=y, then 2x=2y) and adding equations (i.e. the solution to the equations x=y and w=z is the same as the solution to x=y and x+w=y+z).   Gaussian elimination scales and adds equations so as to eliminate each variable from all but one equation, leaving the system in such a state that the solution can just be read off from the equations.

The time complexity of Gaussian elimination on an system of equations is , since for the ith variable we add a scaled copy of the n-term ith row to each of the n-1 other equations. On this problem, however, constants matter. Algorithms that only partially reduce the coefficient matrix and then backsubstitute to get the answer use 50% fewer floating-point operations than the naive algorithm.  

Issues to worry about include:

• Are roundoff errors and numerical stability affecting my solution? -     Implementing Gaussian elimination would be quite straightforward except for round-off errors, which accumulate with each row operation and can quickly wreak havoc on the solution, particularly with matrices that are almost singular.

  To eliminate the danger of numerical errors, it pays to substitute the solution back into each of the original equations and test how close they are to the desired value.   Iterative methods for solving linear systems refine initial solutions to obtain more accurate answers - good linear systems packages will include such routines.

  The key to minimizing roundoff errors in Gaussian elimination is selecting the right equations and variables to pivot on, and to scale the equations so as to eliminate large coefficients.     This is an art as much as a science, which is why you should use one of the many well-crafted library routines described below.

• Which routine in the library should I use? - Selecting the right code is also somewhat of an art. If you are taking your advice from this book, you should start with the general linear system solvers. Hopefully they will suffice for your needs. But search through the manual for more efficient procedures for solving special types of linear systems. If your matrix happens to be one of these special types, the solution time can reduce from cubic to quadratic or even linear.
• Is my system sparse? -   The key to recognizing that you might have a special-case linear system is establishing how many matrix elements you really need to describe A. If there are only a few nonzero elements, your matrix is sparse and you are in luck. If these few nonzero elements are clustered near the diagonal, your matrix is banded and you are in even more luck.   Algorithms for reducing the bandwidth of a matrix are discussed in Section .   Many other regular patterns of sparse matrices can also be exploited, so see the manual of your solver or a better book on numerical analysis for details.
• Will I be solving many systems with the same coefficient matrix? - In certain applications, such as least-squares curve fitting and differential equations, we have to solve repeatedly with different b vectors. For efficiency, we seek to preprocess A to make this easier.       The lower-upper or LU-decomposition of A creates lower- and upper-triangular matrices L and U such that . We can use this decomposition to solve , since

  

  This is efficient since backsubstitution solves a triangular system of equations in quadratic time. Solving and then gives the solution x using two steps instead of one step, once the LU-decomposition has been found in time.

The problem of solving linear systems is equivalent to that of matrix inversion, since , where is the identity matrix.     However, avoid it, since matrix inversion proves to be three times slower than Gaussian elimination. LU-decompositions prove useful in inverting matrices as well as computing determinants (see Section ).

Implementations: The library of choice for solving linear systems is apparently LAPACK, a descendant of LINPACK [DMBS79]. Both of these Fortran codes, as well as many others, are available from Netlib. See Section .        

Algorithm 533 [She78], Algorithm 576 [BS81], and Algorithm 578 [DNRT81] of the Collected Algorithms of the ACM are Fortran codes for Gaussian elimination.   Algorithm 533 is designed for sparse systems, algorithm 576 to minimize roundoff errors, and algorithm 578 to optimize virtual memory performance.   Algorithm 645 [NW86] is a Fortran code for testing matrix inversion programs. See Section for details on fetching these programs.

Numerical Recipes [PFTV86] provides routines for solving linear systems.   However, there is no compelling reason to use these ahead of the free codes described in this section.

C++ implementations of algorithms to solve linear equations and invert matrices are embedded in LEDA (see Section ).    

Notes: Good expositions on algorithms for Gaussian elimination and LU-decomposition include [CLR90] and a host of numerical analysis texts [BT92, CK80, SK93].

Parallel algorithms for linear systems are discussed in [Ort88]. Solving linear systems is one of relatively few problems where parallel architectures are widely used in practice.  

Matrix inversion, and hence linear systems solving, can be done in matrix multiplication time using Strassen's algorithm plus a reduction.   Good expositions on the equivalence of these problems include [AHU74, CLR90].

Certain types of nonsparse systems can be solved efficiently via special algorithms.   In particular, Toeplitz matrices are constructed so that all the elements along a particular diagonal are identical, and Vandermonde matrices   are defined by an n-element vector x where .

Related Problems: Matrix multiplication (see page ), determinant/permanent (see page ).