Programming Praxis

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