独学Common Lisp

簡易行列計算ライブラリ

概要

前回はLU分解について学び、L行列とU行列を用いて連立一次方程式を解く方法を学びました。

当初は今回、産業連関分析を行おうと思っていましたが、いきなり専門性の高い分野に入る前に、Common Lispの「メタ言語」的な特徴を生かし、あらかじめ行列計算のための簡易的なライブラリを作成しておくことにしました。

ただし、前回のLU分解で見たように、行列計算は表面上はシンプルに見えても、正しく行うためには非常に複雑な側面も多く、基本的な方針としてはLLA(Lisp Linear Algebra)とCLNU(cl-num-utils)のよく使う関数を組み合わせて作成する方針にします。

あわせて、実用的なものにするためにCSVから行列を簡単に作成する機能も追加します。CSVの読み取りにcl-csvを使用しますので、あらかじめQuicklispでインストールしておいてください。

主要な関数・マクロの確認

行列計算においてどのような計算を頻繁に使用するかは使用者の状況により異なると思いますので、今回はもっともシンプルな構成を考えます。そのためには、それぞれの計算がどのライブラリのどの関数・マクロに対応するかをあらかじめ知っておかなければなりません。そこで、まず先に実装すべき機能と使えるオペレータを確認します。

行列の加算・減算については、clnuが利用できます。

行列の要素毎の乗算、及び要素毎の除算についてもclnuに関数があります。

行列の積については、LLAのmm関数がBLAS/LAPACKを使用した高速な計算を可能にしてくれます。

行列の行数、列数を確認する関数も基本的なツールとして必要です。今回はCommon Lispに組み込みの標準的な二次元配列を行列として扱いますから、行数・列数の確認はarray-dimension関数で対応できます。

行列の特定の行、あるいは列をベクトルとして取り出したいことはよくあります。これについてはいいオペレータが思い当たらないので、loopマクロを使って自分で実装する事にします。ただし、ベクトルを簡単に作成する関数がclnuにありますので、それを使いましょう。

cl-num-utils.matrix-shorthandというパッケージが出てきたので、簡単に紹介しますと、これはcl-num-utilsライブラリについているサブパッケージです。ベクトルや行列を手軽に作成するための関数やマクロが実装されています。今回は、このパッケージの中で行列の作成を行うmxマクロについても自作ライブラリに含めることにします。なぜわざわざ自作ライブラリに含めるかというと、パッケージ名が長いのと、行列の要素の型を指定する必要があるのですが、行列計算の場合は基本的に'double-floatで行いますので、あらかじめ型を指定しておくためです。

また、このmxマクロは私の感覚からは少し使い方が独特なので、行列風の書き方に対応した別の行列作成オペレータも用意します。ただし、そのオペレータは内部でmxを呼び出すことにします。

行列の作成に関連して、行列のコピーについてもカバーします。Common Lispにおいて配列は内部的に「参照」(ポインタ)が使われますので、数をコピーするようには複製できません。これについては定番のライブラリAlexandriaに配列のコピーに関する関数がありますので、それをそのまま利用します。

行列に関する基本的なオペレータは揃ってきたので、次は行列計算に固有の問題、つまり特定の行列を作成するようなオペレータを用意します。

まず、「転置行列」はclnuに関数があります。転置行列とは、行列の行と列を入れ替えた行列です。行列を転置する目的でも使いますが、ベクトルを行列の列に変換したいときなどにも利用します。

次は、行列の要素が全てゼロの行列です。スカラーにも0がありますから、行列にも0が必要です。また、ゼロ行列はゼロ行列として使うことよりも、特定の行数と列数を持つ行列を事前に準備しておき、あとでその要素を変更して使う場合に利用するケースも多いでしょう。また、単純な「正方行列」、つまり行数と列数が同じ行列を準備したい場合に手軽に利用できるように、列数の指定がなければ正方行列として作成するようにオプショナル引数を活用します。ただし、これはCommon Lisp標準のmake-arrayで簡単に作れますので、ライブラリは使用しません。

さらに、前回紹介したように、行列の1は「単位行列」と呼ばれ、特別な意味を持ちます。そこで、単位行列を作成するオペレータもライブラリに含めましょう。これは、前述のゼロ行列を作成した後、対角要素に全て1をセットするだけですから、自作で簡単に作ることができます。

特別な行列として最後に「逆行列」を計算するオペレータを用意します。逆行列とは、元の行列に対して行列の積を計算すると単位行列になるような行列のことで、スカラーにおける「逆数」に相当します。例えば、2の逆数は1/2であり、2つを掛けると1になります。実は経済効果分析では単位行列と逆行列が共に重要な役割を果たします。逆行列の計算は非常に複雑であり、LU分解をしてから計算するのが実用面では一般的ですが、幸いにもllaに逆行列を計算する関数がありますので、それをそのまま利用します。llaでも一回LU分解をしてから逆行列を計算します。

特定の行列を作成するオペレータは以上ですが、入出力についてもライブラリに含めておきましょう。

行列を表示する際、普通にprintなどをするとあまり綺麗に表示してくれません。行列ですからきちんと行と列が整った状態で表示してほしいと思うのが一般的だと思いますので、綺麗に行列を表示するオペレータも準備します。これも自分で実装すると地味に大変ですが、幸いclnuのサブパッケージに入っていますので、それを使います。

そして最後に、行列の入力をCSVファイルから行うことができるオペレータも準備します。これは様々な方法が考えられますが、今回はシンプルにcl-csvというライブラリを使ってリストとして読み込み、それを前述の行列作成オペレータで行列に変換するように実装します。ただし、cl-csvは標準では全て文字列として読み込んでしまうため、文字列から読み込みを行うCommon Lispの標準関数read-from-stringを使った工夫が必要です。

以上で行列に関するもっとも基本的な操作は網羅できます。巨人の肩の上に立つことで、以上全てをわずか50行程度の自作ライブラリに集約していきます。そうすると様々なライブラリの違いを吸収できます。

ASDFによるライブラリの準備

今回作成する行列計算ライブラリは単純な関数・オペレータの記述の集合としてファイルに保存しておいてもかまいませんが、行列計算一般に使えるものとしてASDF経由で簡単に呼び出せる仕様にしておきます。そのために、Minimum Common Lispで説明した通り、実際の定義記述ファイル以外に2つのファイルを準備し、全てをホームディレクトリの"common-lisp"ディレクトリ内のASDFプロジェクト用サブディレクトリに保存します。

今回作成するライブラリはmy-matrix-utilsと名付け、ニックネームとしてmatを用意することにします。よって、ホームディレクトリの"common-lisp"ディレクトリにmy-matrix-utilsディレクトリを作成し、まずはmy-matrix-utils.asdによりASDFのロードファイルを準備します。

(defsystem :my-matrix-utils
  :components
  ((:file "package")
   (:file "matrix" :depends-on ("package")))
  :depends-on
  (:alexandria
   :cl-num-utils
   :lla
   :cl-csv))

ASDF定義の書き方はMinimum Common Lispと同じですが、今回は依存ライブラリがあるので、それを明記しておきます。依存ライブラリはすでに説明した通り、alexandria, cl-num-utils, lla, cl-csvの4つです。これを:depends-onに記述しておけば、my-matrix-utilsを読み込む前に依存ライブラリを読み込んでくれます。

これはASDFのためのファイルなので、my-matrix-utilsパッケージ用にpackage.lispを用意します。なお、ASDファイルを読めばわかる通り、パッケージの定義がpackage.lispであり、定義本体はmatrix.lispに記述します。

package.lispは以下のように記述します。

(defpackage :my-matrix-utils
  (:nicknames :mat)
  (:use :cl)
  (:export
    #:m+
    #:m-
    #:m*
    #:m/
    #:m.
    #:nrow
    #:ncol
    #:mrow
    #:mcol
    #:mx
    #:mk
    #:mt
    #:mzero
    #:mident
    #:minv
    #:mpp
    #:mload))

今回、4つの外部ライブラリに依存するわけですが、定義を全て取り込む訳ではないので、:useには:clだけを記述します。ここに:llaなどを記述するとLLAなど外部ライブラリの定義についてコロンでアクセスする必要はなく、定義名(関数名・マクロ名)だけでアクセスすることができます。ただし、cl-num-utilsにはalexandriaと重複する定義(meanvariance)が含まれるので両方の定義を同時に:useする場合には注意が必要です。

:exportには外部からアクセス可能なオペレータの名前を列挙します。ここに挙げた定義名はその順番のまま、前節で説明した機能と対応します。実際は定義を書いてからpackage.lispに追加していくという方が多いでしょう。

ライブラリとして記述する場合、定義名にはある程度統一的な規則がある方が望ましいと言えます。今回は一文字目をMatrixのmにし、その後はCommon Lispで作られている数式処理システムMaximaにおける関数名と対応するようなものをイメージしてつけました。ただし、極力短くなるようにしています。なお、nrow, ncolにはmがついていませんが、これもMaximaの関数名にそのまま対応します。また、mppppはPretty Printをイメージして付けました。

それでは、ASDFによるライブラリ構築の準備が整いましたので、次節以降で実際の定義ファイルmatrix.lispを作成していきます。

関数・マクロの定義(1)

matrix.lispの一行目には、独自のパッケージに入ってから定義するためのin-packageを忘れてはいけません。

(in-package :mat)

これを忘れると標準で用意されるパッケージであるCOMMON-LISP-USERパッケージに紐づけられてしまいます。

さて、ここから関数・マクロの定義に入ります。まず、行列の加算・減算と要素毎の乗算・除算については、シンプルにそのまま外部ライブラリを参照します。また、行列の積もLLAを参照します。

(defun m+ (m1 m2) (clnu:e2+ m1 m2))
(defun m- (m1 m2) (clnu:e2- m1 m2))
(defun m* (m1 m2) (clnu:e2* m1 m2))
(defun m/ (m1 m2) (clnu:e2/ m1 m2))
(defun m. (m1 m2) (lla:mm m1 m2))

Maximaでは*は要素毎の乗算を、.は行列の積を示す関数ですので、今回はそれに習って行列の積をm.にしています。

次に、行列の行数と列数の取得も、array-dimensionで簡単に実装できます。

(defun nrow (m) (array-dimension m 0))
(defun ncol (m) (array-dimension m 1))

特定の行、列を取得するには、実装の方法は様々あると思いますが、今回は以下のようにしてみました。

(defun mrow (m row)
  (let ((l (loop
             for i from 0 to (1- (ncol m))
             collect (aref m row i))))
    (apply #'cl-num-utils.matrix-shorthand:vec 'double-float l)))

(defun mcol (m col)
  (let ((l (loop
             for i from 0 to (1- (nrow m))
             collect (aref m i col))))
    (apply #'cl-num-utils.matrix-shorthand:vec 'double-float l)))

先に定義したncol, nrowを使っています。また、リストからベクトルを作成するためのオペレータであるcl-num-utils.matrix-shorthand:vecは関数として実装されていますので、applyのような高階関数を使って引数として渡すことが可能です。これはcl-num-utils.matrix-shorthand:vecの使い方として以下のように定義されているため、リストを直接受け取ることができないからです。

(cl-num-utils.matrix-shorthand:vec 'double-float 1 2 3)
; => #(1.0D0 2.0D0 3.0D0)

要素をリストにして渡すのではなく、ただ列挙して可変長引数として渡す仕様になっています。よって、'(1 2 3)のようなリストはapplyで渡すのがベターです。(より効率的な処理が必要となるケースについては、mloadマクロの箇所で説明しています。)

しかし、ベクトルを作成するvecと同じように行列を作成するためのmxは関数ではなくマクロとして定義されています。マクロとして定義されている、ということは以下の2点に留意しなければなりません。

  1. マクロの引数は評価されてからマクロに渡されるのではなく、そのままの形でマクロの定義に渡されます。そのため、マクロ内部で引数をどのように処理・変形するかを知る必要があります。
  2. マクロの展開は実行時ではなくコンパイル時に行われます。マクロの処理が単純にソースコードの変形を目的としたものだけであれば問題は少ないですが、引数から何らかの情報を得て、その情報を元にソースコードを変形する場合、マクロ展開時に情報を得ることができる状態になっていなければなりません。

要するに、他人が書いたマクロを使うためには、相当な注意を要するということです。マクロのために書かれた優良なドキュメントが存在するか、マクロ定義のソースコードを読むか、何らかの方法でマクロの特性と使い方を知る必要があります。

mxマクロの定義を抜粋すると、以下のようになっています。

(defmacro mx (element-type &body rows)
  (let+ ((rows (map 'vector #'ensure-list rows))
         (nrow (length rows))
         (ncol (length (aref rows 0)))
         ((&once-only element-type)))
    `(make-array (list ,nrow ,ncol)
                 :element-type ,element-type
                 :initial-contents
                 (list
                  ,@(loop for row across rows collect
                             `(list
                               ,@(loop for element in row collect
                                          `(coerce ,element ,element-type))))))))

マクロの結果作成される返還後のソースコードの前に、定義内部の最初の4行(let+以下4行)が存在するので、注意点2に該当することがわかります。つまり、引数から何らかの情報を得て、それをもとにソースコードを生成するタイプのマクロです。

そして、その4行の中身を見ると、可変長引数rowsを受け取って、それをmapでベクタに変換してから、行と列の数に相当する長さをlengthで取得していることが分かります。ensure-listというのはCommon Lisp標準の関数ではないので、その内容はここだけでは不明です。

ここまで読むと、以下のようなmxマクロの特性を理解できます。

このことから、マクロmxは以下のように用いると推定されます。

(cl-num-utils.matrix-shorthand:mx 'double-float (1 2 3) (4 5 6))
; => #2A((1.0D0 2.0D0 3.0D0) (4.0D0 5.0D0 6.0D0))

このmxマクロを今回のライブラリで用いるには、以下のようにします。

(defmacro mx (&body lists)
  `(cl-num-utils.matrix-shorthand:mx 'double-float ,@lists))

単純に(mat:mx ... ...)という短いアクセスで型の指定も含めて展開されるようにしています。しかし、私の感覚では、行列は以下のように指定したいと思っています。

(mat:mk ((1 2 3) (4 5 6)))

これは、cl-csvでCSVファイルを読み込んだ際に、このような二次元リストの形で結果が得られるためでもあります。この場合、mxマクロに引数を渡すためには、最初に外側のリストを取り除き、二次元リストからリストの列挙に変換する必要があります。

また、リストを直接記述して引数として渡すのではなく、変数に束縛されたリストを渡したい場合もあります。マクロは引数が評価される前に引数をキャッチしますから、引数の評価、つまり変数に束縛されたリストを復元する作業はmxマクロに渡す手前で行う必要があります。

そこで、mxマクロの他にmkマクロを作成し、以下のように定義します。

(defmacro mk (list)
  (cond ((listp list) `(mx ,@list))
        ((symbolp list) `(mx ,@(symbol-value list)))))

このマクロの引数は2つのパターンを想定しています。一つは二次元リストで、その場合は,@を使って外側のリストを取り除き、リストの列挙の形に展開してから前述の通り定義したmat:mxに引数を渡します。

もう一つは二次元リストが束縛されたシンボルが渡されてくるパターンで、この場合は引数の評価をマクロ側で行う必要があります。そのため、Common Lisp標準のsymbol-valueを使って束縛された二次元リストを取り出し、,@によって外側のリストを取り除いてmat:maxに渡します。cl-csvを使ってCSVファイルを行列として読み込むことを考えると、この「外側のリストを取り除く」という動作がどうしても必要不可欠になります。

このように定義したmat:mkは以下の2つのパターンで動作します。

(mat:mk ((1 2 3) (4 5 6)))
; => #2A((1.0D0 2.0D0 3.0D0) (4.0D0 5.0D0 6.0D0))

(defvar mat '((1 2 3) (4 5 6)))
(mat:mk mat)
; => #2A((1.0D0 2.0D0 3.0D0) (4.0D0 5.0D0 6.0D0))

cl-num-utilsmxよりは行列風の記述に近いのではないかと思います。(このような行列の指定の仕方はCommon Lisp標準のmake-array:initial-contentsキーワードによる指定と同一でもあります。)

関数・マクロの定義(2)

さて、残りは特定の行列の作成、及び行列の入出力です。

元の行列と全く同じ行列を作成する「行列の複製」は、alexandriaの関数をそのまま参照して定義します。

(defun mcp (m) (alexandria:copy-array m))

転置行列も同じくcl-num-utilsの関数をそのまま参照します。

(defun mt (m) (clnu:transpose m))

ゼロ行列は、行数と列数を引数に取りますが、すでに述べたように正方行列の準備にも利用したいので、列数の指定がない場合は行数と列数が同じ正方行列としてゼロ行列を作成します。定義はmake-arrayを使って自作します。

(defun mzero (row &optional col)
  (when (null col) (setf col row))
  (make-array `(,row ,col)
              :element-type 'double-float
              :initial-element 0.0d0))

whenの行はunlessでも実装できますが、この部分で列数未指定の場合をキャッチします。あとは単純にmake-arrayするだけですが、行列計算は要素の型を'double-floatで統一したいので、今回も型の指定を入れています。

単位行列はゼロ行列を作ったあと、対角成分を全て1.0d0にしていくだけで作ることができます。

(defun mident (n)
  (let ((mat (mzero n n)))
    (dotimes (i n)
      (setf (aref mat i i) 1.0d0))
    mat))

特に難しいところはないように思います。dotimesではなくloopでも実装できますが、配列の要素へのアクセスだけであればdotimesの方がシンプルだと思います。

さらに、逆行列は自作するともっとも計算が複雑だと思いますが、すでに説明した通りllaに関数があるので、それをそのまま使います。

(defun minv (m) (lla:invert m))

あとは入出力です。行列を行列らしく表示する機能もcl-num-utilsをそのまま使いますが、こちらは出力先のストリームを指定する仕様として実装されているため、オプショナル引数を使って、無指定なら標準出力にし、指定があればそのストリームを使うようにします。

(defun mpp (m &optional (out t))
  (cl-num-utils.print-matrix:print-matrix m out))

オプショナル引数は地味に便利ですね。

最後に、cl-csvを用いて行列の入力をサポートします。すでに説明した通り、cl-csvでは標準で文字列として読み取るため、文字列から数値に変更するには読み取ったデータを一度マップする必要があります。cl-csv:read-csv関数はシンプルな使い方の他に高階関数として利用する方法が提供されており、:map-fnというキーワード引数にマッピング用の関数を引数として渡すことができます。この:map-fnは一行内のCSVデータをマップすることができるため、ここにread-from-stringを使うとcl-csv:read-csvの返り値を数値データにすることが可能です。

(defmacro mload (file)
  (let ((l (cl-csv:read-csv (pathname file)
                            :map-fn (lambda (row)
                                      (mapcar #'read-from-string row)))))
    `(mk ,l)))

このmloadマクロは汎用性を重視したものになっており、元のCSVが整数でも自動的に'double-floatに変換します。しかし、次回扱う産業連関表のように、元データがそもそも浮動小数点数である場合や、Excelなどを用いて事前に小数点を付けて'double-floatにすることができる場合は、clnumxマクロが自動的に型の変換を行う仕様になっているため、無駄な処理が発生し、非効率的な計算を行うことになります。実際、次回の190x190の行列でmloadを使ったところ、行列のロードが非常に重い処理になってしまい、LLAを用いて逆行列を計算する部分よりもはるかに長時間の処理を要するようになってしまいました。

そこで、元のCSVが浮動小数点数の形式で保存されており、*read-default-float-format*変数を使って'double-floatで読み込むことを前提として、今回のライブラリでは以下の仕様にして使うことにします。

(defun mload (file)
  (let ((l (cl-csv:read-csv (pathname file)
                            :map-fn (lambda (row)
                                      (mapcar #'read-from-string row)))))
    (make-array `(,(length l) ,(length (first l)))
                :element-type 'double-float
                :initial-contents l)))

こちらはマクロではなく関数として定義していることに注目してください。mxマクロを用いないので、マクロにする必要性はありません。関数として定義できる場合は、基本的には関数として定義した方がエラーなども分かりやすくて便利です。

Common Lisp標準のmake-array関数は:initial-contentsキーワードで二次元リストの形式から行列を作成することができますので、cl-csvの読み込みをそのままmake-arrayに渡すだけで対応できます。この方法により、私の環境では約15倍程度高速化しました。ただし、このライブラリを使う側で、以下の一行を付け加えることを忘れないようにしてください。

(setf *read-default-float-format* 'double-float)

これで小数点数を'double-floatとして読み込むので、変換の必要は無くなります。

ソースコード

以上、紹介した関数・マクロを全て掲載しておきます。clnumatrix-shorthandパッケージについては実際のソースコードではニックネームであるclnu.mxを用いています。それ以外は、ここで紹介したものと同じです。

;;; :my-matrix-utils "matrix.lisp"
(in-package :mat)

(defun m+ (m1 m2) (clnu:e2+ m1 m2))
(defun m- (m1 m2) (clnu:e2- m1 m2))
(defun m* (m1 m2) (clnu:e2* m1 m2))
(defun m/ (m1 m2) (clnu:e2/ m1 m2))
(defun m. (m1 m2) (lla:mm m1 m2))

(defun nrow (m) (array-dimension m 0))
(defun ncol (m) (array-dimension m 1))

(defun mrow (m row)
  (let ((l (loop
             for i from 0 to (1- (ncol m))
             collect (aref m row i))))
    (apply #'clnu.mx:vec 'double-float l)))

(defun mcol (m col)
  (let ((l (loop
             for i from 0 to (1- (nrow m))
             collect (aref m i col))))
    (apply #'clnu.mx:vec 'double-float l)))

(defmacro mx (&body lists)
  `(clnu.mx:mx 'double-float ,@lists))

(defmacro mk (list)
  (cond ((listp list) `(mx ,@list))
        ((symbolp list) `(mx ,@(symbol-value list)))))

(defun mcp (m) (alexandria:copy-array m))

(defun mt (m) (clnu:transpose m))

(defun mzero (row &optional col)
  (when (null col) (setf col row))
  (make-array `(,row ,col)
              :element-type 'double-float
              :initial-element 0.0d0))

(defun mident (n)
  (let ((mat (mzero n n)))
    (dotimes (i n)
      (setf (aref mat i i) 1.0d0))
    mat))

(defun minv (m) (lla:invert m))

(defun mpp (m &optional (out t))
  (cl-num-utils.print-matrix:print-matrix m out))

(defun mload (file)
  (let ((l (cl-csv:read-csv (pathname file)
                            :map-fn (lambda (row)
                                      (mapcar #'read-from-string row)))))
    (make-array `(,(length l) ,(length (first l)))
                :element-type 'double-float
                :initial-contents l)))

Common LispはMATLABのように行列計算を前提とした言語ではなく、汎用のプログラミング言語ですが、非常に柔軟な言語仕様により、問題領域に合わせた機能を追加することが容易です。多くの場合、自分が直面する問題に合わせた「ユーティリティ」を即席で作成し、そのユーティリティでプログラムを書くことが問題解決の最短経路となりえます。

その時に必要なのは十分に吟味され、検証されたライブラリを作成することではなく、「動く」ライブラリを「素早く」作成することです。特に計算問題においては早く計算を始めることが重要ですから、例外処理などはとりあえずあまり考えずに実装するべきだと思います。


Copyright © 2017- satoshiweb.net All rights reserved.