The Common Lisp Cookbook – Multidimensional arrays

Table of Contents

The Common Lisp Cookbook – Multidimensional arrays

📢 New videos: web dev demo part 1, dynamic page with HTMX, Weblocks demo

📕 Get the EPUB and PDF

Common Lisp has native support for multidimensional arrays, with some special treatment for 1-D arrays, called vectors. Arrays can be generalised and contain any type (element-type t), or they can be specialised to contain specific types such as single-float or integer. A good place to start is Practical Common Lisp Chapter 11, Collections by Peter Seibel.

A quick reference to some common operations on arrays is given in the section on Arrays and vectors.

Some libraries available on Quicklisp for manipulating arrays:

This page covers what can be done with the built-in multidimensional arrays, but there are limitations. In particular:

Both of these problems can be solved by using CLOS to define an extended array class, with native arrays as a special case. Some libraries available through quicklisp which take this approach are:

A relatively new but actively developed package is MAGICL, which provides wrappers around BLAS and LAPACK libraries. At the time of writing this package is not on Quicklisp, and only works under SBCL and CCL. It seems to be particularly focused on complex arrays, but not exclusively. To install, clone the repository in your quicklisp local-projects directory e.g. under Linux/Unix:

$ cd ~/quicklisp/local-projects
$ git clone

Instructions for installing dependencies (BLAS, LAPACK and Expokit) are given on the github web pages. Low-level routines wrap foreign functions, so have the Fortran names e.g magicl.lapack-cffi::%zgetrf. Higher-level interfaces to some of these functions also exist, see the source directory and documentation.

Taking this further, domain specific languages have been built on Common Lisp, which can be used for numerical calculations with arrays. At the time of writing the most widely used and supported of these are:

CLASP is a project which aims to ease interoperability of Common Lisp with other languages (particularly C++), by using LLVM. One of the main applications of this project is to numerical/scientific computing.


The function CLHS: make-array can create arrays filled with a single value

* (defparameter *my-array* (make-array '(3 2) :initial-element 1.0))
* *my-array*
#2A((1.0 1.0) (1.0 1.0) (1.0 1.0))

More complicated array values can be generated by first making an array, and then iterating over the elements to fill in the values (see section below on element access).

The array-operations library provides generate, a convenient function for creating arrays which wraps this iteration.

* (ql:quickload :array-operations)
To load "array-operations":
  Load 1 ASDF system:
; Loading "array-operations"


* (aops:generate #'identity 7 :position)
#(0 1 2 3 4 5 6)

Note that the nickname for array-operations is aops. The generate function can also iterate over the array subscripts by passing the key :subscripts. See the Array Operations manual on generate for more examples.

Random numbers

To create an 3x3 array containing random numbers drawn from a uniform distribution, generate can be used to call the CL random function:

* (aops:generate (lambda () (random 1.0)) '(3 3))
#2A((0.99292254 0.929777 0.93538976)
    (0.31522608 0.45167792 0.9411855)
    (0.96221936 0.9143338 0.21972346))

An array of Gaussian (normal) random numbers with mean of zero and standard deviation of one, using the alexandria package:

* (ql:quickload :alexandria)
To load "alexandria":
  Load 1 ASDF system:
; Loading "alexandria"


* (aops:generate #'alexandria:gaussian-random 4)
#(0.5522547885338768d0 -1.2564808468164517d0 0.9488161476129733d0

Note that this is not particularly efficient: It requires a function call for each element, and although gaussian-random returns two random numbers, only one of them is used.

For more efficient implementations, and a wider range of probability distributions, there are packages available on Quicklisp. See CLiki for a list.

Accessing elements

To access the individual elements of an array there are the aref and row-major-aref functions.

The aref function takes the same number of index arguments as the array has dimensions. Indexing is from 0 and row-major as in C, but not Fortran.

* (defparameter *a* #(1 2 3 4))
* (aref *a* 0)
* (aref *a* 3)
* (defparameter *b* #2A((1 2 3) (4 5 6)))
* (aref *b* 1 0)
* (aref *b* 0 2)

The range of these indices can be found using array-dimensions:

* (array-dimensions *a*)
* (array-dimensions *b*)
(2 3)

or the rank of the array can be found, and then the size of each dimension queried:

* (array-rank *a*)
* (array-dimension *a* 0)
* (array-rank *b*)
* (array-dimension *b* 0)
* (array-dimension *b* 1)

To loop over an array nested loops can be used, such as:

* (defparameter a #2A((1 2 3) (4 5 6)))
* (destructuring-bind (n m) (array-dimensions a)
    (loop for i from 0 below n do
      (loop for j from 0 below m do
        (format t "a[~a ~a] = ~a~%" i j (aref a i j)))))

a[0 0] = 1
a[0 1] = 2
a[0 2] = 3
a[1 0] = 4
a[1 1] = 5
a[1 2] = 6

A utility macro which does this for multiple dimensions is nested-loop:

(defmacro nested-loop (syms dimensions &body body)
  "Iterates over a multidimensional range of indices.

   SYMS must be a list of symbols, with the first symbol
   corresponding to the outermost loop.

   DIMENSIONS will be evaluated, and must be a list of
   dimension sizes, of the same length as SYMS.

    (nested-loop (i j) '(10 20) (format t '~a ~a~%' i j))

  (unless syms (return-from nested-loop `(progn ,@body))) ; No symbols

  ;; Generate gensyms for dimension sizes
  (let* ((rank (length syms))
         ;; reverse our symbols list,
         ;; since we start from the innermost.
         (syms-rev (reverse syms))
         ;; innermost dimension first:
         (dims-rev (loop for i from 0 below rank
                         collecting (gensym)))
         ;; start with innermost expression
         (result `(progn ,@body)))
    ;; Wrap previous result inside a loop for each dimension
    (loop for sym in syms-rev for dim in dims-rev do
         (unless (symbolp sym)
           (error "~S is not a symbol. First argument to nested-loop must be a list of symbols" sym))
         (setf result
               `(loop for ,sym from 0 below ,dim do
    ;; Add checking of rank and dimension types,
    ;; and get dimensions into gensym list.
    (let ((dims (gensym)))
      `(let ((,dims ,dimensions))
         (unless (= (length ,dims) ,rank)
           (error "Incorrect number of dimensions: Expected ~a but got ~a" ,rank (length ,dims)))
         (dolist (dim ,dims)
           (unless (integerp dim)
             (error "Dimensions must be integers: ~S" dim)))
         ;; dimensions reversed so that innermost is last:
         (destructuring-bind ,(reverse dims-rev) ,dims

so that the contents of a 2D array can be printed using:

* (defparameter a #2A((1 2 3) (4 5 6)))
* (nested-loop (i j) (array-dimensions a)
      (format t "a[~a ~a] = ~a~%" i j (aref a i j)))

a[0 0] = 1
a[0 1] = 2
a[0 2] = 3
a[1 0] = 4
a[1 1] = 5
a[1 2] = 6

[Note: This macro is available in this fork of array-operations, but not Quicklisp]

Row major indexing

In some cases, particularly element-wise operations, the number of dimensions does not matter. To write code which is independent of the number of dimensions, array element access can be done using a single flattened index via row-major-aref. The array size is given by array-total-size, with the flattened index starting at 0.

* (defparameter a #2A((1 2 3) (4 5 6)))
* (array-total-size a)
* (loop for i from 0 below (array-total-size a) do
     (setf (row-major-aref a i) (+ 2.0 (row-major-aref a i))))
* a
#2A((3.0 4.0 5.0) (6.0 7.0 8.0))

Infix syntax

The cmu-infix library provides some different syntax which can make mathematical expressions easier to read:

* (ql:quickload :cmu-infix)
To load "cmu-infix":
  Load 1 ASDF system:
; Loading "cmu-infix"


* (named-readtables:in-readtable cmu-infix:syntax)

* (defparameter arr (make-array '(3 2) :initial-element 1.0))

* #i(arr[0 1] = 2.0)

* arr
#2A((1.0 2.0) (1.0 1.0) (1.0 1.0))

A matrix-matrix multiply operation can be implemented as:

(let ((A #2A((1 2) (3 4)))
      (B #2A((5 6) (7 8)))
      (result (make-array '(2 2) :initial-element 0.0)))

     (loop for i from 0 to 1 do
           (loop for j from 0 to 1 do
                 (loop for k from 0 to 1 do
                       #i(result[i j] += A[i k] * B[k j]))))

See the section below on linear algebra, for alternative matrix-multiply implementations.

Element-wise operations

To multiply two arrays of numbers of the same size, pass a function to each in the array-operations library:

* (aops:each #'* #(1 2 3) #(2 3 4))
#(2 6 12)

For improved efficiency there is the aops:each* function, which takes a type as first argument to specialise the result array.

To add a constant to all elements of an array:

* (defparameter *a* #(1 2 3 4))
* (aops:each (lambda (it) (+ 42 it)) *a*)
#(43 44 45 46)
* *a*
#(1 2 3 4)

Note that each is not destructive, but makes a new array. All arguments to each must be arrays of the same size, so (aops:each #'+ 42 *a*) is not valid.

Vectorising expressions

An alternative approach to the each function above, is to use a macro to iterate over all elements of an array:

(defmacro vectorize (variables &body body)
  ;; Check that variables is a list of only symbols
  (dolist (var variables)
    (if (not (symbolp var))
        (error "~S is not a symbol" var)))

    ;; Get the size of the first variable, and create a new array
    ;; of the same type for the result
    `(let ((size (array-total-size ,(first variables)))  ; Total array size (same for all variables)
           (result (make-array (array-dimensions ,(first variables)) ; Returned array
                               :element-type (array-element-type ,(first variables)))))
       ;; Check that all variables have the same sizeo
       ,@(mapcar (lambda (var) `(if (not (equal (array-dimensions ,(first variables))
                                                (array-dimensions ,var)))
                                    (error "~S and ~S have different dimensions" ',(first variables) ',var)))
              (rest variables))

       (dotimes (indx size)
         ;; Locally redefine variables to be scalars at a given index
         (let ,(mapcar (lambda (var) (list var `(row-major-aref ,var indx))) variables)
           ;; User-supplied function body now evaluated for each index in turn
           (setf (row-major-aref result indx) (progn ,@body))))

[Note: Expanded versions of this macro are available in this fork of array-operations, but not Quicklisp]

This can be used as:

* (defparameter *a* #(1 2 3 4))
* (vectorize (*a*) (* 2 *a*))
#(2 4 6 8)

Inside the body of the expression (second form in vectorize expression) the symbol *a* is bound to a single element. This means that the built-in mathematical functions can be used:

* (defparameter a #(1 2 3 4))
* (defparameter b #(2 3 4 5))
* (vectorize (a b) (* a (sin b)))
#(0.9092974 0.28224 -2.2704074 -3.8356972)

and combined with cmu-infix:

* (vectorize (a b) #i(a * sin(b)) )
#(0.9092974 0.28224 -2.2704074 -3.8356972)

Calling BLAS

Several packages provide wrappers around BLAS, for fast matrix manipulation.

The lla package in quicklisp includes calls to some functions:

Scale an array

scaling by a constant factor:

* (defparameter a #(1 2 3))
* (lla:scal! 2.0 a)
* a
#(2.0d0 4.0d0 6.0d0)


This calculates a * x + y where a is a constant, x and y are arrays. The lla:axpy! function is destructive, modifying the last argument (y).

* (defparameter x #(1 2 3))
* (defparameter y #(2 3 4))
* (lla:axpy! 0.5 x y)
#(2.5d0 4.0d0 5.5d0)
* x
#(1.0d0 2.0d0 3.0d0)
* y
#(2.5d0 4.0d0 5.5d0)

If the y array is complex, then this operation calls the complex number versions of these operators:

* (defparameter x #(1 2 3))
* (defparameter y (make-array 3 :element-type '(complex double-float)
                                :initial-element #C(1d0 1d0)))
* y
#(#C(1.0d0 1.0d0) #C(1.0d0 1.0d0) #C(1.0d0 1.0d0))

* (lla:axpy! #C(0.5 0.5) a b)
#(#C(1.5d0 1.5d0) #C(2.0d0 2.0d0) #C(2.5d0 2.5d0))

Dot product

The dot product of two vectors:

* (defparameter x #(1 2 3))
* (defparameter y #(2 3 4))
* (lla:dot x y)


The reduce function operates on sequences, including vectors (1D arrays), but not on multidimensional arrays. To get around this, multidimensional arrays can be displaced to create a 1D vector. Displaced arrays share storage with the original array, so this is a fast operation which does not require copying data:

* (defparameter a #2A((1 2) (3 4)))
* (reduce #'max (make-array (array-total-size a) :displaced-to a))

The array-operations package contains flatten, which returns a displaced array i.e doesn’t copy data:

* (reduce #'max (aops:flatten a))

An SBCL extension, array-storage-vector provides an efficient but not portable way to achieve the same thing:

* (reduce #'max (array-storage-vector a))

More complex reductions are sometimes needed, for example finding the maximum absolute difference between two arrays. Using the above methods we could do:

* (defparameter a #2A((1 2) (3 4)))
* (defparameter b #2A((1 3) (5 4)))
* (reduce #'max (aops:flatten
                    (lambda (a b) (abs (- a b))) a b)))

This involves allocating an array to hold the intermediate result, which for large arrays could be inefficient. Similarly to vectorize defined above, a macro which does not allocate can be defined as:

(defmacro vectorize-reduce (fn variables &body body)
  "Performs a reduction using FN over all elements in a vectorized expression
   on array VARIABLES.

   VARIABLES must be a list of symbols bound to arrays.
   Each array must have the same dimensions. These are
   checked at compile and run-time respectively.
  ;; Check that variables is a list of only symbols
  (dolist (var variables)
    (if (not (symbolp var))
        (error "~S is not a symbol" var)))

  (let ((size (gensym)) ; Total array size (same for all variables)
        (result (gensym)) ; Returned value
        (indx (gensym)))  ; Index inside loop from 0 to size

    ;; Get the size of the first variable
    `(let ((,size (array-total-size ,(first variables))))
       ;; Check that all variables have the same size
       ,@(mapcar (lambda (var) `(if (not (equal (array-dimensions ,(first variables))
                                                (array-dimensions ,var)))
                                    (error "~S and ~S have different dimensions" ',(first variables) ',var)))
              (rest variables))

       ;; Apply FN with the first two elements (or fewer if size < 2)
       (let ((,result (apply ,fn (loop for ,indx below (min ,size 2) collecting
                                      (let ,(map 'list (lambda (var) (list var `(row-major-aref ,var ,indx))) variables)
                                        (progn ,@body))))))

         ;; Loop over the remaining indices
         (loop for ,indx from 2 below ,size do
            ;; Locally redefine variables to be scalars at a given index
              (let ,(mapcar (lambda (var) (list var `(row-major-aref ,var ,indx))) variables)
                ;; User-supplied function body now evaluated for each index in turn
                (setf ,result (funcall ,fn ,result (progn ,@body)))))

[Note: This macro is available in this fork of array-operations, but not Quicklisp]

Using this macro, the maximum value in an array A (of any shape) is:

* (vectorize-reduce #'max (a) a)

The maximum absolute difference between two arrays A and B, of any shape as long as they have the same shape, is:

* (vectorize-reduce #'max (a b) (abs (- a b)))

Linear algebra

Several packages provide bindings to BLAS and LAPACK libraries, including:

A longer list of available packages is on CLiki’s linear algebra page.

In the examples below the lla package is loaded:

* (ql:quickload :lla)

To load "lla":
  Load 1 ASDF system:
; Loading "lla"

Matrix multiplication

The lla function mm performs vector-vector, matrix-vector and matrix-matrix multiplication.

Vector dot product

Note that one vector is treated as a row vector, and the other as column:

* (lla:mm #(1 2 3) #(2 3 4))

Matrix-vector product

* (lla:mm #2A((1 1 1) (2 2 2) (3 3 3))  #(2 3 4))
#(9.0d0 18.0d0 27.0d0)

which has performed the sum over j of A[i j] * x[j]

Matrix-matrix multiply

* (lla:mm #2A((1 2 3) (1 2 3) (1 2 3))  #2A((2 3 4) (2 3 4) (2 3 4)))
#2A((12.0d0 18.0d0 24.0d0) (12.0d0 18.0d0 24.0d0) (12.0d0 18.0d0 24.0d0))

which summed over j in A[i j] * B[j k]

Note that the type of the returned arrays are simple arrays, specialised to element type double-float

* (type-of (lla:mm #2A((1 0 0) (0 1 0) (0 0 1)) #(1 2 3)))

Outer product

The array-operations package contains a generalised outer product function:

* (ql:quickload :array-operations)
To load "array-operations":
  Load 1 ASDF system:
; Loading "array-operations"

* (aops:outer #'* #(1 2 3) #(2 3 4))
#2A((2 3 4) (4 6 8) (6 9 12))

which has created a new 2D array A[i j] = B[i] * C[j]. This outer function can take an arbitrary number of inputs, and inputs with multiple dimensions.

Matrix inverse

The direct inverse of a dense matrix can be calculated with invert

* (lla:invert #2A((1 0 0) (0 1 0) (0 0 1)))
#2A((1.0d0 0.0d0 -0.0d0) (0.0d0 1.0d0 -0.0d0) (0.0d0 0.0d0 1.0d0))


* (defparameter a #2A((1 2 3) (0 2 1) (1 3 2)))
* (defparameter b (lla:invert a))
* (lla:mm a b)
#2A((1.0d0 2.220446049250313d-16 0.0d0)
    (0.0d0 1.0d0 0.0d0)
    (0.0d0 1.1102230246251565d-16 0.9999999999999998d0))

Calculating the direct inverse is generally not advisable, particularly for large matrices. Instead the LU decomposition can be calculated and used for multiple inversions.

* (defparameter a #2A((1 2 3) (0 2 1) (1 3 2)))
* (defparameter b (lla:mm a #(1 2 3)))
* (lla:solve (lla:lu a) b)
#(1.0d0 2.0d0 3.0d0)

Singular value decomposition

The svd function calculates the singular value decomposition of a given matrix, returning an object with slots for the three returned matrices:

* (defparameter a #2A((1 2 3) (0 2 1) (1 3 2)))
* (defparameter a-svd (lla:svd a))
* a-svd
   :U #2A((-0.6494608633564334d0 0.7205486773948702d0 0.24292013188045855d0)
          (-0.3744175632000917d0 -0.5810891192666799d0 0.7225973455785591d0)
          (-0.6618248071322363d0 -0.3783451320875919d0 -0.6471807210432038d0))
         :ELEMENTS #(5.593122609997059d0 1.2364443401235103d0
   :VT #2A((-0.2344460799312531d0 -0.7211054639318696d0 -0.6519524104506949d0)
           (0.2767642134809678d0 -0.6924017945853318d0 0.6663192365460215d0)
           (-0.9318994611765425d0 -0.02422116311440764d0 0.3619070730398283d0)))

The diagonal matrix (singular values) and vectors can be accessed with functions:

(lla:svd-u a-svd)
#2A((-0.6494608633564334d0 0.7205486773948702d0 0.24292013188045855d0)
    (-0.3744175632000917d0 -0.5810891192666799d0 0.7225973455785591d0)
    (-0.6618248071322363d0 -0.3783451320875919d0 -0.6471807210432038d0))

* (lla:svd-d a-svd)
   :ELEMENTS #(5.593122609997059d0 1.2364443401235103d0 0.43380279311714376d0))

* (lla:svd-vt a-svd)
#2A((-0.2344460799312531d0 -0.7211054639318696d0 -0.6519524104506949d0)
    (0.2767642134809678d0 -0.6924017945853318d0 0.6663192365460215d0)
    (-0.9318994611765425d0 -0.02422116311440764d0 0.3619070730398283d0))


The Matlisp scientific computation library provides high performance operations on arrays, including wrappers around BLAS and LAPACK functions. It can be loaded using quicklisp:

* (ql:quickload :matlisp)

The nickname for matlisp is m. To avoid typing matlisp: or m: in front of each symbol, you can define your own package which uses matlisp (See the PCL section on packages):

* (defpackage :my-new-code
     (:use :common-lisp :matlisp))

* (in-package :my-new-code)

and to use the #i infix reader (note the same name as for cmu-infix), run:

* (named-readtables:in-readtable :infix-dispatch-table)

Creating tensors

* (matlisp:zeros '(2 2))
  0.000    0.000
  0.000    0.000

Note that by default matrix storage types are double-float. To create a complex array using zeros, ones and eye, specify the type:

* (matlisp:zeros '(2 2) '((complex double-float)))
  0.000    0.000
  0.000    0.000

As well as zeros and ones there is eye which creates an identity matrix:

* (matlisp:eye '(3 3) '((complex double-float)))
  1.000    0.000    0.000
  0.000    1.000    0.000
  0.000    0.000    1.000


To generate 1D arrays there are the range and linspace functions:

* (matlisp:range 1 10)
#<|<SIMPLE-DENSE-TENSOR: (INTEGER 0 4611686018427387903)>| #(9)
 1   2   3   4   5   6   7   8   9

The range function rounds down it’s final argument to an integer:

* (matlisp:range 1 -3.5)
 1.000   0.000   -1.000  -2.000  -3.000
* (matlisp:range 1 3.3)
 1.000   2.000   3.000

Linspace is a bit more general, and the values returned include the end point.

* (matlisp:linspace 1 10)
#<|<SIMPLE-DENSE-TENSOR: (INTEGER 0 4611686018427387903)>| #(10)
 1   2   3   4   5   6   7   8   9   10
* (matlisp:linspace 0 (* 2 pi) 5)
 0.000   1.571   3.142   4.712   6.283

Currently linspace requires real inputs, and doesn’t work with complex numbers.

Random numbers

* (matlisp:random-uniform '(2 2))
  0.7287       0.9480
  2.6703E-2    0.1834
(matlisp:random-normal '(2 2))
  0.3536    -1.291
 -0.3877    -1.371

There are functions for other distributions, including random-exponential, random-beta, random-gamma and random-pareto.

Reader macros

The #d and #e reader macros provide a way to create double-float and single-float tensors:

* #d[1,2,3]
 1.000   2.000   3.000

* #d[[1,2,3],[4,5,6]]
  1.000    2.000    3.000
  4.000    5.000    6.000

Note that the comma separators are needed.

Tensors from arrays

Common lisp arrays can be converted to Matlisp tensors by copying:

* (copy #2A((1 2 3)
            (4 5 6))
        '#.(tensor 'double-float))
  1.000    2.000    3.000
  4.000    5.000    6.000

Instances of the tensor class can also be created, specifying the dimensions. The internal storage of tensor objects is a 1D array (simple-vector) in a slot store.

For example, to create a double-float type tensor:

(make-instance (tensor 'double-float)
    :dimensions  (coerce '(2) '(simple-array index-type (*)))
    :store (make-array 2 :element-type 'double-float))

Arrays from tensors

The array store can be accessed using slots:

* (defparameter vec (m:range 0 5))
* vec
#<|<SIMPLE-DENSE-TENSOR: (INTEGER 0 4611686018427387903)>| #(5)
 0   1   2   3   4
* (slot-value vec 'm:store)
#(0 1 2 3 4)

Multidimensional tensors are also stored in 1D arrays, and are stored in column-major order rather than the row-major ordering used for common lisp arrays. A displaced array will therefore be transposed.

The contents of a tensor can be copied into an array

* (let ((tens (m:ones '(2 3))))
    (m:copy tens 'array))
#2A((1.0d0 1.0d0 1.0d0) (1.0d0 1.0d0 1.0d0))

or a list:

* (m:copy (m:ones '(2 3)) 'cons)
((1.0d0 1.0d0 1.0d0) (1.0d0 1.0d0 1.0d0))

Element access

The ref function is the equivalent of aref for standard CL arrays, and is also setf-able:

* (defparameter a (matlisp:ones '(2 3)))

* (setf (ref a 1 1) 2.0)
* a
  1.000    1.000    1.000
  1.000    2.000    1.000

Element-wise operations

The matlisp-user package, loaded when matlisp is loaded, contains functions for operating element-wise on tensors.

* (matlisp-user:* 2 (ones '(2 3)))
  2.000    2.000    2.000
  2.000    2.000    2.000

This includes arithmetic operators ‘+’, ‘-‘, ‘*’, ‘/’ and ‘expt’, but also sqrt,sin,cos,tan, hyperbolic functions, and their inverses. The #i reader macro recognises many of these, and uses the matlisp-user functions:

* (let ((a (ones '(2 2)))
        (b (random-normal '(2 2))))
     #i( 2 * a + b ))
  0.9684    3.250
  1.593     1.508

* (let ((a (ones '(2 2)))
        (b (random-normal '(2 2))))
     (macroexpand-1 '#i( 2 * a + b )))

Page source: