NumPyはどうやって逆行列を計算しているのか

はじめに

こんにちは！
GA technologiesのAdvanced Innovation Strategy Center（AISC）に所属している三田です。

現在、AISCでは有志で集まって数学の勉強会を開いております。現在のところ線形代数、微分積分について課題図書を決めて毎週輪読していき、疑問点を議論したり問題を解いたりしています。私は数学がとても苦手なので学びの多い会になっております。

最近はその勉強会の参加者で「なにか数学に関する記事を書こう！」となっており、今のところ以下の記事がでております。よろしければこちらもご覧いただければ幸いです。

本記事は線形代数の勉強をする中で抱いた疑問について調べてみた結果を共有するものになります。

線形代数で出てくる概念の計算について

線形代数ではいろいろな概念だったり、それらを求めるアルゴリズムが出てきますが、計算量の多いものもあります。

例えば、$${n}$$次の正方行列$${A}$$があるとして、その逆行列$${ A^{-1} }$$を計算したいとします。計算する方法として、行列式$${ |A| }$$と余因子行列$${ \tilde{A} }$$を使って

$$
A^{-1} = \frac{1}{|A|} \tilde{A}
$$

とする計算方法が線形代数の入門書に書かれていたりしますが、これを愚直に計算すると行列式と余因子行列を計算することになります。行列式を

$$
|A| = \sum_{\sigma \in S_n} \operatorname{sgn}(\sigma)
\prod_{i=1}^n a_{i \sigma(i)}
$$

という式に従って愚直に計算すると、置換$${\sigma}$$の集合$${S_n}$$は$${n!}$$個あるため$${ \sum_{\sigma \in S_n}  }$$は$${ n! }$$回の足し算を意味しており、$${ \prod_{i=1}^n }$$もあるため、計算量が非常に多いことがわかります。余因子行列の計算では行列式をたくさん計算することになるので、上の式のように逆行列を求めると膨大な計算量が必要になりそうです。

…といったことを考えていたとき、「行列計算をするライブラリは愚直にこの計算をしているのだろうか…？」という疑問が湧きました。

機械学習関連の論文を読んでいても、計算を高速化するために行列計算部分を工夫して計算量を減らす話が出てきたりすることがちらほらあります。きっと線形代数のライブラリたちはもっと効率的な計算をしているはず…！

そこで、試しにnumpy（+ついでにscipy）の逆行列の計算の中身を覗いてみることにしました。

numpyの逆行列の実装について

numpyで逆行列を求める関数は numpy.linalg.inv になります。以下のコードのように使われるものになります。

>>> import numpy as np
>>> A = np.array([
...   [1, 2],
...   [3, 4]
... ])
>>> np.linalg.inv(A)
array([[-2. ,  1. ],
       [ 1.5, -0.5]])

numpy.linalg.invのソースコードをたどってみると、LAPACKというパッケージのgesvという関数を呼び出して使っているようです。

LAPACK （Linear Algebra PACKage）は線形代数で使用する基本的な計算を実装したオープンソースのパッケージです。1992年に公開された、長い歴史を持つパッケージになります。

さて、numpyで使われているgesvは何なのかといいますと、LU分解を利用して連立1次方程式を解くサブルーチンになります（LAPACK: dgesv）。

LU分解とは、正方行列$${ A }$$を単位下三角行列$${ L }$$と上三角行列$${ U }$$の積の形

$$
A = LU
=
\left[\begin{array}{cccc}
1 & & & \\
l_{21} & 1 & & \\
\vdots & \ddots & \ddots & \\
l_{n 1} & \cdots & l_{n, n-1} & 1
\end{array}\right]\left[\begin{array}{cccc}
u_{11} & u_{12} & \cdots & u_{1 n} \\
& u_{22} & \cdots & u_{2 n} \\
& & \ddots & \vdots \\
& & & u_{n n}
\end{array}\right]
$$

に分解する手法です（※置換行列$${P}$$も入れて$${ A = PLU }$$と分解するものもあります）。詳しくは後述しますが計算量削減のために使われる手法になります。

また、なぜ逆行列の計算に連立1次方程式が出てくるかというと、連立1次方程式

$$
A X = B
$$

において $${ B }$$ に単位行列をおけば $${ X }$$ を求めることで $${ X = A^{-1} }$$が計算できることを利用して解いているのだと思われます。

LU分解をしておくと、連立1次方程式$${ A\boldsymbol{x} = \boldsymbol{b} }$$は、$${ L \boldsymbol{y} = \boldsymbol{b} }$$と$${ U \boldsymbol{x} = \boldsymbol{y} }$$を使って解くことができます。1回分解したLとUを使いまわして複数の$${ \boldsymbol{x} }$$、$${ \boldsymbol{b} }$$についての解を得ることができるため、上述の $${ X }$$ を効率的に求められます。

なお、scipyの逆行列もソースコードを覗いたところLAPACKが使われている様子で、LU分解を行うgetrfと、LU分解した結果から逆行列を計算するgetriが使われているようでした。

（ちなみに先ほどからgesvやgetrfといったわかりにくい関数名ばかり出てきますが、こういう名前になった理由は昔のFORTRANの規格でルーチン名の上限が6文字だったかららしいです（幸谷, 2016）。当時のプログラマは大変だったでしょうね…）

LU分解を使う場合の計算量

調べてみたところ、LU分解を使って逆行列を求める方法は計算量の観点でかなり効率的なようです。

LU分解に約$${ n^3/3 }$$回の計算が必要で、分解後の逆行列の計算に約$${ 2n^3/3 }$$回の計算をするため、トータルで$${ n^3 }$$程度の計算回数になるようです（伊理 & 藤野, 1985）。

（※線形代数の入門書に出てくることも多いGauss-Jordanの消去法も $${ O(n^3) }$$ ではあるのですが、スパースな（ゼロの多い）行列だった場合にLU分解は分解後のLとUの行列にもスパースな特性が残るためLUのほうが計算量を削減できる利点があり、LU分解が好まれるようです）

冒頭に述べた

$$
A^{-1} = \frac{1}{|A|} \tilde{A}
$$

で愚直に計算する方法と比べてみればかなり効率的に計算されていますね。

numpyの行列式の実装について

ところで、冒頭で行列式を愚直に計算したら計算量が膨大になりそうだと述べました。行列式はnumpyやscipyではどのように計算されているのでしょうか。こちらも見てみます。

numpy.linalg.det も scipy.linalg.det も、ドキュメントを読むとLAPACKのgetrfでLU分解してから計算している旨が書かれています。

$${ |A| }$$ をLU分解して整理すると

$$
\begin{aligned}
|A| &= |LU|\\
&= |L| |U| \quad (\because 行列式の性質のため)\\
&= 1 \times |U| \quad (\because 三角行列の行列式は対角成分の積と等しいため)\\
&= |U|
\end{aligned}
$$

となります。三角行列の行列式は対角成分の積と等しいため、$${ U }$$ の対角成分の積をとれば$${ |A| }$$ が求まることになります。LU分解さえしてしまえば、非常に簡単に行列式が推定できるのですね。

所感

LU分解について、私は存在自体は一応認知していたものの予想以上にさまざまな計算処理の中身にLU分解が使われているのだなと驚きました。また、実際に計算の効率が飛躍的に上がっているので本当にすごい技術だと思います。

ところで、機械学習分野で出てくる連立1次方程式や逆行列の計算の高速化の話ではコレスキー分解や共役勾配法といった単語を見かけた覚えがあるのですが、そういうものはnumpyでは使われていませんでした。

調べたところ、それらの手法は計算量の大幅な削減が可能なものの適用可能な行列が対称行列に限られるなどの制約と引き換えの高速化のようです（杉原 & 室田, 2016）。例えば線形回帰の正規方程式に出てくる $${ X^T X }$$ のような行列であればそういった方法が利用可能なのだろうとは思いますが、numpyのようなライブラリでは一般の正方行列に対して利用可能でなければならないためLU分解を使っているのかなと思いました。
（※一般化共役残差法など非対称でも使えるように拡張した共役勾配法系の手法もあるようですが、制約がゼロではないようですし、共役勾配法などの反復法は結局近似解なのでnumpyでは避けているのかな…？と思っています）

逆に言えば、行列の特性に合わせた計算方法を使うことで、 高速に計算できるようになる可能性がありそうです。今後、大規模な行列計算をすることがあれば気をつけておきたいですし、そのときの自分の引き出しを増やすためにも引き続き数学の勉強を進めていきたいと思います。

参考文献

伊理正夫, & 藤野和建. (1985). 数値計算の常識. 共立出版.
幸谷智紀. (2016). LAPACK/BLAS入門. 森北出版.
杉原正顯, & 室田一雄. (2016). 線形計算の数理. 岩波書店.