Common Lispで車輪の再発明をし続けています。いろいろと統計に使う関数を作っていましたが、そろそろ逆行列を計算したくなってきました。逆行列と固有値は線形代数の基本ですものね。*1
逆行列の計算をする前には、行列式（determinant）という、逆行列が計算できるかを判定するための方法が欲しいところです。逆行列が存在しないのに計算し始めてゼロ除算が生じてしまう、などを避けたいからです。
理系であれば大学1年前期には習うと思いますが、分かりやすい解説が以下のサイトにありました。4×4以上の行列のときの行列式の計算方法はすっかり忘れていたので、参考になりました。いつもはMatlabのdet(A)で一発ですからね。
小行列を得る
2×2行列と3×3行列のときは公式どおりに計算すればいいですが、4×4以上の行列では余因子を使って分割計算していくことになります。余因子の計算には小行列（submatrix）が必要になります。小行列（
小行列 - Wikipedia）は「指定した行と列を取り除いた行列」のことで、たとえば[11 12 13 ; 21 22 23 ; 31 32 33]の1行目・2列目を取り除いた小行列は[21 23 ; 31 33]になります。今回は2次元行列しか考えないことにして、以下のように実装してみました。
(defun submatrix (A r c) "Calculates submatrix (the matrix excluding specified row and column)." (assert (= (array-rank A) 2)) (assert (>= (array-dimension A 0) 2)) (assert (>= (array-dimension A 1) 2)) (let ((res (make-array `(,(1- (array-dimension A 0)) ,(1- (array-dimension A 1))))) (ii 0) (jj 0)) (loop for i from 0 below (array-dimension A 0) do (setf jj 0) (when (/= i r) (loop for j from 0 below (array-dimension A 1) do (when (/= j c) (setf (aref res ii jj) (aref A i j)) (incf jj))) (incf ii))) res))
ループで行と列をそれぞれiとjとして回しています。指定した行・列であるrとcに一致しないときには元の行列Aから小行列resに値をコピーした上で、コピー先の添え字iiとjjを増加させます。一致したときにはコピー先の添え字は増加させていないので、コピー元の添え字であるiとjだけが増加する、ということになります。
これを実行すると以下のようになります。
CL-USER> (submatrix #2A((11 12 13) (21 22 23) (31 32 33)) 0 1) #2A((21 23) (31 33))
そういえば、特定の行だけを取り出すget-rowというような関数を作っていたので、その逆のexclude-rowとかを作れば(exclude-col (exclude-row A 0) 1)みたいにできましたね。そっちのほうが「小さい部品の組み合わせで大きなことをする」というLISP的な思想に合致していました。こんど作り直します。
行列式の計算
2×2以上の正方行列であることが条件なので、assert部分で確認しています。2×2のときはタスキがけの計算式に従って計算しますが、それより大きな行列のときには、小行列をつかって余因子を計算していきます。速度は求めていないので、以下のようにシンプルに書いてみました。detを再帰呼び出しするところが気持ちいいですね。
(defun det (A) "Calculates a determinant of a matrix." (assert (= (array-rank A) 2)) (assert (= (array-dimension A 0) (array-dimension A 1))) (assert (>= (array-dimension A 0) 2)) (if (= (array-dimension A 0) 2) (- (* (aref A 0 0) (aref A 1 1)) (* (aref A 0 1) (aref A 1 0))) (loop for n from 0 below (array-dimension A 0) sum (* (expt -1 (+ n)) (aref A 0 n) (det (submatrix A 0 n))))))
実際に使ってみると以下のようになります。計算結果はMatlabで計算したものと一致していました。
CL-USER> (det #2A((1 2 3) (4 5 6) (7 8 0))) 27 CL-USER> (det #2A((1 7 2 4) (1 5 2 4) (3 0 1 0) (2 1 5 -3))) -134
大きさn×nの行列の行列式を計算するときには(n-1)×(n-1)行列の余因子をn回計算しないといけないので、けっこうな計算量になります。計算オーダーはO(n!)ということなので、大きな行列になると速度やメモリが心配ですが、このていどで遊ぶぶんには十分です。LU分解を経由する方法というのもいずれ試してみたいです。
おまけ
Juliaで計算すると以下のように……。実用上は問題ないんですが、整数にしてほしい。なんか惜しい。
julia> using LinearAlgebra julia> det([1 2 3; 4 5 6; 7 8 0]) 26.999999999999993 julia> det([1 7 2 4 ; 1 5 2 4 ; 3 0 1 0 ; 2 1 5 -3]) -134.00000000000006
さて、こうやって車輪の再発明をしていると、なんとなく週末縄文人さんがやっていることみたいだなぁと思いました。Primitive Technologyの日本語版といった趣で、とても面白いので勝手に紹介します。
