【2023-11-09：入力行列を書き換えながら処理をするようになっていたので、呼び出し元でも行列の内容が書き変わってしまっていました。複製を作って作業するように修正しました。】
昨日LU分解のコードを掲載しましたが、単純なアルゴリズムなので、𝑎₁,₁要素がゼロだとうまく動かないという問題がありました。
それを解消するのがLUP分解です。各列を見ていって、対角要素の絶対値が最大になるように行の入れ替えを行ってからLU分解を行うことで、ゼロ除算を防ぐだけでなく小さい数での除算が少なくなって計算結果が安定するというものです。結果として、正則行列𝐴に対して𝑃𝐴=𝐿𝑈となる𝐿, 𝑈, 𝑃という3つの行列が得られます。
今日は、Cormen, Leiserson, Rivest, and Stein. Introduction to Algorithms, 4ed. MIT Press. (amzn.to/3QK1UFB) の第28章掲載の疑似コードをそのままCommon Lispで実装しました。コードが与えられているので簡単ですが、もともと疑似コードがC言語っぽいので、できあがったLISPコードもCみたいになってしまいますね。
(defun lup-decomp (A) "Apply LUP decomposition on square matrix A such that PA=LU. A has to be 2-dimensional array of size n-by-n. Returns three arrays L, U, and P in a list. Reference: Cormen, Leiserson, Rivest, and Stein. Introduction to Algorithms, 4ed. MIT Press. 2022." (assert (/= (det A) 0) nil "Input matrix has be non-singular.") (let* ((n (array-dimension A 0)) (AA (make-array `(,n ,n))) (L (make-array `(,n ,n) :initial-element 0)) (U (make-array `(,n ,n) :initial-element 0)) (P (make-array `(,n ,n) :initial-element 0)) (perm (make-array n)) (tmp 0) (kp 0)) (loop for i from 0 below n do (loop for j from 0 below n do (setf (aref AA i j) (aref A i j)))) ;make copy of A (loop for i from 0 below n do (setf (aref perm i i)) ;initialize perm to the identity permutation (setf (aref L i i) 1)) ;set 1 for diagonal element of L (loop for k from 0 below n do (setf tmp 0) (loop for i from k below n ;find largest absolute value in column k do (if (> (abs (aref AA i k)) tmp) (setf tmp (abs (aref AA i k)) kp i))) ;row number of the largest found so far (setf tmp (aref perm k) (aref perm k) (aref perm kp) (aref perm kp) tmp) (loop for i from 0 below n do (setf tmp (aref AA k i) ;exahange rows k and kp (aref AA k i) (aref AA kp i) (aref AA kp i) tmp)) (loop for i from (1+ k) below n do (setf (aref AA i k) (/ (aref AA i k) (aref AA k k))) (loop for j from (1+ k) below n do (setf (aref AA i j) ;compute L and U in place in A (- (aref AA i j) (* (aref AA i k) (aref AA k j))))))) (loop for i from 0 below n do (setf (aref P i (aref perm i)) 1)) ;create P (loop for i from 0 below n do (loop for j from 0 below n do (if (> i j) ;extract L and U from A (setf (aref L i j) (aref AA i j)) (setf (aref U i j) (aref AA i j))))) (list L U P)))
同書の例をつかって計算したところ、掲載されていた結果と同じになりました。あたりまえですが。
CL-USER> (lup-decomp #2A((2 0 2 0.6) (3 3 4 -2) (5 5 4 2) (-1 -2 3.4 -1))) (#2A((1 0 0 0) (2/5 1 0 0) (-1/5 1/2 1 0) (3/5 0 0.39999995 1)) #2A((5 5 4 2) (0 -2 2/5 -0.19999999) (0 0 4.0000005 -0.5) (0 0 0 -3.0)) #2A((0 0 1 0) (1 0 0 0) (0 0 0 1) (0 1 0 0)))
Matlabでやってみると、次のようになりました。こちらも同じ結果ですね。
>> [L,U,P] = lu([2 0 2 0.6 ; 3 3 4 -2 ; 5 5 4 2 ; -1 -2 3.4 -1]) L = 1.0000 0 0 0 0.4000 1.0000 0 0 -0.2000 0.5000 1.0000 0 0.6000 0 0.4000 1.0000 U = 5.0000 5.0000 4.0000 2.0000 0 -2.0000 0.4000 -0.2000 0 0 4.0000 -0.5000 0 0 0 -3.0000 P = 0 0 1 0 1 0 0 0 0 0 0 1 0 1 0 0
次は連立方程式を解いてみましょうか。（→ 2023-11-08：連立方程式を解く（Common Lisp） - 丸井綜研）
