Programming Praxis – Solving Systems Of Linear Equations

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.

Tags: , , , , , , , , , , ,

5 Responses to “Programming Praxis – Solving Systems Of Linear Equations”

  1. programmingpraxis Says:

    Glad I made you think. . . .

    The non-stateful version does recursion down the main diagonal, splitting the matrix into four pieces: the current element, the vector to the right of the current element, the vector downward from the current element, and the remaining sub-matrix starting at the next element down the main diagonal, which is processed recursively. Depending on whether you are doing a lower or upper matrix, one of the two vectors becomes 0 in the output matrix, and there is a formula that doesn’t come to mind at the moment for the other vector. CLRS explains this fairly well before they switch from recursion to iteration for their pseudocode.

    The reason matrix algorithms are written in stateful pseudocode is because everybody other than Haskell programmers represents matrices using mutable two-dimensional arrays instead of lists of lists. If you’re going to process large matrices, say a million elements on a side, that’s almost a necessity; otherwise, if you have to generate and then garbage-collect million-by-million lists of lists, your program is going to take a while, and you will run out of space sooner than your imperative friend.

    Does Haskell have standard libraries for mutable matrices?


  2. Remco Niemeijer Says:

    Just did a quick search and didn’t find any mutable matrix libraries, though rolling your own using nested mutable arrays is fairly trivial.
    Here’s a very basic API I wrote while writing this reply:

    import Data.Array.IO
    import Data.Array.MArray
    type Matrix a = IOArray Int (IOArray Int a)
    fromList :: [[a]] -> IO (Matrix a)
    fromList xs = newListArray (1, length xs) =<< mapM (newListArray (1, length $ head xs)) xs
    toList :: Matrix a -> IO [[a]]
    toList m = mapM getElems =<< getElems m
    get :: Int -> Int -> Matrix a -> IO a
    get x y m = flip readArray x =<< readArray m y
    set :: Int -> Int -> a -> Matrix a -> IO ()
    set x y e m = (\a -> writeArray a x e) =<< readArray m y
    main = do m <- fromList [[1,2,3],[4,5,6]]
              print =<< get 3 1 m
              set 2 2 7 m
              print =<< toList m

    Pretty simple, but not nearly as nice as the pure approach since you’re now stuck in the IO monad, and having arrays means you can no longer do all the stuff I used in my solution like lazyness, recursion, using list-based library functions, etc. It may be efficient, but it also makes your program a lot more primitive.

    And of course I know that as soon as you start using big arrays the list-of-lists approach breaks down. But that’s the fun of toy problems: you get to use more creative solutions. Plus, since I’ve started Haskell I’ve come to see state as something that should be avoided if at all practical. It prevents certain classes of bugs and adding state to a pure algorithm, in the rare cases where it’s necessary, is a lot easier than going the other way round. Hence my choice for immutable nested lists for representing matrices.

  3. programmingpraxis Says:

    > Plus, since I’ve started Haskell I’ve come to see state
    > as something that should be avoided if at all practical.

    Agreed. But I draw the “practical” line in a different place than you. If I choose, I can use every tool that you can — recursion, higher-order functions, currying, laziness, list libraries. But equally I can ignore any or all of those tools and still get an elegant, beautiful program. Haskell comes with a strong prejudice, to its detriment.

  4. tinlyx Says:

    Just tried the code. But sometimes it gives errors:
    perm [[0.0,1.0,-1.0],[0.0,1.0,1.0],[0.0,-1.0,-1.0]]
    [[0.0,0.0,1.0],[1.0,0.0,0.0],[1.0,0.0,0.0],*** Exception: Prelude.head: empty list
    Any ideas to fix this?

  5. Remco Niemeijer Says:

    @tinlyx: The problem is that the first column of your matrix consists of zeroes. This results in a division by 0 in the elim function, at which point things go pear-shaped. Either remove the first column entirely (as well as one of the now redundant rows), or enter a non-zero value somewhere in the first column.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: