行列の転置や逆行列の計算については記事を書いたものの、単純な行列の積については記録していませんでしたので、以下にメモしておきます。洗練された高速なアルゴリズムではなく、単純な𝑂(𝑛³)の実装です。
(defun matrix-multiply (A B) "Multiply two matrices (two-dimensional arrays)." (assert (and (= (array-rank A) (array-rank B) 2) (= (array-dimension A 1) (array-dimension B 0))) (A B) "Cannot multiply ~S by ~S." A B) (let* ((Arows (array-dimension A 0)) (Acols (array-dimension A 1)) (Bcols (array-dimension B 1)) (res (make-array `(,Arows ,Bcols) :initial-element 0))) (dotimes (r Arows) (dotimes (c Bcols) (dotimes (n Acols) (setf (aref res r c) (+ (aref res r c) (* (aref A r n) (aref B n c))))))) res))
実際に使ってみると以下のようになりました。
(setq foo #2A((1 2 3) (4 5 6))) (setq bar #2A((1 2) (3 4) (5 6))) (matrix-multiply foo bar) ;=> #2A((22 28) (49 64)) (matrix-multiply bar foo) ;=> #2A((9 12 15) (19 26 33) (29 40 51))
行列の積の応用例として、無向グラフを考えてみます。図のように四つのノードがエッジによって接続されているとします。接続されているノード間では移動が可能とします。
これを行列で表してみます。各行がノードを表しており、各要素は行から列に向かってエッジがつながっているかを0と1で表しています。今回は無向グラフなので対称行列になっています。たとえば、2行3列が1なのでBからCにエッジが接続されている、4行1列は0なのでDからAには接続がない、というように読みます。*1
(setq baz #2A((0 1 1 0) (1 0 1 0) (1 1 0 1) (0 0 1 0)))
さて、すごろくのようにさいころを振って出た目の歩数だけ進めることを考えます。たとえばCから2歩ぴったりで進めるのはAとBで、Dにはたどりつけません。またC→A→C、C→B→C、C→D→Cのように、Cに2歩で戻ってくる方法は3通りあります。
これを行列の乗算で見てみると以下のようになります。行列の2乗をすると、3行3列の要素が3になっています。Cから出発して2歩でCに戻る道筋が3通りあるということを表しています。
(matrix-multiply baz baz) ;;=> #2A((2 1 1 1) ;; (1 2 1 1) ;; (1 1 3 0) ;; (1 1 0 1))
3乗すると以下のようになります。1行1列が2なので、A→B→C→A、A→C→B→Aの2通りでAに戻ってこれることが分かります。2行2列、3行3列も同様に、それぞれBとCから出発して元のノードに戻ってくるのに2通りの道筋があることを示しています。ここから、対角成分の和を6で割ると無向グラフ中の三角形の個数が計算できることが分かります。
(matrix-multiply (matrix-multiply baz baz) baz) ;;=> #2A((2 3 4 1) ;; (3 2 4 1) ;; (4 4 2 3) ;; (1 1 3 0))
行列は応用範囲が広くて楽しいですね。
