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.

About these ads

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

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

    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)
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 629 other followers

%d bloggers like this: