逆行列の計算 - 丸井綜研

Common Lispで逆行列の計算をするための準備をしてきました。

今日はついに逆行列の計算です。とは言っても、効率の良い計算をするわけではなく、単純な𝑂(𝑛³)での実装です。しかもsolve関数を呼ぶたびにLUP分解が実行されるのでかなり効率が悪いです。でも、まずは効率のことは無視して進めましょう。

前回までの記事で、LUP分解を使って連立方程式を解けるようになりました。連立方程式は、行列𝐴とベクトル𝒃に対して𝐴𝒙=𝒃が成立するベクトル𝒙を求めるということになります。ここで、𝐴の逆行列𝑋があったとすると、𝐴𝑋=𝐸ₙが成立します（𝐸ₙは𝑛×𝑛の単位行列）。𝐸ₙの一列目（[1; 0; 0; ...; 0]）は、𝐴と𝑋の1列目の積になります。つまり𝐴𝒙₁=𝒆₁ということなんですが、これはすべての列にあてはまるので𝐴𝒙ₖ=𝒆ₖとなります。この連立方程式を解けば、逆行列𝑋の𝑘列目である𝒙ₖが求まるというわけです。これを𝑛列すべてに適用すれば逆行列全体が求まります。

(defun inv (A)
  "Compute an inverse matrix of A.
Reference: Cormen, Leiserson, Rivest, and Stein. Introduction to Algorithms, 4ed. MIT Press. 2022."
  (assert (/= (det A) 0) nil "Input matrix is singular.")
  (let* ((n (array-dimension A 0))
         (x (make-array n))
         (Ainv (make-array `(,n ,n) :initial-element 0)))
    (loop for k from 0 below n
          do (let ((ek (make-array n :initial-element 0)))
               (setf (aref ek k) 1)
               (setf x (solve A ek)))
             (loop for j from 0 below n
                   do (setf (aref Ainv j k) (aref x j))))
    Ainv))

さっそく使ってみましょう。

CL-USER> (inv #2A((2 5) (1 3)))
#2A((3 -5) (-1 2))

Matlabだと以下のようになります。同じ結果ですね。

>> inv([2 5 ; 1 3])
ans =
     3    -5
    -1     2

さて、試してみて驚いたことに、このプログラムに手を加えずに複素数が入った行列の計算もできるようなのです。

CL-USER> (inv #2A((2.0d0 5.0d0) (1.0d0 #C(3.0d0 4.0d0))))
#2A((#C(0.5384615384615384d0 -0.3076923076923077d0)
     #C(-0.07692307692307693d0 0.6153846153846154d0))
    (#C(-0.015384615384615385d0 0.12307692307692308d0)
     #C(0.03076923076923077d0 -0.24615384615384617d0)))

>> inv([2 5 ; 1 3+4i])
ans =
  0.538461538461538 - 0.307692307692308i -0.076923076923077 + 0.615384615384615i
 -0.015384615384615 + 0.123076923076923i  0.030769230769231 - 0.246153846153846i

Matlabと同じくらいの精度。こんなことってある？　ビックリしました。

前述したように、そもそも計算オーダーは𝑂(𝑛³)ですし、1回しか呼び出さなくていいLUP分解を𝑛回も呼び出しているので無駄が多いです。でも、今回はここまでにします。