独学Common Lisp

LU分解と一次方程式

概要

前回はベクトルの演算を行いました。特に重要なのは内積で、ベクトルの内積は行列の積に応用されます。

今回は行列の積を使い、一次方程式の解法として実用面で一般的なLU分解を学びます。

行列の積

ベクトルはデータの集合ですが、行列はベクトルの集合です。しかし、単純にデータの集合としてグループ化するだけではなく、「代数(algebra)」としての性格を有しています。

「代数」というのは「数字の代わり」という名の通り、方程式で文字(シンボル)をあたかも数であるかのように扱うことを指します。代数で使われるのは文字ではなく、データの関係性をもっとも簡潔に記述できる行列は、計算のためのツールとして「代数学」の領域に属します。

では、つるかめ算の方程式を思い出してください。

$$ \begin{eqnarray} x + y & = & 8 \\ 2x + 4y & = & 26 \end{eqnarray} $$

この方程式を行列に置き換えると、以下のようになります。

$$ \left( \begin{array}{c c} 1 & 1 \\ 2 & 4 \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{c} 8 \\ 26 \end{array} \right) $$

ベクトルの内積と同様に、行列の積は要素毎に掛けて全てを足すという手順で定義されていますが、「行・列」という言葉通り、行と列を掛けるという方法が採用されています。つまり、この行列表記をベクトルで表すとすれば以下のようになります(このコンテンツ内では行列の積にはドットを付けませんが、ベクトルの内積にはドットを付けることにします)。

$$ \begin{eqnarray} \left( \begin{array}{c} 1 \\ 1 \end{array} \right) \cdot \left( \begin{array}{c} x \\ y \end{array} \right) & = & 8 \\ \left( \begin{array}{c} 2 \\ 4 \end{array} \right) \cdot \left( \begin{array}{c} x \\ y \end{array} \right) & = & 26 \end{eqnarray} $$

もしこのベクトルをそのまま行列にしてしまうと、係数の部分が方程式の表記と違ってしまいます。(1 1)(2 4) をどうしても横に並べて書かないと混乱してしまいますから、方程式の書き方と行列の書き方をマッチさせるために、行列の積は行と列を掛けるスタイルになっています。

行列の1: 単位行列

行列を使った計算に入る前に、先に覚えておくべきことが一つあります。それは、行列における1です。

1という数はとても特殊な性質があり、何を掛けても元の数になる特別な数です。よって、1のかけ算はとても簡単です。計算する必要がありませんから。

では行列における1とは何でしょうか。行列における1は「単位行列」と呼ばれる特別な行列で、左上と右下の対角線上に1が並び、あとは全て0の行列です。

試しに単位行列をつるかめ算の係数行列に掛けてみます。

$$ \left( \begin{array}{c c} 1 & 1 \\ 2 & 4 \end{array} \right) \left( \begin{array}{c} 1 & 0 \\ 0 & 1 \end{array} \right) = \left( \begin{array}{c} 1 \cdot 1 + 1 \cdot 0 & 1 \cdot 0 + 1 \cdot 1 \\ 2 \cdot 1 + 4 \cdot 0 & 2 \cdot 0 + 4 \cdot 1 \end{array} \right) = \left( \begin{array}{c c} 1 & 1 \\ 2 & 4 \end{array} \right) $$

このように、単位行列の掛け算はとても計算が簡単で、単位行列であると分かれば計算をする必要さえなくなります。

単位行列はとても単純ですが便利な特性が色々とあるため、例えば逆行列などより複雑な行列の定義にも用いられます。

行列の分解: 三角行列

二次方程式の解を求める方法として、因数分解があります。

$$ x^2 - 5x + 6 = 0 $$

この式におけるxの解はすぐには分かりませんが、因数分解するとすぐに分かります。

$$ (x - 2) (x - 3) = 0 $$

解は$ x = 2, 3 $です。因数分解は「一つの式を二つの掛け算の式に変形する」だけですが、それだけで解を求めやすくなります。

行列の計算も同じで、一つの行列を複数の行列に分解することで計算がしやすくなる場合があります。その代表例が、行列の分解です。

行列の分解も二次方程式の因数分解と同じで、ある行列をより簡単な二つの行列の積で表します。行列の分解には様々なものがありますが、一次方程式の解を求めるために使われるのはLU分解です。

行列で一次方程式の解を求める方法にはいくつかありますが、ほとんどのsolve関数はLU分解を使って実装されているはずです。その理由は計算が早いからで、特に複数の変数パターンに対して式を解く場合に、一度LU分解をしておけばあとは高速に解くことができるためです。

LU分解の発想は単位行列にあります。単位行列のように、0が多いと掛け算が楽なのです。LU分解は、ある行列を0の多い2つの行列に分解する方法です。

LU分解のLとUはそれぞれ別の行列の頭文字であり、「下三角行列(Lower triangular matrix)」と「上三角行列(Upper triangular matrix)」のことを指します。それぞれ名前の通りで、下三角行列は下半分しかなく、残りは0で、上三角行列はその逆です。つまり、ある行列を下三角行列と上三角行列の積に分解するのがLU分解です。

では、やってみましょう。LとUは以下のような形をしています。

$$ L = \left( \begin{array}{c c} 1 & 0 \\ a & 1 \end{array} \right), U = \left( \begin{array}{c c} b & c \\ 0 & d \end{array} \right) $$

つるかめ算の係数行列をAとすると、以下のように変形できるはずです。

$$ \begin{eqnarray} A & = & \left( \begin{array}{c c} 1 & 1 \\ 2 & 4 \end{array} \right) \\ A & = & L U \\ \Leftrightarrow \left( \begin{array}{c c} 1 & 1 \\ 2 & 4 \end{array} \right) & = & \left( \begin{array}{c c} 1 & 0 \\ a & 1 \end{array} \right) \left( \begin{array}{c c} b & c \\ 0 & d \end{array} \right) \\ & = & \left( \begin{array}{c c} 1 \cdot b + 0 \cdot 0 & 1 \cdot c + 0 \cdot d \\ a \cdot b + 1 \cdot 0 & a \cdot c + 1 \cdot d \end{array} \right) \\ & = & \left( \begin{array}{c c} b & c \\ a \cdot b & a \cdot c + d \end{array} \right) \end{eqnarray} $$

LU分解は二次方程式の因数分解と同じで、ゴールとなる行列の形が決まっていますし、10の要素も多いので、それほど難しい処理にはなりません。今回であれば、bcがともに1とわかるので、adもそこから求まります。

$$ \begin{eqnarray} b & = & 1 \\ a \cdot b & = & 2 \\ \Leftrightarrow a \cdot 1 & = & 2 \\ \Leftrightarrow a & = & 2 \\ c & = & 1 \\ a \cdot c + d & = & 4 \\ \Leftrightarrow 2 \cdot 1 + d & = & 4 \\ \Leftrightarrow d & = & 2 \end{eqnarray} $$

これでLとUを求めることができました。

$$ L = \left( \begin{array}{c c} 1 & 0 \\ 2 & 1 \end{array} \right), U = \left( \begin{array}{c c} 1 & 1 \\ 0 & 2 \end{array} \right) $$

LU分解による一次方程式の解法

LU分解を行なったら、係数行列の代わりにLU行列を使ってつるかめ算を解きます。

元のつるかめ算の行列を$ A = LU $を使い、LとUで変形します。

$$ \begin{eqnarray} \left( \begin{array}{c c} 1 & 1 \\ 2 & 4 \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) & = & \left( \begin{array}{c} 8 \\ 26 \end{array} \right) \\ \Leftrightarrow \left( \begin{array}{c c} 1 & 0 \\ 2 & 1 \end{array} \right) \left( \begin{array}{c c} 1 & 1 \\ 0 & 2 \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) & = & \left( \begin{array}{c} 8 \\ 26 \end{array} \right) \end{eqnarray} $$

行列の掛け算は順序を交換することはできませんが、3つ以上の行列の掛け算は左からであればどこから計算しても構いません。そこで、分解したUと(x y)の掛け算を先に計算してさらに変形していきます。

$$ \begin{eqnarray} \left( \begin{array}{c c} 1 & 0 \\ 2 & 1 \end{array} \right) \left( \begin{array}{c} 1 \cdot x + 1 \cdot y \\ 0 \cdot x + 2 \cdot y \end{array} \right) & = & \left( \begin{array}{c} 8 \\ 26 \end{array} \right) \\ \Leftrightarrow \left( \begin{array}{c c} 1 & 0 \\ 2 & 1 \end{array} \right) \left( \begin{array}{c} x + y \\ 2 \cdot y \end{array} \right) & = & \left( \begin{array}{c} 8 \\ 26 \end{array} \right) \\ \Leftrightarrow \left( \begin{array}{c c} {1 \cdot (x + y)} + 0 \cdot (2 \cdot y) \\ {2 \cdot (x + y)} + 1 \cdot (2 \cdot y) \end{array} \right) & = & \left( \begin{array}{c} 8 \\ 26 \end{array} \right) \\ \Leftrightarrow \left( \begin{array}{c c} (x + y) \\ {2 \cdot (x + y)} + (2 \cdot y) \end{array} \right) & = & \left( \begin{array}{c} 8 \\ 26 \end{array} \right) \\ \end{eqnarray} $$

ここまで変形すると、左辺も右辺もベクトルになり、$ (x + y) $ が上にも下にもありますから、$ (x + y) = 8 $を使って先に下の$ y $を求めることができます。

$$ \begin{eqnarray} 2 \cdot 8 + 2 \cdot y & = & 26 \\ \Leftrightarrow y & = & 5 \end{eqnarray} $$

この$ y $で$ x $を求めることができます。

$$ \begin{eqnarray} x + 5 & = & 8 \\ \Leftrightarrow x & = & 3 \end{eqnarray} $$

長い道のりでしたが、鶴が3羽、亀が5匹という答えを求めることができました。

lu関数とピボット選択

このLU分解は大きな行列になっても計算できる一般的な方法ですが、弱点が二つあります。

  1. 対角線上に0があると計算できない
  2. 対角線上の数が小さいと誤差が大きい

対角線上の値を割る数として利用することによる宿命的な弱点です。今回はつるかめ算なので整数の範囲で計算できており、これで全く問題がありませんが、浮動小数点による計算の場合、誤差の問題も無視できません。

そこで、実際のlla:lu関数(及びlla:solve関数)では「ピボット選択」という処理が行われます。

今回の場合、LU分解を行う行列の対角線上の値とは左上の1と右下の4ですが、1の数を見るときに同じ「列」の下の方にある値をみて、1よりも大きいものがないかどうかを調べます。そして、もし自分自身(1)が一番大きければそのまま、もし他の値の方が大きければその行と自身の行を丸ごと入れ替えてしまいます。これは、方程式の並べ替えを行なっているだけです。このように、どの行を最終的に選択したか、という情報をその都度蓄積していくと、元の行列からどのように並べ替えたのかが分かります。

LLAにはLU分解を行う関数lla:luが用意されていますが、必ずこのピボット選択を行いますので、実際には行を並べ替えながら計算を行います。これは、最初に与えられる方程式の順序に関係なく、同一の連立方程式であれば同一のLU分解結果を生じさせ、追加で並べ替えた履歴を残しておく、という処理になります。

試しに、lla:luでLU分解をしてみます。

(lla:lu #2a((1 1) 
            (2 4)))
; =>
#<LU 

 L=#<LOWER-TRIANGULAR-MATRIX element-type DOUBLE-FLOAT
  1.00000       .
  0.50000 1.00000>

 U=#<UPPER-TRIANGULAR-MATRIX element-type DOUBLE-FLOAT
  2.00000  4.00000
        . -1.00000>

  pivot indices=#(2 2)>

LU分解の結果がこのページで説明したものと違うのは、ピボット選択が効いているからです。一番下の行を見てください。

pivot indices=#(2 2)

もしピボット選択が行われなければ、2行の行列なので#(1 2)となるはずです。つまり、「1行目で1行目を選択し、2行目はそのまま2行目を選択した」ということです。今回は#(2 2)となっていますから、「1行目で2行目を選択し、2行目は2行目をそのまま選択した」という意味です。

2行では分かりにくいので、3行3列で試してみます。

(lla:lu #2a((3 5 4) 
            (1 2 3) 
            (4 8 9)))
; =>
#<LU 

 L=#<LOWER-TRIANGULAR-MATRIX element-type DOUBLE-FLOAT
  1.00000        .       .
  0.75000  1.00000       .
  0.25000 -0.00000 1.00000>

 U=#<UPPER-TRIANGULAR-MATRIX element-type DOUBLE-FLOAT
  4.00000  8.00000  9.00000
        . -1.00000 -2.75000
        .        .  0.75000>

  pivot indices=#(3 3 3)>

1行1列は3ですが、3行1列の4の方が大きいので、1行目のピボット選択の結果は「3行目」になり、3行目と1行目を入れ替えます。入れ替えたあとはこのようになります。

((4 8 9)
 (1 2 3)
 (3 5 4))
pivot #(3 ? ?)

次は2行2列を見ますが、後ろはもう入れ替え済みなのでチェックしません。2行2列の2と3行2列の5を比較すると5の方が大きいので、2行目と3行目を入れ替えます。行列はこのように変化します。

((4 8 9)
 (3 5 4)
 (1 2 3))
pivot #(3 3 ?)

最後の行は入れ替えることができませんので、ピボット選択後の行列が最終的に確定し、ピボットインデックスも埋まります。

((4 8 9)
 (3 5 4)
 (1 2 3))
pivot #(3 3 3)

ピボットインデックスを見れば、入れ替えてきた履歴が分かります。今回であれば「1行目と3行目を入れ替え、新しい行列の2行目と3行目を入れ替え、最後はそのまま」ということになります。方程式の場合は右辺の定数部分(つるかめ算の場合の(8 26))も同じように入れ替えなければなりませんから、この情報を元に定数項を入れ替えて方程式を解きます。

なお、行列においてピボットとは左上から右下への対角線を指します。また、マイナスであるかどうかは計算の精度に関係しませんから、大きいかどうかの判定は「絶対値」で行います。

このように行列の計算は実際にはとても様々なことに注意する必要がありますから、信頼性の高いBLAS/LAPACKなどのライブラリを使用するべきです。

行列版つるかめ算

つるかめ算がsolveで解けることはすでに説明済みですが、わざわざLU分解を説明したのは最初に述べた通り繰り返し計算の際にL行列とU行列を使い回すことができるためです。

試しに3連続でつるかめ算を解いてみます。もちろん、LU分解も方程式を解く作業も、全部計算機が処理してくれます。ピボット選択による並べ替えもユーザには隠蔽されていますから、全く気にする必要はありません。

(ql:quickload :lla)

;; 係数行列の定義
(defvar tsurukame #2a((1 1) 
                      (2 4)))

;; LU行列の計算と保存(ピボットインデックスも一緒に保存)
(defvar lu-tk (lla:lu tsurukame))

;; 行列計算版つるかめ関数
(defun mat-tk (head leg)
  (lla:solve lu-tk (make-array '(2 1) 
                               :initial-contents `((,head) (,leg))))) 
;; 頭8, 足26
(mat-tk 8 26)
; => #2A((3.0D0) (5.0D0))

;; 頭38, 足126
(mat-tk 38 126)
; => #2A((13.0D0) (25.0D0))

;; 頭83, 足282
(mat-tk 83 282)
; => #2A((25.0D0) (58.0D0))

mat-tk関数でバッククォートとカンマを使っていますが、これはCommon Lispで一般的なテクニックなので使ってみました。これらは何もマクロだけで使われるわけではありません。テンプレートとして使う場面ではどこでも使うことができます。

複数個の方程式を解く場合は、必ずLU分解を独立して行うようにします。というのも、数学の問題であれば係数が毎回変わるかもしれませんが、つるかめ算のように係数行列が本質的には安定していることはよくあります。例えば最も重要な統計の一つである「産業連関表」はGDP計算の基礎資料になるもので、5年に一回しか改定されませんが(サボっている訳ではなく、5年かけて作っています)、産業間の生産と販売が全て集約された膨大な行列の形式をしています。ニュースでよく聞く「経済効果」というのは産業連関表の「投入係数行列」(を元にしたレオンチェフの逆行列)から求めるのですが、この行列自体は5年に一度しか変更せず、つるかめ算の頭と足に当たる「最終需要」のようなベクトルを変化させて経済波及効果を計算します。そのため、LU分解は最初の1回だけで十分です。

LLAのlla:solve関数は生の行列でもLU分解後のluオブジェクト(L行列+U行列+ピボットインデックス)でも第一引数に取ることができますから、適宜活用してください。

なお、次回は実際の産業連関表を触ってみることにしたいと思っています。産業連関表は方程式を解くための行列ではなく、本来的に行列の形式をしている行列の一つです。

また、このページの代数的行列計算には一部Maximaを使用しました。単純な2行2列の行列でも文章を書きながら計算していると分からなくなってきます。大事なのは計算方法、解法、アルゴリズムであり、数式処理システムのMaximaは変数を変数のまま計算してくれるので抽象的な計算には便利です。一般向けではCommon Lispで作られた最も有名なソフトウェアだと思います。ぜひ使ってみてください(フロントエンドはwxMaximaをお勧めします)。


Copyright © 2017- satoshiweb.net All rights reserved.