先日以下の記事にしたように、LUP分解をすると正則行列𝐴に対して𝑃𝐴 = 𝐿𝑈という分解ができるのでした。𝐿は下三角行列、𝑈は上三角行列、𝑃は並べ替えの情報が入っています。
これを使って連立方程式𝐴𝑥 = 𝑏を解くプログラムを作ります。LU分解のときにお世話になったサイトをふたたび参考にします。
このページの「𝑈𝑥 = 𝑦を𝑥に関して解く」の最後の式（𝑖 < 𝑛の時の𝑥𝑖の計算）の分母にある𝑢𝑖,𝑛は𝑢𝑖,𝑖じゃないといけないので、下記コードではそのようにしています。（コメントが受け付けられていないので、どうやってページの著者に連絡すればいいのやら……）
また、上記解説記事とは違い、今回はLU分解ではなくLUP分解を経由しているので、𝐴𝑥 = 𝑏 → 𝑃𝐴𝑥 = 𝑃𝑏 → 𝐿𝑈𝑥 = 𝑃𝑏 という式変形をしてから計算を開始することになります。そのため、最初に𝑃𝑏の計算（𝑃の値に従って𝑏を並べ替え）が必要になります。
(defun solve (A b) "Solve a system of linear equations of Ax=b for x where A is n-by-n matrix and b is a vector of length n. Reference: https://www.hello-statisticians.com/explain-terms-cat/matrix_factorization5.html" (assert (/= (det A) 0) nil "Input matrix is singular.") (let* ((n (array-dimension A 0)) (x (make-array n)) (y (make-array n)) (bb (make-array n)) (LUP (lup-decomp A)) (L (car LUP)) (U (cadr LUP)) (P (caddr LUP))) (loop for i from 0 below n ;reorder b to match P do (loop for j from 0 below n do (if (= (aref P i j) 1) (setf (aref bb i) (aref b j))))) (setf (aref y 0) (aref bb 0)) ;set initial value (loop for i from 1 below n ;forward substitution do (setf (aref y i) (- (aref bb i) (loop for j from 0 to i sum (* (aref L i j) (aref y j)))))) (setf (aref x (1- n)) (/ (aref y (1- n)) (aref U (1- n) (1- n)))) (loop for i from (- n 2) downto 0 ;backward substitution do (setf (aref x i) (/ (- (aref y i) (loop for j from (1- n) downto i sum (* (aref U i j) (aref x j)))) (aref U i i)))) x))
Common LispとMatlabで同じ連立方程式を解いてみましょう。
CL-USER> (solve #2A((1 2 0) (3 4 4) (5 6 3)) #(3 7 8)) #(-7/5 11/5 3/5)
>> [1 2 0 ; 3 4 4 ; 5 6 3] \ [3 ; 7 ; 8] ans = -1.4000 2.2000 0.6000
分数か少数かの違いはあれど、同じ結果です。
（2023-11-19追記：続いて逆行列を計算しました →逆行列の計算 - 丸井綜研）
