Re: [Sbcl-devel] [Sbcl-bugs] Help requested in the definition of matrix-multiply for n-dimensional matrices. Bugs I do expect! SBCL times reported.

[Stas Boukarev]

Using specialized arrays instead of lists would produce better results.

On Fri, Oct 6, 2023 at 10:34 PM Robert Boyer <robertste...@gmail.com> wrote:

;;; -*- coding: utf-8 -*-

;;; Coded in October 2023.  This software is in the public domain.

;;; Below please find two possibly correct definitions of the addition and
;;; multiplication of n-dimensional matrices, but I am far from sure.  So
;;; please send bug reports to robertste...@gmail.com.

;;; 'AI chips', or so I have heard, do multiplication of three dimensional
;;; matrices very fast.

;;; I have no idea what else AI chips do, but I hope to find out, and would
;;; be most grateful for information on this topic.

;;; What is the definition of the multiplication of two n-dimensional
;;; matrices?  Below is the best answer that I have come up with, but do not
;;; count on its being correct!

;;; Lisp has so much to offer. I define matrix-p thus:

(in-package "COMMON-LISP-USER")

(declaim
 (optimize (speed 3) (safety 0) (debug 0)))

(defun matrix-p (m)
  "(matrix-p m) is a recognizer for n-dimensional matrices built from
  numbers, conses, and NIL."
  (cond ((numberp m) t)
        (t (and (consp m)
                (null (cdr (last m)))
                (loop for mi in m always (matrix-p mi))
                (or
                 (loop for mi in m always (numberp mi))
                 (loop for mi in m always
                       (and (consp mi)
                            (= (length (car m))
                               (length mi)))))))))

(defun random-n-dimensional-matrix-1 (n i r)
  (cond ((= i 1)
         (loop for j from 1 to n collect
               (cond ((complexp r)
                      (complex (random (realpart 4))
                               (random (imagpart r))))
                     (t (random r)))))
        (t
         (loop for j from 1 to n collect
                 (random-n-dimensional-matrix-1
                  n
                  (1- i)
                  r)))))

(defun random-n-dimensional-matrix (n r)
  "(random-n-dimensional-matrix n r) returns a n-dimensional matrix that is
   seeded with random values based upon r."
  (random-n-dimensional-matrix-1 n n r))

;;; foo1, foo2, and foo3 are examples of three dimensional
;;; matrices.  They are printed below.

(defparameter example-n 3)

(defparameter foo1 (random-n-dimensional-matrix example-n 1.0))

(defparameter foo2 (random-n-dimensional-matrix example-n 1.0))

(defparameter foo3 (random-n-dimensional-matrix example-n  #c(1 2)))

;; As a warm up, here is the definition of matrix-add for two n-dimensional
;; matrices.

(defun matrix-add (a b)
  "(matrix add a b) adds two matricies of the same shape."
  (cond ((numberp a)
         (cond ((numberp b)  (+ a b))
               (t (error "~%matrix-add: ~s and ~s must have the same shape."
                         a b))))
        ((or (not (consp a))
             (not (consp b))
             (not (= (length a) (length b))))
         (error "~%matrix-add: ~s and ~s must have the same shape." a b))
        (t (loop for ai in a as bi in b
                 collect (matrix-add ai bi)))))

#|

* foo1
(((0.62944734 0.2709539 0.81158376) (0.6700171 0.2539736 0.93773544)
(0.8267517 0.44206798 0.2647184))
((0.6163341 0.19508076 0.094441175) (0.55699635 0.37676394 0.093762994)
(0.9857626 0.91501355 0.99292254))
((0.929777 0.93538976 0.31522608) (0.45167792 0.9411855 0.96221936)
(0.9143338 0.21972346 0.9707513)))
* foo2
(((0.5962117 0.6005609 0.5940589) (0.2837726 0.009566903 0.84352255)
(0.22492898 0.83147097 0.2795267))
((0.5844146 0.7568612 0.9189848) (0.0073252916 0.31148136 0.59585714)
(0.07142329 0.72258794 0.6982585))
((0.42384863 0.86798644 0.36271906) (0.35747027 0.797477 0.51548016)
(0.4812944 0.48626482 0.94951725)))

(matrix-add foo1 foo2)
(((1.225659 0.8715148 1.4056426) (0.9537897 0.2635405 1.781258)
(1.0516807 1.273539 0.5442451))
((1.2007487 0.95194197 1.013426) (0.56432164 0.6882453 0.68962014)
(1.0571859 1.6376015 1.6911811))
((1.3536257 1.8033762 0.67794514) (0.8091482 1.7386625 1.4776995)
(1.3956282 0.7059883 1.9202685)))

So matrix-add might be ok!

|#

;;; Here are the definitions of 'row' and 'column' for use in the definition
;;; of 'matrix-multiply'.

(defun row (i m)
  (nth i m))

(defun column (i m)
  (loop for ai in m collect (nth i ai)))

;;; Definition of matrix-multiply for two n-dimensional matrices.

(defun matrix-multiply (a b)
  "(matric-multiply a b) multiplies two n-dimensional matrices of the same
   n-dimension.  'All the way down' the number of columns of the first must be
   the number of rows of the second b, a universal requirement."
  (cond ((numberp a) (* a b))
        ((numberp b)
         (error "~%matrix-multiply: unexpected call with a = ~s and b = ~s.~%" a b))
        ((and (numberp (car a)) (numberp (car b)))
         (loop for ai in a as bi in b collect (* ai bi)))
        ((not (= (length (car a)) (length b)))
         (error
          "~%matrix-multiply: # of columns of a = ~s must = # of rows of b = ~s." a b))
        (t (loop for i below (length a)
                 collect
                 (matrix-multiply (row i a) (column i b))))))

#|

(matrix-multiply foo1 foo2)
(((0.37528384 0.15834941 0.34398866) (0.4023861 0.19222277 0.81394166)
  (0.49113917 0.40625376 0.09601841))
 ((0.17489871 0.0014290235 0.033759914) (0.00532873 0.117354944 0.07477383)
  (0.831513 0.54521734 0.5118319))
 ((0.20913379 0.06680862 0.15171655) (0.37555707 0.6800893 0.46789342)
  (0.25558072 0.15342379 0.92174506)))

Here is the best info about matrix multiply that I was able to find:

If  the matrix of T is m×n with entries  a(i, j) (let's call it A),
and the matrix of S is n×p with entries  b(r, s) (let's call it B),
then the matrix of T∘S (let's call it A∘B or AB) is m×p
  with entries c(k, ℓ), where
    c(k, ℓ) =
   (a(k, 1) * b(1, ℓ)) + (a(k, 2) * b(2, ℓ)) + ⋯ + (a(k, n) * b(n, ℓ)).

How well does SBCL do time-wise?

|#

(time (loop for i from 1 to 1000000 do (matrix-multiply foo1 foo2)))

#|

Evaluation took:
  1.997 seconds of real time
  1.995670 seconds of total run time (1.911724 user, 0.083946 system)
  [ Run times consist of 0.031 seconds GC time, and 1.965 seconds non-GC time. ]
  99.95% CPU
  2,185,467,707 processor cycles
  1,199,991,696 bytes consed

Is 2 millionths of a second on a $100 Chromebook for a 3-d matrix multiply
good?  I have no idea.

I don't really know what I am doing, so I would appreciate example matrices
and answers to try tests on.  I searched for hours on the Internet and could
not find any that I could successfully copy and paste.

Multiplying two 7-dimensional matrices 100 times took 7.5 seconds, so
that is 7.5 hundredths of a second per multiply.

Getting past 7 seems to take up more core than I have on my Chromebook.

|#

Thanks so very much for the wonderful SBCL and for the wonderful bug support.

Bob

_______________________________________________
Sbcl-bugs mailing list
Sbcl...@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/sbcl-bugs