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

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:

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.

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

[…] 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, […]

My Haskell solution (see http://bonsaicode.wordpress.com/2010/07/21/programming-praxis-solving-systems-of-linear-equations/ for a version with comments):