Matrix Inverse code in Scheme for Arbitrary Fields

Scheme code for Inverting Matrices in Arbitrary Fields

There is now a repository version of this code. N.B. Full Scheme implementations provide complex numbers, and the alternatives of floating point or rational arithmetic. This Scheme code also allows you to provide your own field. See too OCaml matrix code here, and Algol 68 matrix inversion. See also some sparse matrix code in the context of linear mesh problems.

The following Scheme code takes tools for some field and returns a few procedures to process matrices with elements from that field. Use the definition of trnxx.

(define matpak (lambda (zero zer? one fa fs fm fi) (letrec (
   (tf (lambda (x zp) (if (null? x) (zp (list one))
     (if (zer? (caar x)) (let ((z (tf (cdr x) zp)))
     (cons (car z) (cons (car x) (cdr z)))) x))))
   (tr (trnxx zero zer?))
   (sm (lambda (s a)(map (lambda (w) (fm s w)) a)))
   (fms (lambda (a s b) (if (null? a) (sm (fs zero s) b) (if (null? b) a
     (cons (fs (car a) (fm s (car b))) (fms (cdr a) s (cdr b)))))))
   (deter (lambda (a) (letrec ((p #f)
     (tf (lambda (x) (if (null? x) (list one)
       (if (zer? (caar x)) (let ((z (tf (cdr x))))
       (set! p (not p))
       (cons (car z) (cons (car x) (cdr z)))) x)))))
     (let inx ((d one)(a a))
          (if (null? a) (if p (fs zero d) d) (let* (
            (A (tf a))
            (i (fi (caar A)))
            (b (map (lambda (z) (fm z i)) (cdar A))))
         (inx (fm (caar A) d) (map (lambda (x w) (fms x w b))
          (map cdr (cdr A)) (map car (cdr A))))))))))
   (inv (lambda (a nxp) (let ol ((ut (let inx (
      (a (let pad ((x a)(e (list one)))
         (if (null? x) '() (cons (let ap ((z (car x))(ln a))
           (if (null? ln) e (if (null? z) (cons zero (ap z (cdr ln)))
             (cons (car z)(ap (cdr z)(cdr ln))))))
          (pad (cdr x) (cons zero e))))))
      (np nxp))
          (if (null? a) '() (let* (
            (A (tf a np))
            (i (fi (caar A)))
            (b (map (lambda (z) (fm z i)) (cdar A))))
         (cons b (inx (map (lambda (x w) (fms x w b))
          (map cdr (cdr A)) (map car (cdr A)))
          (lambda (w) (np (cons (fs zero (ip w b)) w))))))))))
      (if (null? ut) '() (cons
        (let eg ((top (car ut))(bod (cdr ut))) (if (null? bod) top
          (eg (fms (cdr top) (car top) (car bod))(cdr bod))))
      (ol (cdr ut)))))))
   (ip (lambda (x y)(if (or (null? x) (null? y)) zero (fa (fm (car x)(car y))
        (ip (cdr x)(cdr y))))))
   (mp (lambda (a b)(let ((b (tr b)))
       (map (lambda(x) (map (lambda (y) (ip x y)) b)) a)))))

(list mp inv ip tr deter))))

The field tools provided to matpak are:

the zero element for the field
a predicate that tests its field argument for being zero
the multiplicative identity of the field
the addition operator that takes two field arguments and returns their sum
the subtraction operator
the field multiplication operator for two field arguments
the multiplicative inverse function

With these inputs matpak returns a list of functions:

Matrix multiplication
Matrix inversion
Inner product of two vectors
Transpose
Determinant

Matrices are composed of lists of lists of field values. Here are comments on the code. The second argument to the inversion routine is an escape procedure of one argument which is called if the matrix is singular. When a singular matrix is passed to the inversion routine, that routine does not return but the escape procedure is called instead with a list of k numbers which are coefficients of a linear combination of the first k columns which is exactly zero. The escape procedure should not itself return. In the examples below we use cc below which depends on the definitions of exit and mer.

(define cc (mer "Singular:"))

Examples

The following supplies Scheme’s standard arithmetic operators which deal with exact rational values, and uses the resulting procedures to invert a real matrix:

(apply (lambda (matm matinv ip tr det)

; A test matrix inversion with rational arithmetic.
(let* ((a '((0 2 4 5) (3 4 5 -2) (7 6 5 3) (4 6 5 7)))
   (b (matinv a cc)))
   (list b (matm a b)(matinv b cc))))

; Here we supply matpak with the field goodies for the rationals.
(matpak 0 zero? 1 + - * /))

For this, both Gambit and MzScheme for the Mac yield:

(((5/36 -7/36 5/12 -1/3)
(-155/288 73/288 -47/96 2/3)
(7/18 1/18 1/6 -1/3)
(5/48 -7/48 1/16 0))

((1 0 0 0) (0 1 0 0) (0 0 1 0) (0 0 0 1))

((0 2 4 5) (3 4 5 -2) (7 6 5 3) (4 6 5 7)))

in perfect perfect agreement with theory.

To invoke these matrix functions interactively one might say

(define ops (matpak 0 zero? 1 + - * /))
(define mm (car ops))
(define inv (cadr ops))
(define det (car(cddddr ops)))

(mm '((0 1)(1 0)) '((2 0)(3 0))) ; => ((3 0) (2 0))
(inv '((2 0)(0 5)) cc) ; => ((1/2 0) (0 1/5))

One must supply enough 0’s on the end of each row to make the right length.

To illustrate the report from a singular matrix:

(inv '((2 3 4) (6 9 3) (4 6 6)) cc)

yields:

("Exception:" "Singular:" (-3/2 1))

Which notes that −3/2 times the first column plus the second column is 0.

Now we adjoin the square root of 2 to the field of rationals where for rationals x and y, (cons x y) denotes x+y√2. This is a new field!

(define ops (matpak
    '(0 . 0)
    (lambda (x) (and (zero? (car x)) (zero? (cdr x))))
    '(1 . 0)
    (lambda (x y) (cons (+ (car x)(car y)) (+ (cdr x)(cdr y))))
    (lambda (x y) (cons (- (car x)(car y)) (- (cdr x)(cdr y))))
    (lambda (x y) (cons (+ (* (car x)(car y)) (* 2 (cdr x)(cdr y)))
               (+ (* (car x)(cdr y)) (* (cdr x)(car y)))))
    (lambda (x)(let* ((a (car x))(b (cdr x))(d (- (* a a) (* 2 b b))))
      (cons (/ a d) (- (/ b d)))))))

(define mm (car ops))
(define inv (cadr ops))
(define det (car (cddddr ops)))

Now we can exactly invert a matrix with irrational elements like (2 . 1) = 2+√2:

(define m '(((2 . 0) (5 . 0) (2 . 0))
   ((2 . 1) (7 . 0) (3 . 0))
   ((5 . 0) (5 . 0) (2 . 0))))
(inv m cc)
(det m)
(det (inv m cc))
(mm m (inv m cc))

Here is a source for the infamous Hilbert matrix. matpak works with finite fields too.