LU Decomposition

LU Decomposition or LU Factorization factors a non singular matrix A as the product of a lower triangular matrix L, and an upper triangular matrix U such that A = LU.
The product may also involve a permutation matrix P, in which case it is PA = LU
LU factorization is mainly used for solving systems of linear equations. It is also used for computing the inverse or the determinant of a matrix. The key reason for using LU decomposition is that linear systems having triangular coefficients matrices are very easy to solve.

LU decomposition essentially records the steps used in Gaussian elimination on a matrix A.

For instance, consider the process of Gaussian Elimination for the matrix

A = [[ 1  -2  3]
     [ 2  -5  12]
     [ 0   2 -10]]

First, we substract 2 times the first row from the second row. We indicate this by a (2).

[[ 1  -2  3]
 [(2) -1  6]
 [ 0   2 -10]]

Then, in the third row, we don't need to do any thing to eliminate the first element. This is denoted by a (0)

[[ 1  -2  3]
 [(2) -5  12]
 [(0)  2 -10]]

Followed by subtraction of -2 times the second row from the third

[[ 1   -2  3]
 [(2)  -5  12]
 [(0) (-2) 2]]

Then, U is the resultant matrix, and L is the matrix representing the operations performed, with ones along the diagonals.

L = [[ 1  0  0]            U = [[ 1 -2  3]]
     [ 2  1  0]                 [ 0 -1  6]
     [ 0 -2  1]]                [ 0  0  2]

and

     [[ 1  0  0]    [[ 1 -2  3]]   [[ 1  -2  3]
LU =  [ 2  1  0]  x  [ 0 -1  6]  =  [ 2  -5  12]  = A
      [ 0 -2  1]]    [ 0  0  2]     [ 0   2 -10]]

But, every matrix does not necessarily have an LU decomposition. Some times we need to perform row exchange operations on the matrix A for it to be decomposable in the product LU. In such cases, A is pre multiplied by P, which is a pivot matrix.
This gives us the equation PA = LU

This is how to perform LU decomposition in core.matrix

(let [{:keys [L U P]} (lu A)] ....)
; => returns all three matrices
(let [{:keys [L]} (lu A (:return [:L]))] ....)
; => returns just L

Solving Ax=b

To solve a system of linear equations Ax=b,
we first multiply both sides by P, the pivot matrix
                         PAx = Pb
this is same as
                         LUx = d        (d = Pb)
This leaves us with two triangular system problems, first substitute Ux with y and solve
                         Ly = d
and then solve
                         Ux = y

Computing inverse of a matrix

If we want to compute the inverse of a square matrix A of size n, We simply compute the solution for Ax_i=b_i for i=1. . .n, where b_i is the i^th coulmn of I_n, an identity matrix of size n.
We get x₁, . . .,x_n which are the columns of A^-1, the inverse of A.

core.matrix specially supports this to make inverting easier. If you only pass the parameter A when calling the solve function, it'll automatically solve for every column of I_n and return the inverse of A.

(let [X (solve A)] ....)
; => X is inverse of A

Determinant of a matrix

If we have the LU decomposition of a matrix A, its determinant can be very easily computed using the knowledge that the determinant of a triangular matrix is the product of its diagonal elements.

det(A) = det(P^-1LU)
           = det(P^-1)det(L)det(U)
           = (-1)^s(∏ⁿ_i=1u_ii)            where s is the number of row exchanges.

(let [{:keys [L U P]} (li/lu c)] 
  (* (det P) 
     (reduce * 1 (diagonal U)) 
     (reduce * 1 (diagonal L))))