## Solving Systems Of Linear Equations

### July 20, 2010

We are following the matrix operations described by Cormen, Leiserson, Rivest and Stein in their textbook Introduction to Algorithms. Here is LU decomposition:

```(define (lu-decomposition a)   (let* ((n (matrix-rows a))          (l (make-matrix n n 0))          (u (make-matrix n n 0)))     (for (i n) (matrix-set! l i i 1))     (for (k n)       (matrix-set! u k k (matrix-ref a k k))       (for (i (+ k 1) n)         (matrix-set! l i k (/ (matrix-ref a i k) (matrix-ref u k k)))         (matrix-set! u k i (matrix-ref a k i)))       (for (i (+ k 1) n)         (for (j (+ k 1) n)           (matrix-set! a i j             (- (matrix-ref a i j)                (* (matrix-ref l i k)                   (matrix-ref u k j)))))))     (values l u)))```

Since the body of the function only processes the “interesting” elements of L and U, the two matrices are initialized to 0, and the first `for` loop, on `i`, sets the elements of the main diagonal of L to 1. Then the outermost of the three nested loops, on `k`, sets the pivot, the middle loop, on `i`, fills in the elements of the pivot row and pivot column, and the innermost of the three loops, on `j`, resets the original matrix for the next turn through the loop. Notice that the input matrix is modified during the decomposition. Here’s an example:

```> (define a     #( #( 2  3  1  5 )        #( 6 13 5 19 )        #( 2 19 10 23 )        #( 4 10 11 31 ))) > (lu-decomposition a) #(#(1 0 0 0) #(3 1 0 0) #(1 4 1 0) #(2 1 7 1)) #(#(2 3 1 5) #(0 4 2 4) #(0 0 1 2) #(0 0 0 3))```

LUP decomposition is similar, but instead of returning two matrices half-filled with zeros, we return the original matrix mutated with both half-matrices occupying their own portions of the matrix; the permutation matrix is returned as a vector of offsets:

```(define (lup-decomposition a)   (let* ((n (matrix-rows a)) (pi (make-vector n 0)))     (for (i n) (vector-set! pi i i))     (for (k n)       (let ((p 0) (k-prime 0))         (for (i k n)           (let ((x (abs (matrix-ref a i k))))             (when (< p x) (set! p x) (set! k-prime i))))         (when (zero? p)           (error 'lup-decomposition "singular matrix"))         (vector-swap! pi k k-prime)         (for (i n) (matrix-swap! a k i k-prime i))         (for (i (+ k 1) n)           (matrix-set! a i k                        (/ (matrix-ref a i k)                           (matrix-ref a k k)))           (for (j (+ k 1) n)             (matrix-set! a i j                          (- (matrix-ref a i j)                             (* (matrix-ref a i k)                                (matrix-ref a k j))))))))     (values a pi)))```

`Lup-decomposition` calls two helper functions:

```(define (vector-swap! v a b)   (let ((t (vector-ref v a)))     (vector-set! v a (vector-ref v b))     (vector-set! v b t)))```

```(define (matrix-swap! m ar ac br bc)   (let ((t (matrix-ref m ar ac)))     (matrix-set! m ar ac (matrix-ref m br bc))     (matrix-set! m br bc t)))```

Here’s an example:

```> (define a     #( #(  2    0    2  3/5 )        #(  3    3    4   -2 )        #(  5    5    4    2 )        #( -1   -2 17/5   -1 ))) > (lup-decomposition a) #(#(5 5 4 2) #(2/5 -2 2/5 -1/5) #(-1/5 1/2 4 -1/2) #(3/5 0 2/5 -3)) #(2 0 3 1)```

Since `lup-decomposition` returns both the L and U matrices intertwined in a single matrix, it is convenient to define functions that untangle them:

```(define (lower lu i j)   (cond ((< i j) 0) ((= i j) 1)     (else (matrix-ref lu i j))))```

```(define (upper lu i j)   (if (< j i) 0     (matrix-ref lu i j)))```

Then the solver is easy: first a forward loop, then a backward loop:

```(define (lup-solve lu pi b)   (let* ((n (matrix-rows lu))          (y (make-vector n))          (x (make-vector n)))     (for (i n)       (vector-set! y i                    (- (vector-ref b (vector-ref pi i))                       (sum-of (* (lower lu i j) (vector-ref y j))                         (j range 0 i)))))     (for (i (- n 1) -1)       (vector-set! x i                    (/ (- (vector-ref y i)                          (sum-of (* (upper lu i j) (vector-ref x j))                            (j range (+ i 1) n)))                       (upper lu i i))))     x))```

Here’s an example:

```> (define a     #( #(1 2 0)        #(3 5 4)        #(5 6 3))) > (define b #(1/10 25/2 103/10)) > (call-with-values     (lambda () (lup-decomposition a))     (lambda (lu p) (lup-solve lu p b))) #(1/2 -1/5 3)```

We used sum-of and the matrix data type from the Standard Prelude. You can run the code at http://programmingpraxis.codepad.org/Eec8zXDy.

Pages: 1 2

### 4 Responses to “Solving Systems Of Linear Equations”

1. Remco Niemeijer said

Since I’m currently at a conference I don’t have as much time to devote to this as I’d like, but here’s my implementation of LU decomposition:

```lu :: Fractional a => [[a]] -> ([[a]], [[a]])
lu = unzip . map unzip . elim where
elim [] = []
elim ~((r:rs):xs) = zip (1 : repeat 0) (r:rs) :
zipWith (:) (map (\(y:_) -> (y / r, 0)) xs)
(elim \$ map sub xs) where
sub ~(y:ys) = zipWith (-) ys \$ map (y / r *) rs
```

If I have some time before Friday I’ll see if I can work on the rest.

Also, you might want to fix the examples in the exercise, since they contain some errors (e.g. in the first example: 1×3 + 3×4 + 1×0 + 0x0 = 15, not 19). The versions used in the solution code are correct.

2. programmingpraxis said

Fixed four errors. I didn’t do that very well, did I?

3. […] Praxis – Solving Systems Of Linear Equations By Remco Niemeijer In yesterday’s Programming Praxis exercise our task is to implement some more matrix-related functions, […]

4. Remco Niemeijer said

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

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

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

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

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)

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

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

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)
```