LU decomposition from scratch

Implementation of the LU factorization in Python.

Description

Linear systems of the form ${\bf Ax}={\bf b}$ are essential in Linear Algebra and they appear in multiple applied contexts. When ${\bf A}$ is an $n \times n$ matrix, we know the $n \times 1$ solution vector ${\bf x}$ is unique, but that doesn’t mean it’s easy to find. We can use Gaussian elimination or find the inverse matrix ${\bf A}^{-1}$, but those methods are not feasible when $n$ is large.

A more efficient strategy (commonly used by computer software) consists of decomposing ${\bf A}$ into the product of a lower triangular matrix ${\bf L}$ and an upper triangular matrix ${\bf U}$. This new system $({\bf LU}){\bf x}={\bf b}$ is much easier to solve, eventhough it requires to first solve ${\bf L}{\bf y}={\bf b}$ and then ${\bf U}{\bf x}={\bf y}$ (See Singh (2013) for details).

Following the theory explained in Singh (2013), I decided to implement the method from scratch. The complete code can be found here

References

Singh, K. (2013). Linear algebra: step by step. Oxford University Press.