创建
- 随机数
访问数组中的元素
- 行索引
- 插入语法
Element-wise 操作
线性代数
Matlisp
- 张量（tensor）的创建
- 位操作

Common Lisp 本身就支持多维数组，比如说一维数组叫做 vectors。一般的数组可以包含任意类型 (element-type t)，或者也可以指定数组只包含指定的类型，如浮点型(single-float)或整型(integer)。Practical Common Lisp Chapter 11, Collections 比较适合新手入门。

本章主要讲解的是数组和向量的一些基本操作。

以下是一些quicklisp上操作数组的包：

array-operations：generate, permute, displace, flatten, split, combine, reshape, each。该库的作者已经不再进行维护了，但有比较活跃的库在维护着。
cmu-infix：多维数组的索引
lla：线性代数库。

本章涉及的是一些内建的对数组的操作，但这也有局限性：

无法与其他语言的数组之间的交互，比如说调用 BLAS, LAPACK 或 GSL 这些库时。
无法进行数组之间的运算，比如说当 a b 至少又一个为数组时，(+ a b) 能够返回正常结果。

以上的问题都可以通过在 CLOS 中定义一个额外的数组类来解决，或者使用 quicklisp 上的一些库：

matlisp：下文会介绍。
MGL-MAT：在机器学习库 MGL 中使用，提供了 BLAS 和 CUDA 的一些绑定。
cl-ana：数据分析库
Antik：GSLL，GNU Scientific Library 库的绑定

同时还一个新出的且很活跃的库 — MAGICL，这个库包括了 BLAS 和 LAPACK 这两个库。目前这个库不支持使用 quicklisp，只在 SBCL 和 CCL 下可以使用。这个库不只是专注与复杂的数组操作。你可以在 quicklisp 中的 local-projects 目录来安装这个库：

$ cd ~/quicklisp/local-projects
$ git clone https://github.com/rigetticomputing/magicl.git

在其相对应的 github web pages 上有这个库的详细介绍。

更深入点，操作特定语言也在 Common Lisp 中实现了，该特定支持多种数组间的操作。同时，也有一些比较收欢迎的库支持这种操作：

CLASP，一个针对 Common Lisp 与其他语言（尤其是C++）之间的交互的项目，由 LLVM 实现。该项目主要应用在数学/科学计算上。

创建

使用函数 CLHS:make-array 进行创建。

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

再复杂的数组的值也可以在创建是进行初始化。

array-operations 库提供了 generate，创建数组时很方便，该函数会将迭代进行打包。

(ql:quickload :array-operations)
To load "array-operations":
  Load 1 ASDF system:
    array-operations
; Loading "array-operations"
(:ARRAY-OPERATIONS)
(aops:generate #'identity 7 :position)
#(0 1 2 3 4 5 6)

其中 array-operations 的别名是 aops。generate 函数可以使用关键词 :subscripts 来指定从数组某个下标进行迭代赋值。具体操作见其 github repository。

随机数

初始化数组时使用随机数进行填充

(aops:generate (lambda () (random 1.0)) '(3 3))
#2A((0.2837726 0.009566903 0.84352255)
    (0.22492898 0.83147097 0.2795267)
    (0.5844146 0.7568612 0.9189848))

另一个使用随机数进行初始化数组的包是：alexandria

(ql:quickload :alexandria)
To load "alexandria":
  Load 1 ASDF system:
    alexandria
; Loading "alexandria"
(:ALEXANDRIA)
(aops:generate #'alexandria:gaussian-random 4)
#(1.6430992263819566d0 -0.5448730501612253d0 -0.027017529457142562d0
  0.974829500480401d0)

访问数组中的元素

aref

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

array-dimensions, array-rank

(array-dimensions *a*)
(4)
(array-dimensions *b*)
(2 3)
(array-rank *a*)
1
(array-dimension *a* 0)
4
(array-rank *b*)
2
(array-dimension *b* 0)
2
(array-dimension *b* 1)
3
;; loop over an array nested loops
(defparameter a #2A((1 2 3) (4 5 6)))
A
(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
NIL

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.
   Example:
    (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))
         (syms-rev (reverse syms)) ; Reverse, since starting with innermost
         (dims-rev (loop for i from 0 below rank collecting (gensym))) ; innermost dimension first
         (result `(progn ,@body))) ; Start with innermost expression
    ;; 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
                     ,result)))
    ;; 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)))
         (destructuring-bind ,(reverse dims-rev) ,dims ; Dimensions reversed so that innermost is last
           ,result)))))

然后可以输出二维数组了：

(defparameter a #2A((1 2 3) (4 5 6)))
A
(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
NIL

行索引

某些场合下，尤其是进行元素相乘的时候，数组的维度不是重要。当要写一些与维度无关的代码时，可以使用 row-major-aref 来获取相对应的元素，访问是就像将整个数组展开成一个一维数组，然后进行访问，该一维数组的大小为 array-total-size

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

插入语法

cmu-infix 提供了别样的语法，让数学表达式更易读：

    (ql:quickload :cmu-infix)
    To load "cmu-infix":
      Load 1 ASDF system:
        cmu-infix
    ; Loading "cmu-infix"
    (:CMU-INFIX)
    (named-readtables:in-readtable cmu-infix:syntax)
    (("COMMON-LISP-USER" . #<NAMED-READTABLE CMU-INFIX:SYNTAX {10030158B3}>)
 ...)
    (defparameter arr (make-array '(3 2) :initial-element 1.0))
    ARR
    #i(arr[0 1] = 2.0)
    2.0
    arr
    #2A((1.0 2.0) (1.0 1.0) (1.0 1.0))

矩阵相乘可以这样实现：

    (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]))))
      result)

Element-wise 操作

each

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

aops:each*

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

向量操作

  (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))
  (defparameter *a* #(1 2 3 4))
  *A*
  (vectorize (*a*) (* 2 *a*))
  #(2 4 6 8)
  (defparameter a #(1 2 3 4)
  A
  (defparameter b #(2 3 4 5))
  B
  (vectorize (a b) (* a (sin b)))
  #(0.9092974 0.28224 -2.2704074 -3.8356972)
  ;; combined with cmu-infix
  (vectorize (a b) #i(a * sin(b)))
  #(0.9092974 0.28224 -2.2704074 -3.8356972)

调用 BLAS

数组的缩放

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

AXPY

  (defparameter x #(1 2 3))
  ;; A
  (defparameter y #(2 3 4))
  ;; B
  (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)

当 y 是复数数组时：

  (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))

点乘

  (defparameter x #(1 2 3))
  (defparameter y #(1 2 3))
  (lla:dot x y)
  ;; => 20.0d0

Reductions

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

array-operations 中包括 flatten

  (reduce #'max (aops:flatten a)

array-storage-vector

  (reduce #'max (array-storage-vector a))
  ;; => 4
  (defparameter a #2A((1 2) (3 4)))
  ;; A
  (defparameter b #2A((1 3) (5 4)))
  ;; B
  (reduce #'max (aops:flatten
                  (aops:each
                    (lambda (a b) (abs (- a b))) a b)))
  ;; 2

vectorize-reduce

  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)))))
         ,result))))

上面的宏只能在这个 array-operations 中有用，而不是 quicklisp 中。

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

比较两个数组的大小时需要数组维度大小相同

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

线性代数

一些包含 BLAS 和 LAPACK 的包有：

更完整的库列表见：CLiki’s linear algebra page
加载 lla 包：

(ql:quickload :lla)
To load "lla":
  Load 1 ASDF system:
    lla
; Loading "lla"
(:LLA)

矩阵相乘

向量点乘

    (lla:mm #(1 2 3) #(2 3 4))
    ;; 20

矩阵与向量相乘（叉乘）

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

矩阵相乘

    (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))

注意，矩阵运算后的结果都是 double-float 类型。

Outer Product

    (ql:quickload :array-operations)
    To load "array-operations":
      Load 1 ASDF system:
        array-operations
    ; Loading "array-operations"
    (:ARRAY-OPERATIONS)
    (aops:outer #'* #(1 2 3) #(2 3 4))
    ;; #2A((2 3 4) (4 6 8) (6 9 12))

A[i j] = B[i] * C[j]

转置

    (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)))
    ;; A
    (defparameter b (lla:invert a))
    ;; B
    (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))
    (defparameter a #2A((1 2 3) (0 2 1) (1 3 2)))
    ;; A
    (defparameter b (lla:mm a #(1 2 3)))
    ;; B
    (lla:solve (lla:lu a) b)
    ;; #(1.0d0 2.0d0 3.0d0)

Singular value decomposition

    * (defparameter a #2A((1 2 3) (0 2 1) (1 3 2)))
    A
    * (defparameter a-svd (lla:svd a))
    A-SVD
    * a-svd
    #S(LLA:SVD
    :U #2A((-0.6494608633564334d0 0.7205486773948702d0 0.24292013188045855d0)
          (-0.3744175632000917d0 -0.5810891192666799d0 0.7225973455785591d0)
          (-0.6618248071322363d0 -0.3783451320875919d0 -0.6471807210432038d0))
    :D #S(CL-NUM-UTILS.MATRIX:DIAGONAL-MATRIX
         :ELEMENTS #(5.593122609997059d0 1.2364443401235103d0
                     0.43380279311714376d0))
    :VT #2A((-0.2344460799312531d0 -0.7211054639318696d0 -0.6519524104506949d0)
           (0.2767642134809678d0 -0.6924017945853318d0 0.6663192365460215d0)
           (-0.9318994611765425d0 -0.02422116311440764d0 0.3619070730398283d0)))
    (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)
   #S(CL-NUM-UTILS.MATRIX:DIAGONAL-MATRIX
   :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))

Matlisp

(ql:quickload :matlisp)
(defpackage :my-new-code
  (:use :common-lisp :matlisp))
#<PACKAGE "MY-NEW-CODE">
(in-package ;my-new-code)
;; use the #i infix reader(note the same name as for cmu-infix)
(named-readtables:in-readtable :infix-dispatch-table)

张量（tensor）的创建

* (matlisp:zeros '(2 2))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
    0.000    0.000
    0.000    0.000
>
* (matlisp:zeros '(2 2) '((complex double-float)))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: (COMPLEX DOUBLE-FLOAT)>| #(2 2)
    0.000    0.000
    0.000    0.000
>
* (matlisp:eye '(3 3) '((complex double-float)))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: (COMPLEX DOUBLE-FLOAT)>| #(3 3)
    1.000    0.000    0.000
    0.000    1.000    0.000
    0.000    0.000    1.000
>

Ranges

* (matlisp:range 1 10)
#<|<SIMPLE-DENSE-TENSOR: (INTEGER 0 4611686018427387903)>| #(9)
1   2   3   4   5   6   7   8   9
>
* (matlisp:range 1 -3.5)
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: SINGLE-FLOAT>| #(5)
1.000   0.000   -1.000  -2.000  -3.000
>
* (matlisp:range 1 3.3)
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: SINGLE-FLOAT>| #(3)
1.000   2.000   3.000
>
* (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)
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(5)
0.000   1.571   3.142   4.712   6.283
>

随机数

* (matlisp:random-uniform '(2 2))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
    0.7287       0.9480
    2.6703E-2    0.1834
>
(matlisp:random-normal '(2 2))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
    0.3536    -1.291
    -0.3877    -1.371
>

Reader macros

  * #d[1,2,3]
  #<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(3)
  1.000   2.000   3.000
  >
  * #d[[1,2,3],[4,5,6]]
  #<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
      1.000    2.000    3.000
      4.000    5.000    6.000
  >

Tensors from arrays

  * (copy #2A((1 2 3)
            (4 5 6))
        '#.(tensor 'double-float))
    #<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
    1.000    2.000    3.000
    4.000    5.000    6.000
    >
  (make-instance (tensor 'double-float)
    :dimensions  (coerce '(2) '(simple-array index-type (*)))
    :store (make-array 2 :element-type 'double-float))

Arrays from tensors

  * (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)
  * (let ((tens (m:ones '(2 3))))
    (m:copy tens 'array))
  #2A((1.0d0 1.0d0 1.0d0) (1.0d0 1.0d0 1.0d0))
  * (m:copy (m:ones '(2 3)) 'cons)
  ((1.0d0 1.0d0 1.0d0) (1.0d0 1.0d0 1.0d0))

访问元素

  * (defparameter a (matlisp:ones '(2 3)))
  * (setf (ref a 1 1) 2.0)
  2.0d0
  * a
  #<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
      1.000    1.000    1.000
      1.000    2.000    1.000
  >

位操作

* (matlisp-user:* 2 (ones '(2 3)))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 3)
    2.000    2.000    2.000
    2.000    2.000    2.000
>
* (let ((a (ones '(2 2)))
      (b (random-normal '(2 2))))
   #i( 2 * a + b ))
 #<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
   0.9684    3.250
   1.593     1.508
 >
  * (let ((a (ones '(2 2)))
      (b (random-normal '(2 2))))
   (macroexpand-1 '#i( 2 * a + b )))
(MATLISP-USER:+ (MATLISP-USER:* 2 A) B)

The CommonLisp Cookbook

06. Arrays

创建

随机数

访问数组中的元素

行索引

插入语法

Element-wise 操作

向量操作

调用 BLAS

数组的缩放

AXPY

点乘

Reductions

线性代数

矩阵相乘

向量点乘

矩阵与向量相乘（叉乘）

矩阵相乘

Outer Product

转置

Singular value decomposition

Matlisp

张量（tensor）的创建

Ranges

随机数

Reader macros

Tensors from arrays

Arrays from tensors

访问元素

位操作