Linear Systems in Lisp

Intro
Matices
Matrix Operations
Elementary Matrices
TODOs

Intro

I'm going through Matrix Analysis & Applied Linear Algebra (Carl Meyer) and implementing the algorithms in ELisp. An excuse to learn ELisp and solidy my linear algebra. There is no emphasis on performance - just on clarity and extensability

Matices

Linear system (ie. a system of equations) are for convenience represented using matrices. Once we define a matrix and all operations on it we can begin to use them for manipulating linear systems. The standard way of representing a matrix that you often see in C is using multidimensional arrays. For a 2D matrix this has 3 values: number of rows number of columns data

data is a long list of size num-row * num-col. When making a new matrix it's useful to be able to make a list of zeroes of the required size

(defun build-list-of-zeroes (size)
  "Recursively builds a list of zeroes of the given size"
  (if (zerop size) 
      '() ;base case
    (cons 0 (build-list-of-zeroes (1- size))))) ; "else"/iterative step

zerop tests if the value is zero

() with a quote is the empty-list

cons attaches the first argument to the second argument (which is normally a list)

Now we can build the actual matrix

(defun matrix-zeroes (number-of-rows number-of-columns)
  "Builds a matrix of zeroes"
  (list 
   number-of-rows 
   number-of-columns 
   (build-list-of-zeroes (* number-of-rows number-of-columns))))

list build a list of values equivalent to (cons number-of-rows (cons number-of-columns (cons (build-list-of-zeroes ...) '() ))) (note the last empty-list element)

Similarly we would like to be able to build matrices of random integers

(defun build-list-of-random-digits (size)
  "Recursively builds a list of zeroes of the given size"
  (if (zerop size) 
      '() ;base case
    (cons (random 10) (build-list-of-random-digits(1- size))))) ; iterative step

(defun matrix-random-digits (number-of-rows number-of-columns)
  "Builds a matrix of zeroes"
  (list 
   number-of-rows 
   number-of-columns 
   (build-list-of-random-digits (* number-of-rows number-of-columns))))

; ex: (matrix-random-digits 2 3)

And matrices from manually entered data values

(defun matrix-from-data-list (number-of-rows number-of-columns data-list)
  "Builds a matrix from a data list"
  (list 
   number-of-rows 
   number-of-columns 
   data-list))

A couple of helper function to get the 3 fields will also improve readability

(defun matrix-rows (matrix)
  "Get the number of rows"
  (nth 0 matrix))
(defun matrix-columns (matrix)
  "Get the number of columns"
  (nth 1 matrix))
(defun matrix-data (matrix)
  "Get the data list from the matrix"
  (nth 2 matrix))

nth gets the nth element of the list

For debugging we need to be able to print out the matrix for inspection

(defun matrix-data-get-first-n-values (data n)
  "Given a list of values, get the first n in a string"
  (if (zerop n)
      "" ;base case
    (concat
     (number-to-string (car data))
     " "
     (matrix-data-get-first-n-values (cdr data) (1- n))))) ;iterative step

(defun matrix-data-print (number-of-rows number-of-columns data)
  "Print out the data list gives the dimension of the original matrix"
  (if (zerop number-of-rows)
      "" ;base case
    (concat
     (matrix-data-get-first-n-values data number-of-columns)
     "\n"
     (matrix-data-print ;iterative step
      (1- number-of-rows)
      number-of-columns
      (nthcdr number-of-columns data )))))

(defun matrix-print (matrix)
  "Print out the matrix"
  (concat "\n" (matrix-data-print
                (matrix-rows matrix)
                (matrix-columns matrix)
                (matrix-data matrix))))
; ex:  (message (matrix-print (matrix-from-values 2 2 '(1 2 3 4))))

cdr returns the list without the first element

Matrix Operations

Next we need to define some basics operations on the martices

Addition

The simplest operation is addition. We need to check the matrices have the right size and then simple add the values lists

(defun matrix-equal-size-p (matrix1 matrix2)
  "Check if 2 matrices are the same size"
  (and
   (equal
    (matrix-rows matrix1)
    (matrix-rows matrix2))
   (equal
    (matrix-columns matrix1)
    (matrix-columns matrix2))))
(defun for-each-pair (list1 list2 operator)
  "Go through 2 lists applying an operator on each pair of elements"
  (if (null list1)
      '()
    (cons
     (funcall operator (car list1) (car list2))
     (for-each-pair (cdr list1) (cdr list2) operator))))

(defun matrix-add (matrix1 matrix2)
  "Add to matrices together"
  (if (check-addition matrix1 matrix2)
      (matrix-from-data-list
       (matrix-rows matrix1)
       (matrix-columns matrix1)
       (for-each-pair
        (matrix-data matrix1)
        (matrix-data matrix2)
        '+))))

(defun matrix-subtract (matrix1 matrix2)
  "Subtract MATRIX2 from MATRIX1"
  (if (check-addition matrix1 matrix2)
      (matrix-from-data-list
       (matrix-rows matrix1)
       (matrix-columns matrix1)
       (for-each-pair
        (matrix-data matrix1)
        (matrix-data matrix2)
        '-))))

funcall applied the first arugment (a function) with the remaining items in the list as arguments

Submatrices

The next fundamental step is we want to be able to extract submatrices

(defun matrix-extract-subrow (matrix row start-column end-column)
  "Get part of a row of a matrix and generate a row matrix from it. START-COLUMN is inclusive,  END-COLUMN is exclusive"
  (let
      ((number-of-columns-on-input (matrix-columns matrix))
       (number-of-columns-on-output (-
                                     end-column 
                                     start-column)))
    (matrix-from-data-list
     1
     number-of-columns-on-output
     (subseq
      (matrix-data matrix)
      (+ (* row number-of-columns-on-input) start-column)
      (+ (* row number-of-columns-on-input) end-column)))))

(defun matrix-append (matrix1 matrix2)
  "Append one matrix (set of linear equations) to another"
  (if (null matrix2)
      matrix1
    (matrix-from-data-list
     (+
      (matrix-rows matrix2)
      (matrix-rows matrix1))
     (matrix-columns matrix1)
     (append
      (matrix-data matrix1)
      (matrix-data matrix2)))))

(defun matrix-submatrix (matrix start-row start-column end-row end-column)
  "Get a submatrix. start-row/column are inclusive. end-row/column are exclusive"
  (if (equal start-row end-row)
      '()
    (matrix-append
     (matrix-extract-subrow matrix start-row start-column end-column)
     (matrix-submatrix
      matrix
      (1+ start-row)
      start-column
      end-row
      end-column))))

Which allows us to build these two very familiar function

(defun matrix-get-row (matrix row)
  "Get a row from a matrix. Index starts are ZERO"
  (matrix-extract-subrow
   matrix
   row
   0
   (matrix-columns matrix)))

(defun matrix-get-column (matrix column)
  "Get a column from a matrix. Index starts are ZERO"
  (matrix-submatrix
   matrix
   0
   column
   (nth 0 matrix)
   (1+ column)))

Multiplication

There are fundamentally two different multiplications for matrices. One is multiplying 2 conformable matrices, and the other is multiplying a matrix by a scalar. First I will deal with the latter b/c it's the simpler case.

(defun matrix-scalar-product (matrix scalar)
  "Multiple the matrix by a scalar. ie. multiply each value by the scalar"
  (matrix-from-values
   (matrix-rows matrix)
   (matrix-columns matrix)
   (mapcar
   (lambda (x) 
     (* scalar x))
   (matrix-data matrix))))

To do matrix multiplcation we need to work in small steps and first define the inner product

(defun matrix-inner-product (row column)
  "Multiply a row times a column and returns a scalar"
  (reduce
   '+
   (for-each-pair
    (matrix-data row)
    (matrix-data column)
    '*)))

reduce works down the list elements-by-element applying the operator on each cumulative result

Once we have the inner product we can calculate each value (it's just row * column) and build up to calculating the whole matrix product.

(defun matrix-product-one-value (matrix1 matrix2 row column)
  "Calculate one value in the resulting matrix of the product of two matrices"
  (matrix-inner-product
   (matrix-get-row matrix1 row )
   (matrix-get-column matrix2 column)))

(defun matrix-product-rec (matrix1 matrix2 row column)
  "A recursive helper function that builds the matrix multiplication's data vector"
  (if (equal (matrix-rows matrix1) row)
      '()
    (if (equal (matrix-columns matrix2) column)
        (matrix-product-rec
         matrix1
         matrix2
         (1+ row)
         0)
      (cons
       (matrix-product-one-value
        matrix1
        matrix2
        row column)
       (matrix-product-rec
        matrix1
        matrix2
        row
        (1+ column))))))

(defun matrix-conformable? (matrix1 matrix2)
  "Check that two matrices can be multiplied"
  (equal
   (matrix-columns matrix1)
   (matrix-rows matrix2)))

(defun matrix-product (matrix1 matrix2)
  "Multiply two matrices"
  (matrix-from-data-list
   (matrix-rows matrix1)
   (matrix-columns matrix2)
   (matrix-product-rec
    matrix1
    matrix2
    0
    0)))

Transposition

Transposition is quick and easy using our existing tools

(defun matrix-transpose (matrix)
  "Get the transpose of a matrix"
  (if (equal (matrix-columns matrix) 1)
    (matrix-from-values
     1
     (matrix-rows matrix)
     (matrix-data matrix))
    (matrix-append
     (matrix-from-values
      1
      (matrix-rows matrix)
      (matrix-data (matrix-get-column matrix 0)))
     (matrix-transpose
      (matrix-submatrix
       matrix
       0
       1
       (matrix-rows matrix)
       (matrix-columns matrix))))))

Elementary Matrices

Identity Matrix

The simplest matrix is the identity matrix I. For any matrix A * I = A. This is the matrix of all 1's on the diagonal.

(defun matrix-build-identity-rec (rank row column)
"Helper function that build the data vector of the identity matrix"
(if (equal column rank) ; time to build next row
    (if (equal row (1- rank))
        '() ; we're done
      (matrix-build-identity-rec 
       rank
       (1+ row)
       0))
  (if (equal row column)
      (cons 
       1
       (matrix-build-identity-rec
        rank 
        row 
        (1+ column)))
    (cons 
       0
       (matrix-build-identity-rec
        rank 
        row 
        (1+ column))))))

(defun matrix-identity (rank)
"Build an identity matrix of the given size/rank"
(matrix-from-values rank rank (matrix-build-identity-rec rank 0 0 )))

Unit Column

The next simplest matrix is the unit column

(defun matrix-unit-column (row size)
  "Build a unit column. ROW is where you want the 1 to be placed. SIZE is the overall length"
  (letrec
      ((matrix-unit-column-data-rec
       (lambda (row size)
         (if (equal size 0)
             '()
           (if (equal row 0)
               (cons
                1
                (funcall matrix-unit-column-data-rec
                 (1- row)
                 (1- size)))
             (cons
              0
              (funcall matrix-unit-column-data-rec
               (1- row)
               (1- size))))))))
  (matrix-from-values
   size
   1
   (funcall matrix-unit-column-data-rec
            row
            size))))

Here I'm just trying out a new notation. With letrec we can hide the recursive helper function inside the function that uses it.

Matrix Elementary Transformations

Taking the identity matrix we can change it into any matrix we want. The manipulation of matrices can be broken down into 3 fundamental elementary operation

Type I - Row/Column Interchange

Interchaning rows (or columns) i and j

(defun matrix-elementary-interchange (rowcol1 rowcol2 rank)
  "Make an elementary row/column interchange matrix for ROWCOL1 and ROWCOL2 (ZERO indexed)"
  (let ((u
         (matrix-subtract
          (matrix-unit-column rowcol1 rank)
          (matrix-unit-column rowcol2 rank))))
  (matrix-subtract
   (matrix-identity rank)
   (matrix-product
    u
    (matrix-transpose u)))))

Type II - Row/Column Multiple

Multiplying row (or column) i by α

(defun matrix-elementary-multiply (rowcol scalar rank)
  "Make an elementary row/column multiple matrix for a given ROWCOL (ZERO indexed)"
  (let ((elementary-column
         (matrix-unit-column rowcol rank)))
  (matrix-subtract
   (matrix-identity rank)
   (matrix-product
    elementary-column
    (matrix-scalar-product
     (matrix-transpose elementary-column)
     (- 1 scalar))))))

Type III - Row/Column Addition

Adding a multiple of a row (or column) i to row (or column) j

(defun matrix-elementary-addition (rowcol1 rowcol2 scalar rank)
  "Make an elementary row/column product addition matrix. Multiply ROWCOL1 (ZERO indexed) by SCALAR and add it to ROWCOL2 (ZERO indexed)"
  (matrix-add
   (matrix-identity rank)
   (matrix-scalar-product
    (matrix-product
     (matrix-unit-column rowcol2 rank)
     (matrix-transpose
      (matrix-unit-column rowcol1 rank)))
    scalar)))

TODOs

Remove zero matrix and random values matrix. Makes things longer and confusing and isn't necessary at all

This webpage is generated from an org-document (at ./index.org) that also generates all the files described.

Once opened in Emacs:

C-c C-e h h generates the webpage

C-c C-v C-t exports the code blocks into the appropriate files