(The linear equation solver library is hosted at http://github.com/LeventErkok/conjugateGradient. Comments, bug reports, and patches are always welcome.)
Sparse matrix linear-equation solver, using the conjugate gradient algorithm. Note that the technique only applies to matrices that are symmetric and positive definite. See http://en.wikipedia.org/wiki/Conjugate_gradient_method for details.
The conjugate gradient method can handle very large sparse matrices, where direct methods (such as LU decomposition) are way too expensive to be useful in practice. Such large sparse matrices arise naturally in many engineering problems, such as in ASIC placement algorithms and when solving partial differential equations.
Here's an example usage, for the simple system:
4x + y = 1 x + 3y = 2
let a = SM (2, fromList [(0, SV (fromList [(0, 4), (1, 1)])), (1, SV (fromList [(0, 1), (1, 3)]))]) :: SM Double
let b = SV (fromList [(0, 1), (1, 2)]) :: SV Double
let g = mkStdGen 12345
let x = solveCG g a b
putStrLn $ showSolution 4 a b xA | x = b --------------+---------------- 4.0000 1.0000 | 0.0909 = 1.0000 1.0000 3.0000 | 0.6364 = 2.0000
- newtype SV a = SV (IntMap a)
- newtype SM a = SM (Int, IntMap (SV a))
- lookupSV :: Num a => Int -> SV a -> a
- lookupSM :: Num a => (Int, Int) -> SM a -> a
- addSV :: Num a => SV a -> SV a -> SV a
- subSV :: Num a => SV a -> SV a -> SV a
- dotSV :: Num a => SV a -> SV a -> a
- normSV :: RealFloat a => SV a -> a
- sMulSV :: Num a => a -> SV a -> SV a
- sMulSM :: Num a => a -> SM a -> SM a
- mulSMV :: Num a => SM a -> SV a -> SV a
- solveCG :: (RandomGen g, RealFloat a, Random a) => g -> SM a -> SV a -> SV a
- showSolution :: RealFloat a => Int -> SM a -> SV a -> SV a -> String
A sparse vector containing elements of type
a. Only the indices that contain non-
0 elements should be given
for efficiency purposes. (Nothing will break if you put in elements that are
0's, it's just not as efficient.)
A sparse matrix is essentially an int-map containing sparse row-vectors:
- The first element,
n, is the number of rows in the matrix, including those with all
- The matrix is implicitly assumed to be
nxn, indexed by keys
- When constructing a sparse-matrix, only put in rows that have a non-
0element in them for efficiency.
- Note that you have to give all the non-0 elements: Even though the matrix must be symmetric for the algorithm
to work, the matrix should contain all the non-
0elements, not just the upper (or the lower)-triangle.
- Make sure the keys of the int-map is a subset of
[0 .. n-1], both for the row-indices and the indices of the vectors representing the sparse-rows.
Norm of a sparse vector. (Square-root of its dot-product with itself.)
Multiply a sparse matrix (nxn) with a sparse vector (nx1), obtaining a sparse vector (nx1).
|:: (RandomGen g, RealFloat a, Random a)|
The seed for the random-number generator.
|-> SM a|
|-> SV a|
|-> SV a|
Conjugate Gradient Solver for the system
Ax=b. See: http://en.wikipedia.org/wiki/Conjugate_gradient_method.
NB. Assumptions on the input:
Amatrix is symmetric and positive definite.
- All non-
0rows are present. (Even if the input is assumed symmetric, all rows must be present.)
- The indices start from
0and go consecutively up-to
n-1. (Only non-
0value/row indices has to be present, of course.)
For efficiency reasons, we do not check that these properties hold of the input. (If these assumptions are violated, the algorithm will still produce a result, but not the one you expected!)
We perform either
10^6 iterations of the Conjugate-Gradient algorithm, or until the error
factor is less than
1e-10. The error factor is defined as the difference of the norm of
the current solution from the last one, as we go through the iterative solver. See
for a discussion on the convergence properties of this algorithm.
The solver can throw an error if it does not converge by
10^6 iterations. This is typically an indication
that the input matrix is not well formed, i.e., not symmetric positive-definite.
|:: RealFloat a|
Precision: Use this many digits after the decimal point.
|-> SM a|
|-> SV a|
|-> SV a|
Display a solution in a human-readable form. Needless to say, only use this method when the system is small enough to fit nicely on the screen.