LU decomposition in three languages

This post is about implementing an algorithm which finds the lower-upper decomposition or LU decomposion of a matrix \(\mathrm{A} \in \mathbb{R}^{n\times n}\). The letter \(n\) will always refer to the shape of the square matrix \(\mathrm{A}\).

Pseudo code

We use the pseudo-code from the following lecture IIT notes.

Using python syntax, the code for the factorization can be written

U, L = A, I
for k in range(0, n-1):
    for j in range(k+1, n):
        L[j,k] = U[j,k] / U[k,k]
        U[j, k:n] = U[j, k:n] - L[j,k]U[k, k:n]

where I stands for the identity \(\mathrm{I}_n\).

First example

Let us pick and factorize a matrix for \(n=4\). Our example matrix is

$$ \mathrm{A} = \begin{bmatrix} -1.076 & 0.657 & -1.222 & -0.467 \\ -1.101 & 0.805 & -1.247 & -0.62 \\ -0.594 & -0.112 & -0.327 & 0.13 \\ 0.412 & 0.604 & 0.398 & 0.472 \end{bmatrix} $$

Our python code provides the factorization (rounded to the third decimal place):

$$ \mathrm{L} = \begin{bmatrix} 1. & 0. & 0. & 0. \\ 1.023 & 1. & 0. & 0. \\ 0.552 & -3.576 & 1. & 0. \\ -0.383 & 6.446 & -0.255 & 1. \end{bmatrix} \;;\; \mathrm{U} = \begin{bmatrix} -1.076 & 0.657 & -1.222 & -0.467 \\ 0. & 0.133 & 0.003 & -0.142 \\ 0. & 0. & 0.36 & -0.121 \\ 0. & 0. & 0. & 1.179 \end{bmatrix} $$

The error in approximating \(\mathrm{A} \approx \mathrm{LU}\) (due to floating point operations) can be computed as $$ E(\mathrm{A}) = \sum_{i,j = 1}^{n}\left\lvert \mathrm{A}_{ij} – \left[\mathrm{LU}\right]_{ij}\right\rvert \approx 6.939 \cdot 10^{-17}$$

which is extremely small.

Python code and performance

Let's analyze our python program. We want to estimate two things:

Numpy routine

To estimate this, we sample for each \(n\) a number \(m=100\) of randomly generated matrices (using numpy's rand.randn routine) and run our factorization routine which operates on numpy arrays.

To be added.

JAX routine

To be added.

Julia routine

In julia, the identity matrix I is built so that algebraic properties of the identity are used to simplify computations without storing the \(n\times n\) matrix in memory. Our algorithm however needs an identity matrix which can be modified. With the initialization lower = I(n), the following error occurs:

ERROR: LoadError: ArgumentError: cannot set off-diagonal entry (2, 1) to a nonzero value (2.331270057426189)

this is because julia's I is a constant. In order to have a modifiable copy of I, we use Matrix{Float64}(I, n, n) instead.

Our function to compute the LU decomposition is thus:

using LinearAlgebra

function lu_factorize(matrix)

    n = size(matrix)[1]
    upper, lower = deepcopy(matrix), Matrix{Float64}(I, n, n)
    for k = 1:(n-1)
        for j = (k+1):n
            lower[j,k] = upper[j,k] / upper[k,k]
            upper[j, k:n] = upper[j, k:n] - lower[j,k] * upper[k, k:n]
        end
    end
    return lower, upper
end

R routine

To be added.