行列の乗算 (Common Lisp)

行列の転置や逆行列の計算については記事を書いたものの、単純な行列の積については記録していませんでしたので、以下にメモしておきます。洗練された高速なアルゴリズムではなく、単純な𝑂(𝑛³)の実装です。

(defun matrix-multiply (A B)
  "Multiply two matrices (two-dimensional arrays)."
  (assert (and (= (array-rank A) (array-rank B) 2)
               (= (array-dimension A 1) (array-dimension B 0)))
          (A B)
          "Cannot multiply ~S by ~S." A B)
  (let* ((Arows (array-dimension A 0))
         (Acols (array-dimension A 1))
         (Bcols (array-dimension B 1))
         (res (make-array `(,Arows ,Bcols) :initial-element 0)))
    (dotimes (r Arows)
      (dotimes (c Bcols)
        (dotimes (n Acols)
          (setf (aref res r c)
                (+ (aref res r c)
                   (* (aref A r n) (aref B n c)))))))
    res))

実際に使ってみると以下のようになりました。

(setq foo #2A((1 2 3) (4 5 6)))
(setq bar #2A((1 2) (3 4) (5 6)))
(matrix-multiply foo bar)   ;=> #2A((22 28) (49 64))
(matrix-multiply bar foo)   ;=> #2A((9 12 15) (19 26 33) (29 40 51))

行列の積の応用例として、無向グラフを考えてみます。図のように四つのノードがエッジによって接続されているとします。接続されているノード間では移動が可能とします。

これを行列で表してみます。各行がノードを表しており、各要素は行から列に向かってエッジがつながっているかを0と1で表しています。今回は無向グラフなので対称行列になっています。たとえば、2行3列が1なのでBからCにエッジが接続されている、4行1列は0なのでDからAには接続がない、というように読みます。*1

(setq baz #2A((0 1 1 0)
              (1 0 1 0)
              (1 1 0 1)
              (0 0 1 0)))

さて、すごろくのようにさいころを振って出た目の歩数だけ進めることを考えます。たとえばCから2歩ぴったりで進めるのはAとBで、Dにはたどりつけません。またC→A→C、C→B→C、C→D→Cのように、Cに2歩で戻ってくる方法は3通りあります。

これを行列の乗算で見てみると以下のようになります。行列の2乗をすると、3行3列の要素が3になっています。Cから出発して2歩でCに戻る道筋が3通りあるということを表しています。

(matrix-multiply baz baz)
;;=> #2A((2 1 1 1)
;;       (1 2 1 1)
;;       (1 1 3 0)
;;       (1 1 0 1))

3乗すると以下のようになります。1行1列が2なので、A→B→C→A、A→C→B→Aの2通りでAに戻ってこれることが分かります。2行2列、3行3列も同様に、それぞれBとCから出発して元のノードに戻ってくるのに2通りの道筋があることを示しています。ここから、対角成分の和を6で割ると無向グラフ中の三角形の個数が計算できることが分かります。

(matrix-multiply (matrix-multiply baz baz) baz)
;;=> #2A((2 3 4 1)
;;       (3 2 4 1)
;;       (4 4 2 3)
;;       (1 1 3 0))

行列は応用範囲が広くて楽しいですね。

*1:これはグラフ理論では隣接行列 (adjacency matrix) と呼ばれるグラフの表現方法です。よりスパースなグラフ（エッジが少ないグラフ）の場合には隣接リスト (adjacency list) を使うほうがよいです。

信号の包絡線を計算する (Common Lisp)

信号の包絡線 (signal envelope) は「信号の頂点にシーツをかぶせたときのような概形」を持つ波形です。

MatlabやRで書いた信号包絡線の計算プログラムをCommon Lispに移植しました。信号包絡線の計算方法は振幅の絶対値に対して、(1)低域フィルタをかける、(1')Savitzky-Golayフィルターがよい、(2)ピーク検出をしてスプラインでつなぐ、(2')いやいや秋間補間がよい、(3)ヒルベルト変換を使う、(3')結局それに低域フィルタかけるじゃろ、(4)スペクトル包絡には線形予測符号も使われてるよね、などなどいくつも提案できるのですが、ここではヒルベルト変換を用いて解析信号を計算する方法でいきます。解析信号の絶対値を計算すると信号包絡線になります。

まずは解析信号の実装。フーリエ変換をして、負の周波数成分をゼロにして、逆フーリエ変換します。

(defun hilbert (x)
  "Compute analytic signal of real signal vector X using Hilbert transform."
  (let* ((xx (napa-fft:fft (zeropad-pow2 x)))
         (yy (make-array (length xx) :initial-element #C(0 0)))
         (zz nil)
         (out (make-array (length x))))
    (setf (aref yy 0) (aref xx 0))
    (loop for n from 1 below (/ (length xx) 2)
          do (setf (aref yy n) (* 2 (aref xx n))))
    (setf (aref yy (/ (length xx) 2)) (aref xx (/ (length xx) 2)))
    (loop for n from (1+ (/ (length xx) 2)) below (length xx)
          do (setf (aref yy n) 0))
    (setq zz (napa-fft:ifft yy))
    (loop for n from 0 below (length out)
          do (setf (aref out n) (aref zz n)))
    out))

Matlabのhilbert()と比較してみましょう。

>> hilbert([1 2 -3 4])
ans =
   1.0000 + 1.0000i   2.0000 + 2.0000i  -3.0000 - 1.0000i   4.0000 - 2.0000i

CL-USER> (hilbert #(1 2 -3 4))
#(#C(1.0d0 1.0d0) #C(2.0d0 2.0d0) #C(-3.0d0 -1.0d0) #C(4.0d0 -2.0d0))

同じ結果になっていますね。ヨシ。

次に信号包絡線を計算します。

(defun envelop (x)
  "Calculate amplitude envelop of a signal vector X."
  (let* ((xh (hilbert x))
         (xe (make-array (length xh))))
    (dotimes (n (length xh))
      (setf (aref xe n) (abs (aref xh n))))
    xe))

こちらもMatlabと比較してみましょう。

>> abs(hilbert([1 2 -3 4]))
ans =
    1.4142    2.8284    3.1623    4.4721

CL-USER> (envelope #(1 2 -3 4))
#(1.4142135623730951d0 2.8284271247461903d0 3.1622776601683795d0
  4.47213595499958d0)

本物の音響信号でやってみましょう。スネアドラムを一発叩いた音の録音があるので、信号包絡線を計算してみます。wavreadは以前作ったプログラムです。（→Common Lisp用のwavread - 丸井綜研）

(ql:quickload :clgplot)
(let* ((snd (wavread "snare04-short.wav"))
       (x (first snd))
       (fs (second snd))
       (xe (envelop x))
       (xel (filt (biquad-lpf 100 (sqrt 1/2)) xe))
       (tt (list->vector (mapcar #'(lambda (v) (/ v fs)) (range (length xe))))))
  (clgp:plot x :x-seq tt)
  (clgp:plot xel :x-seq tt))

結果は以下のようになりました。左側が元の信号で右側が包絡線です。横軸は時間（秒）で縦軸は振幅（プラスマイナス1の範囲に正規化）。少しスムーズにするために包絡線にカットオフ周波数100 Hzの低域フィルタをかけてあげています。

ピンクノイズを生成するLISP関数

以前PythonやJSFX用のものも作ったのですが、ピンクノイズはいろいろなところで活躍するので、Common Lisp用にも作っておきます。

• ピンクノイズを作るPython関数 - 丸井綜研
• REAPERのJSFXでピンクノイズを生成してみた - 丸井綜研

Voss-McCartney法

DSP generation of Pink (1/f) Noiseには、(1)ホワイトノイズにフィルタをかけてピンクにする方法、(2)Voss-McCartney法、の二つがあると紹介されています。Voss-McCartney法はホワイトノイズを低域に重みをつけて重ね合わせていく方法です。

具体的には、いくつかの乱数生成器を用意して、それらで1, 2, 4, 8, 16, ...サンプルごとに新しい乱数を作らせ、すべてを加算して出力するというものです。たとえば4サンプルごとの生成をする乱数生成器は、4サンプルに一度だけ乱数を更新しますが、それ以外の3サンプルではその前に生成した乱数を返します。より低い頻度（ここでは1/4の標本化周波数）での乱数生成をしていることになりますね。つまり、1, 2, 4, 8, 16, ...サンプルごと生成して加算するということは、オクターブ低いホワイトノイズを次々と作って足し合わせてピンクノイズを作っていくことになります。2ⁿが大きくなればなるほど低い周波数まで伸びるピンクノイズになります。

Common Lispでの実装

とりあえず標本化周波数を決めておきます。

(defparameter +samprate+ 48000)

Voss-McCartney法の本体です。指定した秒数ぶんのピンクノイズを生成します。返り値ですが、多値を使ってピンクノイズ振幅値の列と各サンプルの時刻の列の二つを返しています。どちらも1次元配列（simple-vector）です。

(defun synth-pinknoise (dur &optional (fs +samprate+))
  "Generates pinknoise (wideband noise having -3dB/oct slope) of specified duration (seconds). Uses Voss-McCartney algorithm. http://www.firstpr.com.au/dsp/pink-noise/

Example:
(wavwrite \"pinknoise.wav\" (synth-pinknoise 2.0))"
  (let* ((lowfreq 20)
         (levels (ceiling (log (/ fs lowfreq) 2)))
         (tt (sample-time dur fs))
         (out (make-array (length tt)))
         (x (make-randn-array levels)))
    (dotimes (n (length tt))
      (dotimes (m levels)
        (if (= (1- (expt 2 m)) (mod n (expt 2 (1+ m))))
            (setf (aref x m) (randn)))
        (setf (aref out n) (+ (randn) (vecsum x)))))
    (values (normalize out) tt)))

sample-timeは、指定した秒数の各サンプルの時刻を計算する関数です。横軸に時間・縦軸に振幅をとってグラフを書く時などに使います。標本化周波数が48000 Hzのときに1秒間の時刻列は、#(0 1/48000 2/48000 ... 47999/48000)というようになります。以下のように定義しています。rangeはPythonの、linspaceはMatlabの関数をそれぞれ真似しています。

(defun sample-time (dur &optional (fs +samprate+))
  "Make a vector of time values of each sample in seconds for specified duration."
  (coerce (mapcar #'(lambda (x) (/ x fs))
                  (range (round (* dur fs)))) 'vector))

(defun range (n)
  "Generates a list of integers from 0 to N-1."
  (linspace 0 (1- n) n))

(defun linspace (x1 x2 &optional (n 100))
  "Generates a list of N (100 by default) equally spaced numbers between X1 and X2."
  (let ((x (loop for x from 0 to (1- n)
                 collect (/ x (1- n)))))
    (mapcar #'(lambda (v) (+ v x1))
            (mapcar #'(lambda (v) (* v (- x2 x1)))
                    x))))

randnとmake-randn-arrayは標準正規分布に従う乱数を作る関数で、以下のエントリに書いたものです。

marui.hatenablog.com

vecsumはベクトル内の数の総和を計算します。

(defun vecsum (vec)
  (loop as x across vec
        sum x))

ピンクノイズは音波形として使用したいので、出力時には±1の範囲に収めたいです。これをノーマライズといい、以下のnormalizeで実行しています。±1ビチビチに入れてしまうと（0 dBFSにすると）、オーディオ編集ソフトによっては「レベルオーバーしてるよ」という表示になることがあるので-1 dBFSをデフォルトにしていますが、-3 dBFSとかのほうが安全かもしれません。

(defun normalize (x &optional (dB -1.0))
  "Makes a sampled signal normalized to specified level (in dB)."
  (let ((res (make-array (array-dimensions x)
                         :initial-element 0.0))
        (gain (db->amp dB))
        (tmpmax (loop for n from 0 to (1- (length x))
                      maximize (abs (aref x n)))))
    (dotimes (n (length x))
      (setf (aref res n) (* gain (/ (aref x n) tmpmax))))
    res))

(defun db->amp (dB &optional (p0 1.0))
  (* (expt 10 (/ dB 20)) p0))

できあがったピンクノイズをオーディオ編集ソフトに取り込んで30分にまとめたのが以下の動画です。YouTubeにアップロードするときに聴覚符号化で圧縮されてしまうので、本来の音とは少し違いますが……。

youtu.be

重回帰分析 (Lisp Advent Calendar 2023)

Lisp Advent Calendar 2023の第1日目の記事です。Lispについては全くのヘナチョコなので1日目は恐縮なのですが、カレンダーの先頭が空になっているのも幸先が悪いので……。

さて、車輪の再発明を続けて、ようやく逆行列の計算までたどりつきました。今回は逆行列を使った応用をしたいと思います。

重回帰分析

逆行列の応用例として、重回帰分析を取り上げます。

重回帰分析は、ひとつの量的な応答変数を、複数の量的な説明変数で予測・説明するための線形モデルです。例えば、アパートの家賃（応答変数、単位：円）を、駅からの徒歩時間（説明変数、単位：分）と延べ床面積（説明変数、単位：m²）から推測するようなときに用いることができます。数式は 𝑦 = 𝛽₀ + 𝛽₁ 𝑥₁ + 𝛽₂ 𝑥₂ + ... + 𝛽ₙ 𝑥ₙ + 𝜖 のように書かれます。家賃の例の場合、家賃 = 係数0 + 係数1×徒歩時間 + 係数2×面積 + 誤差のようになるでしょうか。係数0は賃貸物件としてのベースラインの金額、係数1は徒歩時間が1分増えるごとに増える（あるいは減る）金額、係数2は延べ床面積が1 m²増えるごとに上がる金額を表します。また、誤差には時間と面積だけでは説明しきれない要素（築年数、階数、インターネット回線の品質など）や物件ごとの個体差が含まれます。多くのデータをもとに係数の値を求めるのが重回帰分析です。係数0と係数1だけを求めればいい場合（説明変数がひとつしかないとき）を単回帰分析と区別することもありますが、やっていることに違いはありません。

係数が決まったとして、そこから計算できる予測値は実際の値からは多少ずれます。そのずれの総量を最小化すれば最適な回帰式が得られるはずです。その発想で係数を求めるのが最小二乗法です。イメージとしては、下のグラフのように各データ点から近似直線までの距離（つまり誤差）の二乗和を最小にするように近似直線を定めるかんじです。

誤差の二乗和を最小にするためには微分を使って計算を進めることになります。やたらとめんどくさそうですが、実は簡単な行列計算に帰着します。英語版WikipediaのLinear Regressionの項に詳しく書かれていますが、まず説明変数𝑥₁ ... 𝑥ₙを以下のような行列にまとめます。ここで𝑛は説明変数の数、𝑚はデータの数となっています。行列の1列目に並んでいる1は𝛽₀と掛け算されるために追加されています。このようにすると 𝑦 = 𝛽₀1 + 𝛽₁ 𝑥₁ + 𝛽₂ 𝑥₂ + ... + 𝛽ₙ 𝑥ₙ + 𝜖 を 𝑦 = 𝑋𝛽 + 𝜖 と書け、𝑋𝛽 の部分は行列とベクトルの積で一発計算ができるようになるわけです。

ここで予測誤差である 𝜖 = 𝑦 - 𝑋𝛽 の総和を最小にしたいです。𝜖 の値は正にも負にもなるので、単純に総和を計算するとゼロになってしまいます。そこで、𝜖 の二乗和を使います。つまり、𝑋𝛽 で予測した値と、データで与えられている応答変数 𝑦 との差の二乗和ということです。以下の式を微分して最小になる 𝛽 を探ります。

ここでは導出は省略しますが、実は以下のような行列の転置、逆行列、行列の積という操作だけで 𝛽 の最適値が求まってしまいます。導出については前述のWikipedia記事を参照してください。

Lisp実装

この式を単純にCommon Lispコードにしたのが以下です。

(defun mra (y X)
  (let ((XX (make-array `(,(array-dimension X 0) ,(1+ (array-dimension X 1))) :initial-element 1)))
    (loop for i from 0 below (array-dimension X 0)
          do (loop for j from 0 below (array-dimension X 1)
                   do (setf (aref XX i (1+ j)) (aref X i j))))
    (matrix-multiply (matrix-multiply (inv (matrix-multiply (matrix-transpose XX) XX)) (matrix-transpose XX)) y)))

なんの工夫もない、単純なコードですね。まず、第1列に1が並び、それ以降にXの内容が入っているXX行列を作り、それに上記式を適用した結果を返すだけです。matrix-multiplyとmatrix-transposeはこのブログでは紹介したことはありませんが、二次元配列で与えられた行列についての積と転置の単純な実装になっています。

実行例

南風原朝和『心理統計学の基礎 https://amzn.to/48K5fKY』（有斐閣アルマ，2002年）の第8章に掲載されている例を使って実行してみましょう。この本はt検定、分散分析、回帰分析、因子分析などの背景について、数式を使いつつも直感的に理解できるように書かれている良書です。

p.226に掲載されている「母親の協調性価値得点 (𝑥₁)」「通年年数 (𝑥₂)」「子どもの協調性得点 (𝑦)」を抜き出しました。まずは、Rで実行してみます。重回帰分析はモデル式をlmで評価させるだけです。

kachi <- c(12, 12, 7, 17, 14, 9, 10, 13, 15, 12, 12, 15, 11, 14, 17, 17, 16, 15, 15, 10, 12, 9, 12, 12, 19, 11, 14, 15, 15, 15, 16, 15, 12, 10, 11, 12, 15, 13, 15, 12, 12, 12, 13, 17, 13, 11, 14, 16, 12, 12)
nensu <- c(2, 2, 2, 3, 2, 2, 3, 3, 3, 1, 3, 3, 2, 2, 4, 2, 4, 3, 4, 2, 2, 1, 2, 2, 4, 2, 3, 2, 3, 3, 2, 3, 2, 2, 3, 1, 2, 3, 2, 2, 2, 3, 3, 3, 2, 3, 2, 4, 2, 2)
kyouchou <- c(6, 11, 11, 13, 13, 10, 10, 15, 11, 11, 16, 14, 10, 13, 12, 15, 16, 14, 14, 8, 13, 12, 12, 11, 16, 9, 12, 13, 13, 14, 12, 15, 8, 12, 11, 6, 12, 15, 9, 13, 9, 11, 14, 12, 13, 9, 11, 14, 16, 8)
res <- lm(kyouchou ~ kachi + nensu)
summary(res)

以下のように、子どもの協調性得点 = 5.1716 + 0.3027×母親の協調性価値得点 + 1.1259×通園年数という結果になりました。

Call:
lm(formula = kyouchou ~ kachi + nensu)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.0563 -1.4780 -0.0563  1.4751  4.9437 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   5.1716     1.6647   3.107   0.0032 **
kachi         0.3027     0.1467   2.064   0.0446 * 
nensu         1.1259     0.4713   2.389   0.0210 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.12 on 47 degrees of freedom
Multiple R-squared:  0.3137,    Adjusted R-squared:  0.2845 
F-statistic: 10.74 on 2 and 47 DF,  p-value: 0.000144

今回のLisp実装でこの3つの係数が得られれば成功ということにします。

(defparameter *kachi* #(12 12 7 17 14 9 10 13 15 12 12 15 11 14 17 17 16 15 15 10 12 9 12 12 19 11 14 15 15 15 16 15 12 10 11 12 15 13 15 12 12 12 13 17 13 11 14 16 12 12))
(defparameter *nensu* #(2 2 2 3 2 2 3 3 3 1 3 3 2 2 4 2 4 3 4 2 2 1 2 2 4 2 3 2 3 3 2 3 2 2 3 1 2 3 2 2 2 3 3 3 2 3 2 4 2 2))
(defparameter *kyouchou* #(6 11 11 13 13 10 10 15 11 11 16 14 10 13 12 15 16 14 14 8 13 12 12 11 16 9 12 13 13 14 12 15 8 12 11 6 12 15 9 13 9 11 14 12 13 9 11 14 16 8))

(mra (matrix-transpose *kyouchou*) (cbind (matrix-transpose *kachi*) (matrix-transpose *nensu*)))

matrix-transposeは行列の転置ですが、入力がsimple vectorの場合はn×1のarrayにします。cbindは二つの2次元arrayを横方向につなげる関数です。この結果、#2A((64133/12401) (11263/37203) (13962/12401))が得られました。64133÷12401 = 5.171599、11263÷37203 = 0.3027444、13962÷12401 = 1.1258769なので、Rでの計算結果と同じになりました。

ここでは決定係数（R²）の計算はしていませんが、𝑋𝛽の分散を𝑦の分散で割ると得られます。今回のデータでは0.31367722です。

逆行列の計算

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

• 行列式の計算（Common Lisp） - 丸井綜研
• LU分解（Common Lisp） - 丸井綜研
• LUP分解（Common Lisp） - 丸井綜研
• 連立方程式を解く（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分解を𝑛回も呼び出しているので無駄が多いです。でも、今回はここまでにします。

連立方程式を解く

先日以下の記事にしたように、LUP分解をすると正則行列𝐴に対して𝑃𝐴 = 𝐿𝑈という分解ができるのでした。𝐿は下三角行列、𝑈は上三角行列、𝑃は並べ替えの情報が入っています。

marui.hatenablog.com

これを使って連立方程式𝐴𝑥 = 𝑏を解くプログラムを作ります。LU分解のときにお世話になったサイトをふたたび参考にします。

www.hello-statisticians.com

このページの「𝑈𝑥 = 𝑦を𝑥に関して解く」の最後の式（𝑖 < 𝑛の時の𝑥_𝑖の計算）の分母にある𝑢_𝑖,𝑛は𝑢_𝑖,𝑖じゃないといけないので、下記コードではそのようにしています。（コメントが受け付けられていないので、どうやってページの著者に連絡すればいいのやら……）

また、上記解説記事とは違い、今回は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追記：続いて逆行列を計算しました　→逆行列の計算 - 丸井綜研）

LUP分解

【2023-11-09：入力行列を書き換えながら処理をするようになっていたので、呼び出し元でも行列の内容が書き変わってしまっていました。複製を作って作業するように修正しました。】

昨日LU分解のコードを掲載しましたが、単純なアルゴリズムなので、𝑎₁,₁要素がゼロだとうまく動かないという問題がありました。

marui.hatenablog.com

それを解消するのが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） - 丸井綜研）