LU Decomposition: Gaussian reduction

Ax=b

(ns morelinear.gauss
  (:require [clojure.core.matrix :refer :all]
            [clojure.core.matrix.linear :as linear]))
(set-current-implementation :vectorz)

Ax=b

This section is a placeholder of sorts and fills in some gaps in functionality in core.matrix. For a detailed explanation of the LU decomposition, see elinear

Note: The core.matrix library is embarassingly inadequate here. We now need to enable one of the backends (which I am loading at the top with the :vectorz) b/c it turns out the default core.matrix backend doesn't implement parts of its own API. (the lu in this case)

Gaussian Elimination gives us the form PA=LU and the core.matrix library gives us a method to do this clojure.core.matrix.linear/lu. We now need to add functions to do forward and backward substitution to actually solve the system. The functions are there but I have no idea how to access them b/c they're wrapped in some macro and not visible on any namespace I could poke at..

(defn- forward-substitution-rec
  ""
  [lower-triangular output-vector partial-input-vector step-number]

  (let [substitution-sum (scalar (inner-product (subvector (get-row lower-triangular
                                                                    0)
                                                           0
                                                           step-number)
                                                (subvector partial-input-vector
                                                           0
                                                           step-number)))

        new-input-value (/ (- (mget output-vector
                                    step-number)
                              substitution-sum)
                           (mget lower-triangular
                                 0
                                 step-number))]
    (if (= step-number
           (dec (ecount partial-input-vector)))
      (mset partial-input-vector 
            step-number
            new-input-value)
      (recur (matrix (submatrix lower-triangular
                                1
                                (dec (row-count lower-triangular))
                                0
                                (column-count lower-triangular)))
             output-vector
             (mset partial-input-vector 
                   step-number
                   new-input-value)
             (inc step-number)))))

(defn forward-substitution
  "Lx=b by forward subsitution, where L is lower triangular"
  [lower-triangular output-vector]
  (forward-substitution-rec lower-triangular
                            output-vector
                            (zero-array (shape output-vector))
                            0))

(defn backward-substitution
  "Ux=b by backward subsitution, where U is upper triangular"
  [upper-triangular output-vector]
  (let[rank (column-count upper-triangular)
       U (submatrix upper-triangular
                    0
                    rank
                    0
                    rank)
       b-short (subvector output-vector 0 rank)
       flipped-to-L (reshape (reverse (to-vector U))
                             (shape U))
       flipped-b (reshape (reverse (to-vector b-short))
                          (shape b-short))
       flipped-input (forward-substitution flipped-to-L
                                           flipped-b)]
    (reshape (reverse (to-vector flipped-input))
             (shape b-short))))
;; flips the matrix around to make it lower triangular..
;; then reuses the forward-substitution code
;; finally flips the result. A bit ugly, but short and easier to debug

Putting it all together to solve the square Ax=b system

(defn solve-with-lu
  "Solve Ax=b directly using Gaussian elimination with backward/forward substitution"
  [A b]
  (let [{:keys [L
                U
                P]} (linear/lu A)
        Pb (mmul P b)
        y (forward-substitution L Pb)]
    (backward-substitution U y)))

Test example:

(let [skinny-A [[ 1.0 2.0 ]
                [ 0.0 3.0 ]
                [ 0.0 0.0 ]]
      long-b [ 5.0 6.0 999.0 ]]
  (pm (backward-substitution skinny-A
                             long-b)))
;;[1.000 2.000]

(let [A [[0.4 0.6 0.8 0.9]
         [0.4 0.6 0.6 0.8]
         [0.2 0.2 0.4 0.2]
         [0.5 0.3 0.3 0.8]]
      {L :L
       U :U
       P :P} (linear/lu A)
      b [2.0 5.0 8.0 1.0]]
  (pm (mmul L
            (forward-substitution L
                                  b)))
  ;;[2.000 5.000 8.000 1.000]
  (pm (mmul U
            (backward-substitution U
                                   b)))
  ;;[2.000 5.000 8.000 1.000]
  (pm (mmul A (solve-with-lu A (mmul P b)))))
  ;;[2.000 5.000 8.000 1.000]

LU Decomposition: Gaussian reduction

Table of Contents

Ax=b