Posts Tagged ‘matrix’

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.

Programming Praxis – Matrix Operations

June 22, 2010

In today’s Programming Praxis exercise our task is to implement four common matrix operations. The provided Scheme solution has 44 lines. Let’s see what we can do to reduce that a bit.

There are many ways to represent matrices. The most common one in Haskell (if not the most efficient if you need random access) is to use a list of lists, which is what we’ll be using here.

For addition, all we need to do is add the corresponding numbers in the two matrices together.

add :: Num a => [[a]] -> [[a]] -> [[a]]
add = zipWith $ zipWith (+)

Scaling (multiplying by a constant) is even simpler: just apply the multiplication to every element.

scale :: Num a => a -> [[a]] -> [[a]]
scale = map . map . (*)

For transposition I could’ve just used the definition from Data.List, but that would be cheating. Instead, I figured I would just use a right fold.

transpose :: [[a]] -> [[a]]
transpose [] = []
transpose xs = foldr (zipWith (:)) (repeat []) xs

Multiplying two matrices requires multiplying rows and columns. Since nested lists don’t have easy access to columns, we can use the transpose function we just defined to fix that.

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

And of course we need to test if everything works correctly.

main :: IO ()
main = do let m = [[1..3],[4..6]]
          let n = [[1..4], [2..5], [3..6]]
          print $ add m [[2..4],[3..5]] == [[3,5,7], [7,9,11]]
          print $ scale 2 m == [[2,4,6], [8,10,12]]
          print $ mult m n == [[14,20..32], [32,47..77]]
          print $ transpose m == [[1,4], [2,5], [3,6]]

Yep. And with a 89% reduction in line count, that’s good enough for me.