トップ > LISP > LU分解

LU分解

LISP

10月29日(日)にShibuya.lispのもくもく会がありました。残念ながら雑談の途中で中座しないといけなくなってしまいましたが、とても楽しかったです。この日のもくもく会の目標は「逆行列を計算できるようにすること」でした。実際は途中までしかできなかったのですが、そこまでをまとめておきます。

さて、逆行列を計算する方法をChatGPT 3.5に「正則行列の逆行列を計算する方法を列挙してください。方法の詳細には立ち入らないで大丈夫です。」と聞いてみたところ、以下の5つが回答されました。

  • 行列式と余因子行列を用いる方法
  • ガウス・ジョルダン法を用いる方法
  • LU分解を用いる方法
  • 行列の特異値分解（SVD）を用いる方法
  • 行列のクラメール法を用いる方法（特に2x2行列の場合）

ガウス＝ジョルダンの掃き出し法は学部1年生のときの授業に習いました。「たしかにそうなるけど、こんなにシンプルでええの？」と思った記憶があります。シンプルですが計算オーダーもそれなりに大きいです。

ChatGPTがあげた5つの方法の中で、広い範囲の行列に適用できて計算速度もそれなりに速そうなLU分解を用いる方法で逆行列計算を実装していくことにしました。とはいえ今日は逆行列はやらず、LU分解だけです。

少し検索すると、先達はあらまほしき事なり、以下のような記事を発見しました。

www.hello-statisticians.com

LU分解は、正方行列𝐴を下三角行列𝐿と上三角行列𝑈の積として分解するという操作です。Matlabにはluという関数がありますが、置換行列𝑃も入った𝑃𝐴=𝐿𝑈という形にするLUP分解になっています。上記記事では純粋に𝐴=𝐿𝑈にする方法が説明されていますので、それに従ってやっていこうと思います。

(defun lu-decomp (A)
  "Apply LU decomposition on square matrix A. A has to be 2-dimensional array of size n-by-n. Returns L and U in a list, with L having ones in the diagonal elements.
Reference: https://www.hello-statisticians.com/explain-terms-cat/matrix_factorization4.html"
  (assert (= (array-dimension A 0) (array-dimension A 1)) nil "Input has to be a square matrix.")
  (let* ((n (array-dimension A 0))
         (L (make-array `(,n ,n) :initial-element 0.0))
         (U (make-array `(,n ,n) :initial-element 0.0)))
    (loop for i from 0 below n
          do (loop for j from 0 below n
                   do (cond ((= i 0)
                             (if (= j 0)
                                 (setf (aref L i j) 1.0))
                             (setf (aref U i j) (aref A i j)))
                            ((= j 0)
                             (setf (aref L i j) (/ (aref A i j) (aref U 0 0))))
                            ((= i j)
                             (setf (aref L i j) 1)
                             (setf (aref U i j)
                                   (- (aref A i j)
                                      (loop for k from 0 below j
                                            sum (* (aref L i k) (aref U k j))))))
                            ((> i j)
                             (setf (aref L i j)
                                   (/ (- (aref A i j)
                                         (loop for k from 0 below j
                                               sum (* (aref L i k) (aref U k j))))
                                      (aref U j j))))
                            ((< i j)
                             (setf (aref U i j)
                                   (- (aref A i j)
                                      (loop for k from 0 below i
                                            sum (* (aref L i k) (aref U k j)))))))))
    (list L U)))

行列[1 0 2 ; -1 2 2 ; 1 2 0]をLU分解してみましょう。

CL-USER> (lu-decomp #2A((1 0 2) (-1 2 2) #(1 2 0)))
(#2A((1.0 0.0 0.0) (-1 1 0.0) (1 1 1)) #2A((1 0 2) (0.0 2 4) (0.0 0.0 -6)))

一列に表示されてしまうのでわかりにくいですが、以下のように分解されています。

L = [ 1  0  0
     -1  1  0
      1  1  1 ]

U = [ 1  0  2
      0  2  4
      0  0 -6 ]

LU分解を使うと、少なめの計算量で連立方程式を解くことができるようになりますし、それを経由して逆行列の計算をすることもできます。ただ、今回実装した単純なLU分解には、𝑎1,1成分が0だとゼロ除算が発生してしまうという問題があり、このままでは使いものになりません。そのために使われるのがLUP分解です。いずれLUP分解を実装してみましょうね。（参考：Cormen, Leiserson, Rivest, and Stein. Introduction to Algorithms, 4ed. MIT Press. (amzn.to/3QK1UFB) 第28章）

【2023-11-06追記】LUP分解を実装しました→LUP分解（Common Lisp） - 丸井綜研