Every square matrix 



A


{\displaystyle A}

 can be decomposed into a product of a lower triangular matrix 



L


{\displaystyle L}

 and a upper triangular matrix 



U


{\displaystyle U}

,
as described in LU decomposition.





A
=
L
U


{\displaystyle A=LU}


It is a modified form of Gaussian elimination.
While the Cholesky decomposition only works for symmetric,
positive definite matrices, the more general LU decomposition
works for any square matrix.

There are several algorithms for calculating L and U.
To derive Crout's algorithm for a 3x3 example,
we have to solve the following system:





A
=


(




a

11





a

12





a

13







a

21





a

22





a

23







a

31





a

32





a

33





)


=


(




l

11




0


0





l

21





l

22




0





l

31





l

32





l

33





)




(




u

11





u

12





u

13






0



u

22





u

23






0


0



u

33





)


=
L
U


{\displaystyle A={\begin{pmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\\\end{pmatrix}}={\begin{pmatrix}l_{11}&0&0\\l_{21}&l_{22}&0\\l_{31}&l_{32}&l_{33}\\\end{pmatrix}}{\begin{pmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\end{pmatrix}}=LU}


We now would have to solve 9 equations with 12 unknowns. To make the system uniquely solvable, usually the diagonal elements of 



L


{\displaystyle L}

 are set to 1






l

11


=
1


{\displaystyle l_{11}=1}







l

22


=
1


{\displaystyle l_{22}=1}







l

33


=
1


{\displaystyle l_{33}=1}


so we get a solvable system of 9 unknowns and 9 equations.





A
=


(




a

11





a

12





a

13







a

21





a

22





a

23







a

31





a

32





a

33





)


=


(



1


0


0





l

21




1


0





l

31





l

32




1



)




(




u

11





u

12





u

13






0



u

22





u

23






0


0



u

33





)


=


(




u

11





u

12





u

13







u

11



l

21





u

12



l

21


+

u

22





u

13



l

21


+

u

23







u

11



l

31





u

12



l

31


+

u

22



l

32





u

13



l

31


+

u

23



l

32


+

u

33





)


=
L
U


{\displaystyle A={\begin{pmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\\\end{pmatrix}}={\begin{pmatrix}1&0&0\\l_{21}&1&0\\l_{31}&l_{32}&1\\\end{pmatrix}}{\begin{pmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\end{pmatrix}}={\begin{pmatrix}u_{11}&u_{12}&u_{13}\\u_{11}l_{21}&u_{12}l_{21}+u_{22}&u_{13}l_{21}+u_{23}\\u_{11}l_{31}&u_{12}l_{31}+u_{22}l_{32}&u_{13}l_{31}+u_{23}l_{32}+u_{33}\end{pmatrix}}=LU}


Solving for the other 



l


{\displaystyle l}

 and 



u


{\displaystyle u}

, we get the following equations:






u

11


=

a

11




{\displaystyle u_{11}=a_{11}}







u

12


=

a

12




{\displaystyle u_{12}=a_{12}}







u

13


=

a

13




{\displaystyle u_{13}=a_{13}}







u

22


=

a

22


−

u

12



l

21




{\displaystyle u_{22}=a_{22}-u_{12}l_{21}}







u

23


=

a

23


−

u

13



l

21




{\displaystyle u_{23}=a_{23}-u_{13}l_{21}}







u

33


=

a

33


−
(

u

13



l

31


+

u

23



l

32


)


{\displaystyle u_{33}=a_{33}-(u_{13}l_{31}+u_{23}l_{32})}


and for 



l


{\displaystyle l}

:






l

21


=


1

u

11





a

21




{\displaystyle l_{21}={\frac {1}{u_{11}}}a_{21}}







l

31


=


1

u

11





a

31




{\displaystyle l_{31}={\frac {1}{u_{11}}}a_{31}}







l

32


=


1

u

22




(

a

32


−

u

12



l

31


)


{\displaystyle l_{32}={\frac {1}{u_{22}}}(a_{32}-u_{12}l_{31})}


We see that there is a calculation pattern, which can be expressed as the following formulas, first for 



U


{\displaystyle U}








u

i
j


=

a

i
j


−

∑

k
=
1


i
−
1



u

k
j



l

i
k




{\displaystyle u_{ij}=a_{ij}-\sum _{k=1}^{i-1}u_{kj}l_{ik}}


and then for 



L


{\displaystyle L}








l

i
j


=


1

u

j
j




(

a

i
j


−

∑

k
=
1


j
−
1



u

k
j



l

i
k


)


{\displaystyle l_{ij}={\frac {1}{u_{jj}}}(a_{ij}-\sum _{k=1}^{j-1}u_{kj}l_{ik})}


We see in the second formula that to get the 




l

i
j




{\displaystyle l_{ij}}

 below the diagonal, we have to divide by the diagonal element (pivot) 




u

j
j




{\displaystyle u_{jj}}

, so we get problems when 




u

j
j




{\displaystyle u_{jj}}

 is either 0 or very small, which leads to numerical instability.

The solution to this problem is pivoting 



A


{\displaystyle A}

, which means rearranging the rows of 



A


{\displaystyle A}

, prior to the 



L
U


{\displaystyle LU}

 decomposition, in a way that the largest element of each column gets onto the diagonal of 



A


{\displaystyle A}

. Rearranging the rows means to multiply 



A


{\displaystyle A}

 by a permutation matrix 



P


{\displaystyle P}

:





P
A
⇒

A
′



{\displaystyle PA\Rightarrow A'}


Example:







(



0


1




1


0



)




(



1


4




2


3



)


⇒


(



2


3




1


4



)




{\displaystyle {\begin{pmatrix}0&1\\1&0\end{pmatrix}}{\begin{pmatrix}1&4\\2&3\end{pmatrix}}\Rightarrow {\begin{pmatrix}2&3\\1&4\end{pmatrix}}}


The decomposition algorithm is then applied on the rearranged matrix so that





P
A
=
L
U


{\displaystyle PA=LU}




Task description

The task is to implement a routine which will take a square nxn matrix 



A


{\displaystyle A}

 and return a lower triangular matrix 



L


{\displaystyle L}

, a upper triangular matrix 



U


{\displaystyle U}

 and a permutation matrix 



P


{\displaystyle P}

,
so that the above equation is fulfilled.

You should then test it on the following two examples and include your output.



Example 1

A

1   3   5
2   4   7
1   1   0

L

1.00000   0.00000   0.00000
0.50000   1.00000   0.00000
0.50000  -1.00000   1.00000

U

2.00000   4.00000   7.00000
0.00000   1.00000   1.50000
0.00000   0.00000  -2.00000

P

0   1   0
1   0   0
0   0   1

Example 2

A

11    9   24    2
 1    5    2    6
 3   17   18    1
 2    5    7    1

L

1.00000   0.00000   0.00000   0.00000
0.27273   1.00000   0.00000   0.00000
0.09091   0.28750   1.00000   0.00000
0.18182   0.23125   0.00360   1.00000

U

11.00000    9.00000   24.00000    2.00000
 0.00000   14.54545   11.45455    0.45455
 0.00000    0.00000   -3.47500    5.68750
 0.00000    0.00000    0.00000    0.51079

P

1   0   0   0
0   0   1   0
0   1   0   0
0   0   0   1



 
module Main 
 
import Data.Vect
 
Matrix : Nat -> Nat -> Type -> Type
Matrix m n t = Vect m (Vect n t)
 
-- Creates list from 0 to n (not including n) 
upTo : (m : Nat) -> Vect m (Fin m)
upTo Z = []
upTo (S n) = 0 :: (map FS (upTo n))
 
-- Creates list from 0 to n-1  (not including n-1)
upToM1 : (m : Nat) -> (sz ** Vect sz (Fin m))
upToM1 m = case (upTo m) of
                  (y::ys) => (_ ** init(y::ys)) 
                  [] => (_ ** [])
 
-- Creates list from i to n (not including n)
fromUpTo : {n : Nat} -> Fin n -> (sz ** Vect sz (Fin n))
fromUpTo {n} m = filter (>= m) (upTo n)
 
-- Creates list from i+1 to n (not including n)
fromUpTo1 : {n : Nat} -> Fin n -> (sz ** Vect sz (Fin n))
fromUpTo1 {n} m with (fromUpTo m)
  | (_ ** xs) = case xs of 
                  (y::ys) => (_ ** ys)
                  [] => (_ ** [])
 
 
-- Create Zero Matrix of size m by n 
zeros : (m : Nat) -> (n : Nat) -> Matrix m n Double
zeros m n = replicate m (replicate n 0.0)
 
replaceAtM : (Fin m, Fin n) -> t -> Matrix m n t -> Matrix m n t
replaceAtM (i, j) e a = replaceAt i (replaceAt j e (index i a)) a
 
-- Create Identity Matrix of size m by m
eye : (m : Nat) -> Matrix m m Double
eye m = map create1Vec (upTo m)
  where
    set1 : Vect m Double -> Fin m -> Vect m Double
    set1 a n = replaceAt n 1.0 a
 
    create1Vec : Fin m -> Vect m Double
    create1Vec n = set1 (replicate m 0.0) n
 
 
indexM : (Fin m, Fin n) -> Matrix m n t -> t
indexM (i, j) a = index j (index i a)
 
 
-- Obtain index for the row containing the
-- largest absolute value for the given column
colAbsMaxIndex : Fin m -> Fin m -> Matrix m m Double -> Fin m
colAbsMaxIndex startRow col a {m} with (fromUpTo startRow)
  | (_ ** xs) = 
    snd $ foldl (\(absMax, idx), curIdx => 
          let curAbsVal = abs(indexM (curIdx, col) a) in
            if (curAbsVal > absMax) 
              then (curAbsVal, curIdx)
              else (absMax, idx)
        ) (0.0, startRow) xs
 
 
-- Swaps two rows in a given matrix
swapRows : Fin m -> Fin m -> Matrix m n t -> Matrix m n t
swapRows r1 r2 a = replaceAt r2 tempRow $ replaceAt r1 (index r2 a) a
  where tempRow = index r1 a
 
 
-- Swaps two individual values in a matrix
swapValues : (Fin m, Fin m) -> (Fin m, Fin m) -> Matrix m m Double -> Matrix m m Double
swapValues (i1, j1) (i2, j2) m = replaceAtM (i2, j2) v1 $ replaceAtM (i1, j1) v2 m
  where
      v1 = indexM (i1, j1) m
      v2 = indexM (i2, j2) m
 
-- Perform row Swap on Lower Triangular Matrix
lSwapRow : Fin m -> Fin m -> Matrix m m Double -> Matrix m m Double
lSwapRow row1 row2 l {m} with (filter (< row1) (upTo m))
  | (_ ** xs) =  foldl (\l',col => swapValues (row1, col) (row2, col) l') l xs
 
 
rowSwap : Fin m -> (Matrix m m Double,  Matrix m m Double, Matrix m m Double) -> 
                        (Matrix m m Double, Matrix m m Double, Matrix m m Double)
rowSwap col (l,u,p) = (lSwapRow col row l, swapRows col row u, swapRows col row p) 
      where row = colAbsMaxIndex col col u
 
 
calc : (Fin m) -> (Fin m) -> (Matrix m m Double, Matrix m m Double) -> 
                                (Matrix m m Double, Matrix m m Double)
calc i j (l, u) {m} = (l', u')
   where 
         l' : Matrix m m Double
         l' = replaceAtM (j, i) ((indexM (j, i) u) / indexM (i, i) u) l
 
         u'' : (Fin m) -> (Matrix m m Double) -> (Matrix m m Double)
         u'' k u = replaceAtM (j, k) ((indexM (j, k) u) - 
                  ((indexM (j, i) l') * (indexM (i, k) u))) u
 
         u' : (Matrix m m Double)
         u' with (fromUpTo i) | (_ ** xs) = foldl (\curU, idx => u'' idx curU) u xs
 
 
-- Perform a single iteration of the algorithm for the given column
iteration : Fin m -> (Matrix m m Double, Matrix m m Double, Matrix m m Double) ->
                        (Matrix m m Double, Matrix m m Double, Matrix m m Double) 
iteration i lup {m} = iterate' (rowSwap i lup) 
 
          where 
                modify : (Matrix m m Double, Matrix m m Double) -> 
                            (Matrix m m Double, Matrix m m Double)
                modify lu with (fromUpTo1 i) | (_ ** xs) = 
                                            foldl (\lu',j => calc i j lu') lu xs 
 
                iterate' : (Matrix m m Double, Matrix m m Double, Matrix m m Double) ->
                              (Matrix m m Double, Matrix m m Double, Matrix m m Double)
                iterate' (l, u, p) with (modify (l, u)) | (l', u') = (l', u', p) 
 
 
-- Generate L, U, P matricies from a given square matrix.
-- Where L * U = A, and P is the permutation matrix
luDecompose : Matrix m m Double -> (Matrix m m Double, Matrix m m Double, Matrix m m Double) 
luDecompose a {m} with (upToM1 m) 
  | (_ ** xs) = foldl (\lup,idx => iteration idx lup) (eye m,a,eye m) xs 
 
 
 
ex1 : (Matrix 3 3 Double, Matrix 3 3 Double, Matrix 3 3 Double)
ex1 = luDecompose [[1, 3, 5], [2, 4, 7], [1, 1, 0]] 
 
ex2 : (Matrix 4 4 Double, Matrix 4 4 Double, Matrix 4 4 Double)
ex2 = luDecompose [[11, 9, 24, 2], [1, 5, 2, 6], [3, 17, 18, 1], [2, 5, 7, 1]]
 
printEx : (Matrix n n Double, Matrix n n Double, Matrix n n Double) -> IO ()
printEx (l, u, p) = do
  putStr "l:"
  print l
  putStrLn "\n"
 
  putStr "u:"
  print u
  putStrLn "\n"
 
  putStr "p:"
  print p
  putStrLn "\n"
 
main : IO()
main = do 
  putStrLn "Solution 1:"
  printEx ex1
  putStrLn "Solution 2:"
  printEx ex2
 