Matrix multiplication (functional approach in Python)

Matrix multiplication

The last entry looked at Gaussian elimination modeling matrices and vectors as tuples and using some idioms of functional programming. The linear algebra theme is continued here with programs to compute matrix/vector and matrix/matrix multiplications.

First some utilities to warm up on. Functions to project the rows or columns of a matrix.

def proj_r (A, i): return A[i]
def proj_c (A, j): return tuple (map (lambda t: (t[j],), A))

The function that scales a row vector by a constant is here again along with its symmetric partner for scaling a column vector.

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

def scale_c (c, k):
  def f (e):
    (xi,) = e
    return k*xi,
  return tuple (map (f, c))

For example, if c=((1,),(2,)), then scale_c (c, 2) is ((2,), (4,)).

The utility flatten reinterprets a sequence of column vectors as a matrix of row vectors.

def flatten (t):
  def row(i):
    def f (acc, y): 
        return acc + y[i]
    return functools.reduce (f, t, ())
  return tuple (map (row, range (0, len (t[0]))))

Here is an example. If t1=(1,),(2,) and t2=(3,),(4,) then s=(t1, t2) is a column vector sequence and flatten (s) is the matrix ((1, 3), (2, 4)).

Matrix vector multiplication comes next.

def mat_vec_mul (A, x):
  def f (i):
      (xi,) = x[i]
      return scale_c (proj_c (A, i), xi)
  M = tuple (map (f, range (0, len (x))))
  return tuple (map (lambda r : (sum (r), ), flatten (M)))

This computation utilizes the interpretation of the product as the sum of the column vectors of A weighted by the components of x.

Matrix multiplication follows immediately by seeing it as the flattening of a sequence of matrix/vector multiplications (that is of A on the individual columns of B).

def mat_mat_mul (A, B):
  f = lambda j : mat_vec_mul (A, proj_c (B, j))
  return flatten (tuple (map (f, range (0, len (B[0])))))

The transcript below shows the function at work in the Python top-level.

>>> A = ((2, 3), (4, 0) )
>>> B = ((1,  2, 0), (5, -1, 0) )
>>> print (str (mat_mat_mul (A, B)))
((17, 1, 0), (4, 8, 0))

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.

Term matching

Term matching

Extending one substitution with another can only work if the two substitutions are consistent with each other.

exception Term_match_exc

let extend_subst s1 s2 =
  let f acc (x, t)= 
    try 
      let u = List.assoc x acc in
      if t = u then acc
      else raise Term_match_exc
    with Not_found -> (x, t)::acc
in List.fold_left f s1 s2

A pattern match of a term u by a term t is a substitution that transforms t into u. If t is a variable s then the substitution (s, u) is an obvious pattern match for example. More generally though, if t is of the form f (t₁, ..., t_n), then

• u must be of the form f (u₁, ..., u_n);
  
• for each i there must exist a match of u_i with t_i;
  

and all the matches must be mutually compatible for a solution to exist.

let term_match (t1, t2) =
  let rec term_match' subst = function
    | (Var v, t) -> extend_subst [v, t] subst
    | (t, Var v) -> raise Term_match_exc
    | (Term (f, l), Term (g, r)) ->
      if f = g then 
        List.fold_left term_match' subst (List.combine l r)
      else raise Term_match_exc
  in term_match' [] (t1, t2)

For example, a pattern match of h(y(z, k(x, s))) and h(y(l(m, n,o()), k(t(), u))) is {(x,t()),(z, l(m,n,o()), (s, u)}.