Lisp functions for elementary vector and matrix products

;; Lisp functions for elementary vector and matrix products

;; (dot U V) == dot product of two vectors == U.V
;; (transpose A) == transpose matrix A == A^T
;; (mult A B) === multiply two matrices AxB
;; (outer f U V) == outer product of vectors U and V wrt function f
;;   == matrix of elements [f(U[i], V[j])]
;; (inner f g U V) == generalized inner product of vectors U and V
;;   == fold(f, [g(U[i],V[i])])
;; (poly U x) == evaluate polynomial w/ coefficients U at point x

;; The goal here is to define these operations as simply as possible,
;; in order to demonstrate basic principles,
;; so we don't do any input validation or error catching.

;; N.B. to use the following examples you must define RANGE and IOTA
;; as described here: http://pastebin.com/h1dNCaiQ
;; If this file is called range.lsp then do this:
(load "range.lsp")

;; vector is just a list
;; example vectors
(defvar U '(1 3 -5))
(defvar V '(4 -2 -1))

;; dot product of two vectors U.V
(defun dot (U V) (apply #'+ (mapcar #'* U V)))

;; example dot product
(dot U V) ; U.V == 3

;; matrix is a list of lists
;; example matrices
(defvar A '((1 3 -5) (4 -2 -1) (0 7 2)))
(defvar B '((-3 2 5) (1 3 -7) (-1 -2 8)))

;; transpose matrix: exchange rows and columns
(defun transpose (M) (apply #'mapcar #'list M))

;; matrix multiplication
(defun mult (A B)
  (mapcar (lambda (row)
    (mapcar (lambda (col)
      (dot row col))
    (transpose B)))
  A))

;; example matrix multiplication
(mult A B) ;; AxB == ((5 21 -56) (-13 4 26) (5 17 -33))

;; remember, vector is list and matrix is list of lists
;; so to mult vector by matrix, nest vector in parens like this...
(mult (list U) B) ; == 1x3 matrix: ((5 21 -56))
(mult A (transpose (list V))) ; == 3x1 matrix: ((3) (21) (-16))

;; dot product U.V == UxV^T
;; except we'll get nested parens, b/c this is matrix, not vector
(mult (list U) (transpose (list V))) ; == ((3))

;; outer product of two vectors
(defun outer (fn U V)
  (mapcar (lambda (x)
    (mapcar (lambda (y)
      (funcall fn x y))
    V))
  U))

;; example outer product
(outer #'* (iota 12) (iota 12)) ; 1-12 times tables

;; inner product of two vectors
(defun inner (f g U V) (apply f (mapcar g U V)))

;; same as (dot U V)
(inner #'+ #'* U V) ; == 3

;; (1^1 2^2 3^3) == (1 4 27)
(inner #'list #'expt (iota 3) (iota 3))

;; 1^1 + 2^2 + 3^3 == 32
(inner #'+ #'expt (iota 3) (iota 3))

;; 1^1 * 2^2 * 3^3 == 108
(inner #'* #'expt (iota 3) (iota 3))

;; polynomial evaluation
;; coefficient list a == (a_0 a_1 a_2 ... a_n)
;; point to eval at == x
;; want: a_0 + a_1 * x + a_2 * x^2 + ... + a_n * x^n
(defun poly (a x)
  (dot a (inner #'list #'expt (make-list (length a) :initial-element x) (range 0 (1- (length a))))))

;; 4*(-3)^0 - 2*(-3)^1 - 1*(-3)^2 == 1
(poly V -3)

;; 1*(1/2)^0 + 3*(1/2)^1 - 5*(1/2)^2 == 1.25 == 5/4
(poly U 1/2)