The Common Lisp Cookbook – 多次元配列

Table of Contents

The Common Lisp Cookbook – 多次元配列

📢 🎓 ⭐ Learn Common Lisp efficiently in videos, by the Cookbook's main contributor. Learn more.

🖊️ Discover a new Common Lisp and Coalton editor for beginners: mine and a new VSCode extension for Common Lisp: OLIVE.

Common Lisp は多次元配列をネイティブにサポートしており、1 次元配列には「ベクトル」という特別な扱いがあります。配列は汎用化されており、任意の型(element-type t)の要素を含められます。また、single-floatinteger など、特定の型の要素だけを含むよう特殊化することもできます。出発点としては、Peter Seibel による Practical Common Lisp 第11章「コレクション」 がよいでしょう。

配列に対する一般的な操作の早見表は、配列とベクトル の節にあります。

配列を操作するために Quicklisp で利用できるライブラリには、次のものがあります。

このページでは組み込みの多次元配列でできることを扱いますが、制限もあります。特に次の点です。

これら 2 つの問題は、ネイティブ配列を特殊なケースとして持つ拡張配列クラスを CLOS で定義することで解決できます。この手法を採用し、Quicklisp 経由で利用できるライブラリには、次のものがあります。

比較的新しく活発に開発されているパッケージとして MAGICL があります。これは BLAS と LAPACK ライブラリのラッパーを提供します。執筆時点ではこのパッケージは Quicklisp に含まれておらず、SBCL と CCL でのみ動きます。特に複素配列に焦点を当てているようですが、それだけに限定されているわけではありません。インストールするには、たとえば Linux/Unix では Quicklisp の local-projects ディレクトリにリポジトリをクローンします。

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

依存ライブラリ(BLAS、LAPACK、Expokit)のインストール手順は GitHub のページにあります。低水準ルーチンは外部関数をラップしているため、たとえば magicl.lapack-cffi::%zgetrf のような Fortran の名前を持ちます。これらの関数の一部には高水準インターフェイスも存在します。ソースディレクトリドキュメントを参照してください。

さらに、Common Lisp 上に構築されたドメイン固有言語を、配列を使った数値計算に利用できます。執筆時点で、これらの中でも最も広く使われ、活発にサポートされているものは次のとおりです。

CLASP は、LLVM を使って Common Lisp と他の言語(特に C++)との相互運用を容易にすることを目指すプロジェクトです。このプロジェクトの主要な用途の 1 つが、数値計算と科学技術計算です。

作成

関数 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 を渡すことで配列の添字も反復できます。さらに多くの例は array-operations の generate マニュアル を参照してください。

乱数

一様分布から得た乱数を含む 3 × 3 配列を作るには、generate から Common Lisp の random 関数を呼び出します。

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

Alexandria パッケージを使って、平均 0、標準偏差 1 のガウス分布(正規分布)に従う乱数の配列を作る例です。

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

(:ALEXANDRIA)

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

これは特に効率的ではないことに注意してください。要素ごとに関数呼び出しが必要であり、gaussian-random は 2 つの乱数を返しますが、そのうち 1 つしか使われません。

より効率的な実装や、より広い範囲の確率分布を扱うために、Quicklisp で利用できるパッケージがあります。一覧は CLiki を参照してください。

要素へのアクセス

配列の個々の要素にアクセスするには、arefrow-major-aref 関数があります。

aref 関数は、配列の次元数と同じ数のインデックス引数を取ります。インデックスは 0 から始まります。要素の格納順序は C と同じ行優先ですが、Fortran とは異なります。

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

配列を反復処理するには、次のような入れ子のループを使えます。

* (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))
         ;; 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
                     ,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)))
         ;; dimensions reversed so that innermost is last:
         (destructuring-bind ,(reverse dims-rev) ,dims
           ,result)))))

これにより、2 次元配列の内容を次のように表示できます。

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

[注: このマクロは array-operations のこのフォークで利用できますが、Quicklisp には収録されていません]

行優先順でのインデックス指定

場合によっては、特に要素ごとの演算では、次元数は重要ではありません。次元数に依存しないコードを書くには、row-major-aref で一次元化した単一のインデックスを使い、配列要素へアクセスできます。配列全体の要素数は array-total-size で得られ、一次元化したインデックスは 0 から始まります。

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

別の行列乗算の実装については、下の線形代数の節を参照してください。

要素ごとの演算

同じサイズの 2 つの数値配列を掛け合わせるには、array-operations ライブラリの 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)

each は破壊的ではなく、新しい配列を作ることに注意してください。each に渡す配列はすべて同じサイズでなければならないため、(aops:each #'+ 42 *a*) は正しくありません。

式のベクトル化

上の each 関数に代わる方法として、配列のすべての要素を反復するマクロを使えます。

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

[注: このマクロの拡張版は array-operations のこのフォークで利用できますが、Quicklisp には収録されていません]

これは次のように使えます。

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

式の本体(vectorize 式の 2 番目のフォーム)の中では、シンボル *a* は単一の要素に束縛されます。つまり、組み込みの数学関数を使えるということです。

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

さらに cmu-infix と組み合わせられます。

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

BLAS を呼び出す

高速な行列操作のために、いくつかのパッケージが BLAS のラッパーを提供しています。

Quicklisp の lla パッケージは、いくつかの BLAS 関数を呼び出す機能を提供します。

配列をスケーリングする

定数倍する演算です。

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

AXPY

これは a * x + y を計算します。ここで a は定数、xy は配列です。lla:axpy! 関数は破壊的で、最後の引数(y)を変更します。

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

内積

2 つのベクトルの内積です。

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

畳み込み

reduce 関数は、ベクトル(1 次元配列)を含むシーケンスに対して動作しますが、多次元配列には動作しません。これを回避するには、多次元配列を参照する 1 次元の displaced 配列を作成できます。displaced 配列は元の配列と記憶領域を共有するため、データをコピーせずに高速に処理できます。

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

array-operations パッケージには flatten があり、displaced 配列、つまりデータをコピーせず元の配列と記憶領域を共有する配列を返します。

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

SBCL 拡張の array-storage-vector は、同じことを実現する効率的ですが移植性のない方法を提供します。

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

ときには、より複雑な畳み込みが必要です。たとえば、2 つの配列の差の絶対値の最大値を求める場合です。上の方法を使うと次のようにできます。

* (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 と同様に、割り当てを行わないマクロを次のように定義できます。

(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 には収録されていません]

このマクロを使うと、任意の形状の配列 A における最大値を次のように求められます。

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

形状が同じであれば、任意の形状の 2 つの配列 A と B について、差の絶対値の最大値を次のように求められます。

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

線形代数

いくつかのパッケージが BLAS と LAPACK ライブラリへのバインディングを提供しています。たとえば次です。

利用可能なパッケージのより詳しい一覧は CLiki の線形代数ページにあります。

下の例では LLA パッケージがロードされています。

* (ql:quickload :lla)

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

行列の乗算

lla の関数 mm は、ベクトル同士、行列とベクトル、行列同士の乗算を行います。

ベクトルの内積

一方のベクトルは行ベクトルとして、もう一方は列ベクトルとして扱われることに注意してください。

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

これは A[i j] * x[j]j について総和したものです。

行列同士の乗算

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

これは A[i j] * B[j k]j について総和したものです。

返される配列は単純配列で、要素型が double-float に特殊化されていることに注意してください。

* (type-of (lla:mm #2A((1 0 0) (0 1 0) (0 0 1)) #(1 2 3)))
(SIMPLE-ARRAY DOUBLE-FLOAT (3))

外積

array-operations パッケージには、一般化された外積関数が含まれます。

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

これは新しい 2 次元配列 A[i j] = B[i] * C[j] を作成しています。outer 関数は任意個の入力を受け取れ、多次元の入力にも対応します。

逆行列

密行列の逆行列は 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)))
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))

逆行列を直接計算することは、特に大きな行列では一般におすすめできません。代わりに LU 分解を計算し、複数回の逆行列計算に使えます。

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

特異値分解

svd 関数は、与えられた行列の特異値分解を計算し、得られる 3 つの行列を格納するスロットを持つオブジェクトを返します。

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

Matlisp は科学技術計算ライブラリです。BLAS と LAPACK 関数のラッパーを含み、配列に対する高性能な演算を提供します。Quicklisp でロードできます。

* (ql:quickload :matlisp)

matlisp のニックネームは m です。各シンボルの前に matlisp:m: と入力せずに済むよう、Matlisp を使用する独自のパッケージを定義できます(PCL のパッケージ節を参照)。

* (defpackage :my-new-code
     (:use :common-lisp :matlisp))
#<PACKAGE "MY-NEW-CODE">

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

また、#i 中置記法リーダー(cmu-infix と同じ名前であることに注意)を使うには、次を実行します。

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

テンソルの作成

* (matlisp:zeros '(2 2))
#<|<BLAS-MIXIN SIMPLE-DENSE-TENSOR: DOUBLE-FLOAT>| #(2 2)
  0.000    0.000
  0.000    0.000
>

デフォルトでは、行列の格納型は double-float です。zerosoneseye を使って複素配列を作るには、型を指定します。

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

zerosones に加えて、単位行列を作る eye があります。

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

数列

1 次元配列を生成するには rangelinspace 関数があります。

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

range 関数は最後の引数を整数に切り下げます。

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

linspace はもう少し汎用的で、返される値に終点を含みます。

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

現在、linspace への入力は実数でなければならず、複素数には対応していません。

乱数

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

他の確率分布に対応する関数もあり、random-exponentialrandom-betarandom-gammarandom-pareto などが含まれます。

リーダーマクロ

#d#e リーダーマクロを使うと、double-floatsingle-float のテンソルを作成できます。

* #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
>

カンマ区切りが必要であることに注意してください。

配列からテンソルへ

Common Lisp 配列は、コピーして Matlisp のテンソルへ変換できます。

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

次元を指定して tensor クラスのインスタンスを作成することもできます。tensor オブジェクトの内部記憶領域は、スロット store 内の 1 次元配列(simple-vector)です。

たとえば、double-float 型のテンソルを作るには次のようにします。

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

テンソルから配列へ

テンソルの内部データにはスロットを使ってアクセスできます。

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

多次元テンソルも 1 次元配列に格納されますが、その順序は Common Lisp 配列で使われる行優先ではなく列優先です。したがって、displaced 配列として参照すると行列は転置されます。

テンソルの内容は配列へコピーできます。

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

要素へのアクセス

ref 関数は標準 Common Lisp 配列に対する aref に相当し、setf も可能です。

* (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 とともにロードされる matlisp-user パッケージには、テンソルの要素ごとに演算する関数が含まれます。

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

これには算術演算子 +-*/expt だけでなく、sqrtsincostan、双曲線関数、およびそれらの逆関数も含まれます。#i リーダーマクロはこれらの多くを認識し、matlisp-user の関数を使います。

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

Page source: ja/arrays.md

T
O
C