# LU Decomposition: Gaussian reduction

## Table of Contents

(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]