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.