Gaussian elimination (functional approach in Python)

Gaussian elimination

The goal here is to implement simple Gaussian elimination in Python, in a functional style just using tuples.

We view (a, b, c) a row vector and interpret ((a,),(b,),(c,)) as a column vector.

These first two elementary operations (scaling a row by a scalar and subtracting one row from another) come easily.

def scale_row (u, lam): 
    return tuple (
      map (lambda ui: lam*ui, u))

def subtract_rows (u, v) :
  return tuple (
    map (lambda i : u[i] - v[i], range (0, len (u))))

These alone enable us to proceed directly to forward elimination of an augmented coefficient matrix.

def elimination_phase (t):
  def elimination_step (tp, k): 
    pr = tp[k]
    def g (e):
      (i, r) = e
      if i <= k : 
        return r
      else:
        lam = r[k]/pr[k]
        return subtract_rows (r, scale_row (pr, lam))
    return tuple (
             map (g, tuple (
                      zip (range (0, len (tp)), tp))))
  return functools.reduce (
              elimination_step, range (0, len (t)), t)

Vector dot products are a straight forward fold.

def dot (u, v):
  def f (acc, i):
    (a, (b,)) = (u[i], v[i])
    return acc + a*b
  return functools.reduce (f, range (0, len (u)), 0.0)

With this, the back substitution phase can be written like this.

def back_substitution_phase (t):
  n = len (t)
  def back_substitution_step (x, k):
    bk = (t[k])[n]
    akk = (t[k])[k]
    xk = bk/akk if k == n-1 \
      else (bk - dot ((t[k][k+1:])[0:n-k-1], x))/akk
    return ((xk,),) + x
  return functools.reduce (
    back_substitution_step, range (len (t)-1, -1, -1), ())

Can we test it? Yes we can!

#Solve Ax=b with 
#
#        6   -4     1
#  A =  -4    6    -1
#        1   -4     6
#
#and b = (-14, 36, 6)^T.

A = (
    ( 6.,  -4.,  1.,  -14.)
  , (-4.,   6., -4.,   36.)
  , ( 1.,  -4.,  6.,    6.))

print (str (back_substitution_phase (elimination_phase (A))))

#Should give (10, 22, 14)^T.