# Shae's Ramblings

Stuff I find slightly meaningful and close to achievable

## 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:

• Approximation error $\mathrm{A} \approx \mathrm{LU}$ as a function of $n$
• Time complexity as a function of $n$

### 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.

## 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