lisp308

A collection of matrix utility functions.

I wrote most of this while taking Math 308, Linear Algebra, at the University of Washington. Programming isn’t required by the course, but it’s good fun!

Features & Examples

Echelon Form Reduction

Can reduce matrices to echelon or reduced echelon form, and records the operations taken to get there.

(reduce-ef '((1 2 3)
             (4 5 6)
             (7 8 9)))

1	2	3
0	-3	-6
0	0	0

(reduce-ref '((1 2 3)
              (4 5 6)
              (7 8 9)))

1	0	-1
0	1	2
0	0	0

The raw return values, including row operations (in reverse order):

((1 0 -1) (0 1 2) (0 0 0))
((:SWAP 0 2) (:* 1 2) (:* -1/3 1) (:+ 2/3 1 2) (:SWAP 0 2) (:+ -2 1 2)
 (:+ -7 0 2) (:+ -4 0 1))

So, for example, the first row operation performed was to subtract 4 of row 0 from row 1, as represented by (:+ -4 0 1). These operations can be used to easily calculate determinants, inverses, etc.

Matrix algebra

(setq a '((1 2 3)
          (4 5 6)
          (7 8 9)))
(setq b '((2 0 0)
          (0 2 0)
          (0 0 2)))
(mat* b a)

2	4	6
8	10	12
14	16	18

Equivalently:

(mat* '((1 2 3)
        (4 5 6)
        (7 8 9))
      (mat-scalar* (make-identity 3) 2))

We may also invert:

(invert '((1 2 3)
          (4 5 6)
          (9 2 7)))

-23/36	2/9	1/12
-13/18	5/9	-1/6
37/36	-4/9	1/12

The default invert implementation uses row operations, but there is another one included that uses cofactors and determinant.

Other matrix algebra utilities include:

(mat+): Addition
(mat-expt): Exponentiation
(transpose): Exchange rows with columns

Spaces

(column-space-basis '((1 2 3)
                      (4 5 6)
                      (7 8 9)))

((1 4 7) (2 5 8))

Here it properly detects that the third column can be represented by the first two.

(nullspace-basis '((1 2 3)
                   (4 5 6)
                   (7 8 9)))

((1 -2 1))

This is the correct one-dimensional nullspace.

Determinants

(determinant (make-identity 3))

(determinant (mat-scalar* (make-identity 3) 5))

(determinant '((1 2 3)
               (4 5 6)
               (7 8 9)))

The default determinant implementation uses row operations. Alternate implementations are provided that calculate the determinant using the permutation expansion (generate all permutation matrices, multiply corresponding values from original matrix, then sum) and the Laplace cofactor expansion.

Eigen{values,vectors}

(eigenvalues '((1 2 3)
               (4 5 6)
               (7 8 9)))

(eigenvalues '((0 1)
               (1 1)))

1.618034

-0.618034

(characteristic-polynomial '((0 1)
                             (1 1)))

-1

The polynomial is outputted from least to greatest power, so this is read as x^2-x-1.

Eigenvalues are not calculated in a very reliable way. For 1x1 and 2x2 matrices, I simply use the quadratic formula, which does always find the results. But for larger matrices, a trial-and-error search is conducted on the inegers from -1000 to 1000; larger or fractional eigenvalues will not be found. We only search for real eigenvalues in general.

(eigenspace-basis '((0 1)
                    (1 1))
                  1.618034)

((0.618034 1))

In this case the eigenspace was one-dimensional.

(diagonalize '((0 1)
               (1 1)))

((1.618034 0)
 (0       -0.618034))

((0.618034 -1.618034)
 (1         1))

The two matrices are D and P respectively, where PDP^(-1) = the original matrix.

Other functions

Most of these are tiny utilities that re-use the functionality from the main functions above.

(rank)
(nullity)
()

markasoftware/lisp308