## Saturday, November 2, 2013

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