next up previous contents index CD CD Algorithms
Next: Bandwidth Reduction Up: Numerical Problems Previous: Numerical Problems

Solving Linear Equations



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

Problem description: What is the vector x such that tex2html_wrap_inline26852 ?

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 tex2html_wrap_inline26874 system of equations is tex2html_wrap_inline26876 , 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:

The problem of solving linear systems is equivalent to that of matrix inversion, since tex2html_wrap_inline26894 , where tex2html_wrap_inline26896 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 gif).

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 gif.        

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 gif 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 tex2html_wrap_inline26898 algorithms to solve linear equations and invert matrices are embedded in LEDA (see Section gif).    

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 tex2html_wrap_inline26900 .

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

next up previous contents index CD CD Algorithms
Next: Bandwidth Reduction Up: Numerical Problems Previous: Numerical Problems

Mon Jun 2 23:33:50 EDT 1997