#+TITLE: Linear Systems in Lisp
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="https://geokon-gh.github.io/static/worg.css" />
#+options: num:nil
# This will export a README.org file for Github, so that people that land in my repo know where to find the relevant webpage
#+BEGIN_SRC org :tangle README.org :exports none :eval never
  see description [[http://geokon-gh.github.io/linearsystems/index.html][here]]
#+END_SRC

* Intro
I'm going through [[matrixanalysis.com][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
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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
#+END_SRC
#+BEGIN_QUOTE
~zerop~ tests if the value is zero
#+END_QUOTE
#+BEGIN_QUOTE
~()~ with a quote is the /empty-list/ 
#+END_QUOTE
#+BEGIN_QUOTE
~cons~ attaches the first argument to the second argument (which is normally a list)
#+END_QUOTE

Now we can build the actual matrix
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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))))
#+END_SRC
#+BEGIN_QUOTE
~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)*
#+END_QUOTE
Similarly we would like to be able to build matrices of random integers
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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)
#+END_SRC
And matrices from manually entered data values
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el

  (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))
#+END_SRC
A couple of helper function to get the 3 fields will also improve readability
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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))
#+END_SRC
#+BEGIN_QUOTE
~nth~ gets the nth element of the list
#+END_QUOTE
For debugging we need to be able to print out the matrix for inspection
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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))))
#+END_SRC
#+BEGIN_QUOTE
~cdr~ returns the list without the first element
#+END_QUOTE
* 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
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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)
          '-))))
#+END_SRC
#+BEGIN_QUOTE
~funcall~ applied the first arugment (a function) with the remaining items in the list as arguments
#+END_QUOTE
** Submatrices
The next fundamental step is we want to be able to extract submatrices
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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))))

#+END_SRC
Which allows us to build these two very familiar function
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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)))
#+END_SRC

** 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.
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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))))
#+END_SRC
To do matrix multiplcation we need to work in small steps and first define the inner product
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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)
      '*)))
#+END_SRC
#+BEGIN_QUOTE
~reduce~ works down the list elements-by-element applying the operator on each cumulative result
#+END_QUOTE
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.
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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)))
#+END_SRC

** Transposition
Transposition is quick and easy using our existing tools
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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))))))
#+END_SRC
** TESTS :noexport:
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (ert-deftest matrix-test-operations ()
    "Testing - Matrix Operations"
    (let ((matrix1 '(2 2 (1 2 3 4)))
          (matrix2 '(2 2 (5 6 7 8))))
    (should (equal
             (matrix-extract-subrow '(2 2 (1 2 3 4)) 1 0 2)
             '(1 2 (3 4))))
    (should (equal
             (matrix-equal-size-p matrix1 matrix2)
             't))
    (should (equal
             (matrix-add matrix1 matrix2)
             '(2 2 (6 8 10 12))))
    (should (equal
             (matrix-subtract matrix1 matrix2)
             '(2 2 (-4 -4 -4 -4))))
    (should (equal
             (matrix-scalar-product
              (matrix-identity 3)
              7)
             '(3 3 (7 0 0 0 7 0 0 0 7))))))

#+END_SRC
* 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.
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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 )))
#+END_SRC
** Unit Column
The next simplest matrix is the unit column
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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))))
#+END_SRC
#+BEGIN_QUOTE
Here I'm just trying out a new notation. With ~letrec~ we can hide the recursive helper function inside the function that uses it.
#+END_QUOTE
** 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/
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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)))))
#+END_SRC

*** Type II - Row/Column Multiple
Multiplying row (or column) /i/ by /\alpha/
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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))))))
#+END_SRC

*** Type III - Row/Column Addition
Adding a multiple of a row (or column) /i/ to row (or column) /j/
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (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)))
#+END_SRC

** TESTS :noexport:
#+BEGIN_SRC emacs-lisp :results output silent :session :tangle matrix.el
  (ert-deftest matrix-test-elementary-operation ()
    "Testing - Elementary Matrix Transformations"
    (let ((matrix1 '(2 2 (1 2 3 4)))
          (matrix2 '(2 2 (5 6 7 8))))
      (should (equal
               (matrix-identity 3)
               '(3 3 (1 0 0 0 1 0 0 0 1))))
      (should (equal
               (matrix-unit-column 3 5)
               '( 5 1 (0 0 0 1 0))))
      (should (equal
               (matrix-elementary-interchange 0 1 3)
               '(3 3 (0 1 0 1 0 0 0 0 1))))
      (should (equal
               (matrix-elementary-multiply 1 7 3)
               '(3 3 (1 0 0 0 7 0 0 0 1))))
      (should (equal
               (matrix-elementary-addition 0 2 7 3)
               '(3 3 (1 0 0 0 1 0 7 0 1))))))
#+END_SRC

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

#+BEGIN_QUOTE
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\\
#+END_QUOTE