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.
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.