Posts Tagged ‘linear’

Programming Praxis – Solving Systems Of Linear Equations

July 21, 2010

In yesterday’s Programming Praxis exercise our task is to implement some more matrix-related functions, specifically LU and LUP decomposition and a solver for systems of linear equations. The provided Scheme solution comes in at 69 lines. Let’s see what we can do about that.

Some imports:

import Control.Arrow
import Data.List
import qualified Data.List.Key as K

We’re going to need the matrix multiplication function we defined last time.

mult :: Num a => [[a]] -> [[a]] -> [[a]]
mult a b = [map (sum . zipWith (*) r) $ transpose b | r <- a]

Since we will be using Gauss elimination, we need a way to eliminate a single row.

elim :: Fractional a => [a] -> [a] -> [a]
elim ~(x:xs) ~(y:ys) = zipWith (-) ys $ map (y / x *) xs

Also, we need to construct identity matrices.

identity :: Num a => [[a]] -> [[a]]
identity = zipWith (zipWith const) (iterate (0: ) (1 : repeat 0))

With that out of the way, we can defined LU decomposition.

lu :: Fractional a => [[a]] -> ([[a]], [[a]])
lu = unzip . map unzip . f where
    f []      = []
    f ~(x:xs) = zip (1 : repeat 0) x :
                zipWith (:) (map (\(y:_) -> (y / head x, 0)) xs)
                            (f $ map (elim x) xs)

LUP decomposition is the same as LU decomposition, except that the matrix is permuted so that the pivot at each step is maximized. Therefore, we need a way to compute the correct permutation matrix.

perm :: (Fractional a, Ord a) => [[a]] -> [[a]]
perm m = f $ zip (identity m) m where
    f [] = []
    f xs = a : f (map (second $ elim b) $ delete (a,b) xs)
           where (a,b) = K.maximum (abs . head . snd) xs

As mentioned, once we have the permutation matrix, LUP decomposition is trivial.

lup :: (Fractional a, Ord a) => [[a]] -> ([[a]], [[a]], [[a]])
lup xs = (perm xs, l, u) where (l,u) = lu $ mult (perm xs) xs

And finally, the function to solve systems of equations. On a related note, I’d really appreciate it if algorithm descriptions came with pseudocode that wasn’t stateful. This one took some thinking to convert to a functional style.

lupsolve :: (Fractional a, Ord a) => [[a]] -> [a] -> [a]
lupsolve a b = f y u where
    (p,l,u) = lup a
    y = foldl (\x (l', pb') -> x ++ [pb' - sum (zipWith (*) x l')])
              [] (zip l (concat . mult p $ map return b))
    f _ [] = []
    f ~(y':ys) ~((r:rs):us) = (y' - sum (zipWith (*) rs z)) / r : z
        where z = (f ys $ map tail us)

As usual, a test to see if everything’s working correctly:

main :: IO ()
main = do print $ lu [[ 2, 3, 1, 5]
                     ,[ 6,13, 5,19]
                     ,[ 2,19,10,23]
                     ,[ 4,10,11,31]]
          print $ lup [[ 2, 0,   2,3/5]
                      ,[ 3, 3,   4, -2]
                      ,[ 5, 5,   4,  2]
                      ,[-1,-2,17/5, -1]]
          print $ lupsolve [[1,2,0],[3,5,4],[5,6,3]]
                           [1/10, 25/2, 103/10]

Yep (for more accurate results, use Ratio Ints instead of the default Floats). And at 20 lines, that’s less than a third of the Scheme solution. Not too bad.