numpy scipy gaussian elimination using LU decomposition with pivoting

Gaussian_elimination.py

#! /usr/bin/env python

"""
Solve linear system using LU decomposition and Gaussian elimination
"""

import numpy as np
from scipy.linalg import lu, inv

def gausselim(A,B):
    """
    Solve Ax = B using Gaussian elimination and LU decomposition.

    A = LU   decompose A into lower and upper triangular matrices
    LUx = B  substitute into original equation for A

    Let y = Ux and solve:
    Ly = B --> y = (L^-1)B  solve for y using "forward" substitution
    Ux = y --> x = (U^-1)y  solve for x using "backward" substitution

    :param A: coefficients in Ax = B
    :type A: numpy.ndarray of size (m, n)
    :param B: dependent variable in Ax = B
    :type B: numpy.ndarray of size (m, 1)
    """
    # LU decomposition with pivot
    pl, u = lu(A, permute_l=True)
    # forward substitution to solve for Ly = B
    y = np.zeros(B.size)
    for m, b in enumerate(B.flatten()):
        y[m] = b
        # skip for loop if m == 0
        if m:
            for n in xrange(m):
                y[m] -= y[n] * pl[m,n]
        y[m] /= pl[m, m]

    # backward substitution to solve for y = Ux
    x = np.zeros(B.size)
    lastidx = B.size - 1  # last index
    for midx in xrange(B.size):
        m = B.size - 1 - midx  # backwards index
        x[m] = y[m]
        if midx:
            for nidx in xrange(midx):
                n = B.size - 1  - nidx
                x[m] -= x[n] * u[m,n]
        x[m] /= u[m, m]
    return x

if __name__ == '__main__':
    x = gausselim(np.array([[3, 2], [1, -4]]), np.array([[5], [10]]))
    print x

Isn't it the same as solve_lu by scipy?

Hi @Wikunia,

Yes they're probably functionally the same, but my goal here was to understand Gaussian elimination using LU decomposition simply using pure Python. If you read my blog post, you'll see this was just for fun, to understand it for my own education. I was not and would not ever recommend anyone to use this Gist over the existing SciPy implementation. But if anyone were interested in learning how it works, then they might find this Gist to be useful. I personally thought that this Gaussian elimination algorithm was super cool, and very clever, and I really liked how LU decomposition just means Lower/Upper. Math is so fun! And I love learning!

The scipy.linalg.solve_lu function calls LAPACK and uses the double getrs FORTRAN function. This then calls several BLAS functions. See netlib.org to explore the HTML docs. This of course is awesome for solving systems of linear equations, because these methods are tested, robust, efficient, well-documented, and established.

Thanks for this clarification! 😁

Hi @mikofski,
actually I can't remember why I commented here. I appreciate people who code things from scratch to improve their understanding and I do the same on my blog. Enjoy the rest of your weekend!

Hello @mikofski, I am a new Python learner. I am trying to do Gaussian elimination using LU decomposition using Python as well but I am trying to do it with test matrices are stored in the adjacency list (in each row of the file we have three numbers) something like this:
23 3 0.000001370542294
4 4 0.107816040610854
7 4 0.022782277293175
11 4 -0.00921782470662
And file might have 25 or 50 rows.
Can you give me advice on how I could read .txt file and implement the code?
Thanks

Hi @sbykl , I think what you're looking for is either numpy.loadtxt or pandas.read_csv. Read the docs, hopefully you'll find usage straightforward. 😃