matrices.lsp

   1 ;;;; matrices -- Basic matrix operations
   2 ;;;;
   3 ;;;; Copyright (c) 1991, by Luke Tierney. Permission is granted for
   4 ;;;; unrestricted use.
   5
   6 (provide "matrices")
   7
   8 ;;;;
   9 ;;;; Package Setup
  10 ;;;;
  11
  12 #+:CLtL2
  13 (in-package lisp-stat-basics)
  14 #-:CLtL2
  15 (in-package 'lisp-stat-basics)
  16
  17 (export '(matrixp num-rows num-cols matmult identity-matrix diagonal
  18           row-list column-list inner-product outer-product cross-product
  19           transpose bind-columns bind-rows))
  20
  21 (defun matrixp (x)
  22 "Args: (x)
  23 Returns T if X is a matrix, NIL otherwise."
  24   (and (arrayp x) (= (array-rank x) 2)))
  25
  26 (defun num-rows (x)
  27 "Args: (x)
  28 Returns number of rows in X."
  29   (array-dimension x 0))
  30
  31 (defun num-cols (x)
  32 "Args: (x)
  33 Returns number of columns in X."
  34   (array-dimension x 1))
  35
  36 (defun matmult (a b &rest args)
  37 "Args: (a b &rest args)
  38 Returns the matrix product of matrices a, b, etc. If a is a vector it is
  39 treated as a row vector; if b is a vector it is treated as a column vector."
  40   (let ((rtype (cond ((and (matrixp a) (matrixp b)) 'matrix)
  41                      ((and (sequencep a) (sequencep b)) 'number)
  42                      ((sequencep a) (if (consp a) 'list 'vector))
  43                      ((sequencep b) (if (consp b) 'list 'vector)))))
  44
  45     (if (sequencep a)
  46       (setf a (vector-to-array (coerce a 'vector) (list 1 (length a)))))
  47     (if (sequencep b)
  48       (setf b (vector-to-array (coerce b 'vector) (list (length b) 1))))
  49     (if (not (= (array-dimension a 1) (array-dimension b 0)))
  50       (error "dimensions do not match"))
  51     (if args
  52       (reduce #'matmult args :initial-value (matmult a b))
  53       (let* ((n (array-dimension a 0))
  54              (m (array-dimension b 1))
  55              (p (array-dimension a 1))
  56              (c (make-array (list n m)))
  57              x)
  58         (declare (fixnum n m p))
  59         (dotimes (i n)
  60                  (declare (fixnum i))
  61                  (dotimes (j m)
  62                           (declare (fixnum j))
  63                           (setq x 0)
  64                           (dotimes (k p)
  65                                    (declare (fixnum k))
  66                                    (setq x (+ x
  67                                               (* (aref a i k) (aref b k j)))))
  68                           (setf (aref c i j) x)))
  69         (case rtype
  70           (matrix c)
  71           (number (aref c 0 0))
  72           (t (coerce (compound-data-seq c) rtype)))))))
  73
  74 (defun identity-matrix (n)
  75 "Args: (n)
  76 Returns the identity matrix of rank N."
  77   (let ((result (make-array (list n n) :initial-element 0)))
  78     (dotimes (i n result)
  79              (declare (fixnum i))
  80              (setf (aref result i i) 1))))
  81
  82 ;; this thing is not very efficient at this point - too much coercing
  83 (defun diagonal (x)
  84 "Args: (x)
  85 If X is a matrix, returns the diagonal of X. If X is a sequence, returns a
  86 diagonal matrix of rank (length X) with diagonal elements eq to the elements
  87 of X."
  88   (cond ((matrixp x)
  89          (let* ((n (min (num-rows x) (num-cols x)))
  90                 (result (make-array n)))
  91            (dotimes (i n (coerce result 'list))
  92                     (setf (aref result i) (aref x i i)))))
  93         ((sequencep x)
  94          (let* ((x (coerce x 'vector))
  95                 (n (length x))
  96                 (result (make-array (list n n) :initial-element 0)))
  97            (dotimes (i n result)
  98                     (setf (aref result i i) (aref x i)))))
  99         (t (error "argument must be a matrix or a sequence"))))
 100
 101 (defun row-list (x)
 102 "Args: (m)
 103 Returns a list of the rows of M as vectors"
 104   (check-matrix x)
 105   (let ((m (num-rows x))
 106         (n (num-cols x))
 107         (result nil))
 108     (declare (fixnum m n))
 109     (flet ((get-row (k)
 110              (declare (fixnum k))
 111              (let ((row (make-array n)))
 112                (dotimes (i n row)
 113                  (declare (fixnum i))
 114                  (setf (aref row i) (aref x k i))))))
 115       (dotimes (i m result)
 116         (declare (fixnum i))
 117         (setf result (cons (get-row (- m i 1)) result))))))
 118
 119 (defun column-list (x)
 120 "Args: (m)
 121 Returns a list of the columns of M as vectors"
 122   (check-matrix x)
 123   (let ((m (num-rows x))
 124         (n (num-cols x))
 125         (result nil))
 126     (declare (fixnum m n))
 127     (flet ((get-col (k)
 128              (declare (fixnum k))
 129              (let ((col (make-array m)))
 130                (dotimes (i m col)
 131                  (declare (fixnum i))
 132                  (setf (aref col i) (aref x i k))))))
 133       (dotimes (i n result)
 134         (declare (fixnum i))
 135         (setf result (cons (get-col (- n i 1)) result))))))
 136
 137 (defun inner-product (x y)
 138 "Args: (x y)
 139 Returns inner product of sequences X and Y."
 140   (check-sequence x)
 141   (check-sequence y)
 142   (let ((n (length x))
 143         (cx (make-next-element x))
 144         (cy (make-next-element y))
 145         (result 0))
 146     (declare (fixnum n))
 147     (if (/= n (length y)) (error "sequence lengths do not match"))
 148     (dotimes (i n result)
 149       (declare (fixnum i))
 150       (setf result
 151             (+ result (* (get-next-element cx i) (get-next-element cy i)))))))
 152
 153 (defun outer-product (x y &optional (f #'*))
 154 "Args: (x y &optional (fcn #'*))
 155 Returns the generalized outer product of x and y, using fcn. Tat is, the result
 156 is a matrix of dimension ((length x) (length y)) and the (i j) element of the
 157 result is computed as (apply fcn (aref x i) (aref y j))."
 158   (let* ((x (coerce x 'vector))
 159          (y (coerce y 'vector))
 160          (m (length x))
 161          (n (length y))
 162          (a (make-array (list m n))))
 163     (declare (fixnum m n))
 164     (dotimes (i m a)
 165       (declare (fixnum i))
 166       (dotimes (j n)
 167         (declare (fixnum j))
 168         (setf (aref a i j) (funcall f (aref x i) (aref y j)))))))
 169
 170 (defun cross-product (x)
 171 "Args: (x)
 172 If X is a matrix returns (matmult (transpose X) X). If X is a vector returns
 173 (inner-product X X)."
 174   (check-matrix x)
 175   (let* ((n (num-rows x))
 176          (p (num-cols x))
 177          (c (make-array (list p p))))
 178     (declare (fixnum n p))
 179     (dotimes (i p c)
 180       (declare (fixnum i))
 181       (dotimes (j (+ i 1))
 182         (declare (fixnum j))
 183         (let ((val 0))
 184           (dotimes (k n)
 185             (declare (fixnum k))
 186             (incf val (* (aref x k i) (aref x k j))))
 187           (setf (aref c i j) val)
 188           (setf (aref c j i) val))))))
 189
 190 (defun transpose-list (x)
 191   (let ((m (length (first x))))
 192     (dolist (next x)
 193       (if (not (consp next)) (error "not a list - ~a" x))
 194       (if (/= m (length next)) (error "sublists not the same length")))
 195     (do* ((cx (copy-list x))
 196           (result (make-list m))
 197           (next result (cdr next)))
 198          ((null next) result)
 199       (setf (first next) (mapcar #'first cx))
 200       (do ((next cx (cdr next)))
 201           ((null next))
 202         (setf (first next) (rest (first next)))))))
 203
 204 (defun transpose (x)
 205 "Args: (m)
 206 Returns the transpose of the matrix M."
 207   (cond
 208    ((consp x) (transpose-list x))
 209    (t
 210     (check-matrix x)
 211     (let* ((m (num-rows x))
 212            (n (num-cols x))
 213            (tx (make-array (list n m))))
 214       (declare (fixnum m n))
 215       (dotimes (i m tx)
 216         (declare (fixnum i))
 217         (dotimes (j n)
 218           (declare (fixnum j))
 219           (setf (aref tx j i) (aref x i j))))))))
 220
 221 (defun bind-columns (&rest args)
 222 "Args (&rest args)
 223 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
 224 along their columns.
 225 Example: (bind-columns #2a((1 2)(3 4)) #(5 6)) returns #2a((1 2 5)(3 4 6))"
 226   (flet ((check-arg (x)
 227            (if (not (or (sequencep x) (matrixp x)))
 228                (error "bad argument type")))
 229          (arg-cols (x) (if (sequencep x) 1 (num-cols x)))
 230          (arg-rows (x) (if (sequencep x) (length x) (num-rows x))))
 231     (dolist (x args) (check-arg x))
 232     (let ((m (arg-rows (first args)))
 233           (n (arg-cols (first args))))
 234       (declare (fixnum m n))
 235       (dolist (x (rest args))
 236         (if (/= m (arg-rows x)) (error "column lengths do not match"))
 237         (incf n (arg-cols x)))
 238       (do* ((result (make-array (list m n)))
 239             (args args (rest args))
 240             (firstcol 0)
 241             (x (first args) (first args)))
 242            ((null args) result)
 243         (cond
 244          ((sequencep x)
 245           (let ((cx (make-next-element x)))
 246             (dotimes (i m)
 247               (setf (aref result i firstcol) (get-next-element cx i)))))
 248          (t
 249           (let ((k (arg-cols x)))
 250             (dotimes (i m)
 251               (dotimes (j k)
 252                 (setf (aref result i (+ firstcol j)) (aref x i j)))))))
 253         (incf firstcol (arg-cols x))))))
 254
 255 (defun bind-rows (&rest args)
 256 "Args (&rest args)
 257 The ARGS can be matrices, vectors, or lists. Arguments are bound into a matrix
 258 along their rows.
 259 Example: (bind-rows #2a((1 2)(3 4)) #(5 6)) returns #2a((1 2)(3 4)(5 6))"
 260   (flet ((check-arg (x)
 261            (if (not (or (sequencep x) (matrixp x)))
 262                (error "bad argument type")))
 263          (arg-cols (x) (if (sequencep x) (length x) (num-cols x)))
 264          (arg-rows (x) (if (sequencep x) 1 (num-rows x))))
 265     (dolist (x args) (check-arg x))
 266     (let ((m (arg-rows (first args)))
 267           (n (arg-cols (first args))))
 268       (declare (fixnum m n))
 269       (dolist (x (rest args))
 270         (if (/= n (arg-cols x)) (error "row lengths do not match"))
 271         (incf m (arg-rows x)))
 272       (do* ((result (make-array (list m n)))
 273             (args args (rest args))
 274             (firstrow 0)
 275             (x (first args) (first args)))
 276            ((null args) result)
 277         (cond
 278          ((sequencep x)
 279           (let ((cx (make-next-element x)))
 280             (dotimes (i n)
 281               (setf (aref result firstrow i) (get-next-element cx i)))))
 282          (t
 283           (let ((k (arg-rows x)))
 284             (dotimes (i n)
 285               (dotimes (j k)
 286                 (setf (aref result (+ firstrow j) i) (aref x j i)))))))
 287         (incf firstrow (arg-rows x))))))