BETA

安易に逆行列を数値計算するのはやめよう

投稿日:2019-11-15
最終更新:2019-11-15

逆行列を使った計算というのは機械学習ではそれなりに出てきます。
例えば、最小二乗法では
$$ x = (X^T X) ^{-1} Xb$$
の形の式を計算する必要がありますし、正規分布の分散を扱うときにも逆行列が出てきます。
こういうときにnp.linalg.invを使って逆行列を求めて、その後にベクトルとの積を求めるは簡単にできますから、特に何も考えずにそういうふうにしたくなります。

でもそれって本当に逆行列の計算が必要ですか?

多くの問題では逆行列の値そのものよりも、$x=A^{-1}b$のような逆行列とベクトルとの積が必要になります。そのような場合、実は計算はもっと速くできますよ、というのが今日のお話です。

ただし今回は式を深く追うことはしませんので、細かい計算量などが気になる方は別途どこかの講義資料などの参照をお願いします。

逆行列を求めるための計算量

逆行列を求めるための方法として多くの人が思いつくのが、おそらく線形代数の教科書に載っている掃き出し法でしょう。掃き出し法は逆行列を求めたい行列$A$に対して操作をおこない、単位行列にしていくやり方ですね。
行列$A$のサイズを$n \times n$としたとき、掃き出し法に必要な乗除算は$n^3$回、引き算は$n(n-1)^2$回です。
また別途、行列$b$との積を計算する場合には乗算が$n^2$回、足し算が$n(n-1)$回かかることに注意してください。

実際にはnp.linalg.invはこの方法ではなく、後述する方法を利用して(半ば無理やり?)逆行列を求めますが、そうしても計算量は上記と同じ程度になります。

連立一次方程式を解く方法

$x=A^{-1}b$の計算は、$Ax=b$の形をした連立一次方程式とみなすことができます($x=A^{-1}b$の両辺に左から$A$を掛けるとわかりますね)。よって、連立一次方程式が解ければ、逆行列を求める必要はないということです。

以下ではnp.linalg.solveでもおこなわれている、LU分解と前進後退代入を使った連立一次方程式の解き方について述べます。

LU分解

行列$A$に対してLU分解をおこなうことを考えます。LU分解というのは下三角行列$L$と上三角行列$U$の積に行列$A$を分解することを指します。つまり、$$A = LU$$が成り立つような$L$と$U$を求めます。

LU分解の計算量は乗除算が$(n-1)(n^2+n+3)/3$回で引き算が$n(n-1)(2n-1)/6$回です。ここまでは先程出てきた逆行列を求めるための計算量よりも大分少ない計算量です。

もちろんLU分解だけでは連立一次方程式は解けず、次の前進後退代入をおこなう必要があります。

前進後退代入

LU分解が済んでいるとすると、$Ax=b$は$LUx=b$とあらわせます。$y=Ux$とおいてあげると、
$$Ax=LUx= Ly=b$$
となりますので、$Ly=b$の連立一次方程式が出てきます。これを$y$について解くと次に
$$Ux = y$$
の連立一次方程式があらわれます。最後にこれを$x$について解くことで、ようやく欲しかった$x$が求まります。

$Ly=b$と$Ux=y$という連立一次方程式を解くなんて計算が重そうだ!と思うかもしれません。
しかしながら、$L$は下三角行列、$U$は上三角行列であるということを考慮するとそれほど計算量は多くなりません。実際、

  • $Ly=b$を求める計算(前進代入):乗算$n(n-1)/2$回、加減算$n(n-1)/2$回
  • $Ux=y$を求める計算(後退代入):乗除算$n(n+1)/2$回、加減算$n(n-1)/2$回
  • 上2つの計算量の和:乗除算$n^2$回、加減算$n(n-1)$回

となります。なんとこれは前述した$A^{-1}$を$b$に掛けるときの計算量と等しいです!
一見大変そうな計算をしているのに、実は行列とベクトルの積と同じ計算量だなんて驚きです。

LU分解と前進後退代入から逆行列を求める方法

np.linalg.invでは連立一次方程式の計算を利用して逆行列を求めるといいました。これは単位行列$E$を右辺とした連立一次方程式を解くことを指しています。つまり以下の方程式です(右辺と解$X$が行列になりますが、単純に列の分だけ解くべき方程式が増えたと思えばOKです)。
$$A X = E.$$
この方程式を解くと、$X = A^{-1}$となるのがわかりますね。

この方法の前進後退代入の計算量は乗除算$n(2n^2+1)/3$回、加減算$n(n-1)(4n-5)/6$回となります(この計算量の計算は結構大変…)。
LU分解の計算量との合計は乗除算が$n^3 + n- 1$回、加減算が$n(n-1)^2$回となります。掃き出し法と比べて乗除算が$n-1$回増えますが、$n$が大きくなれば無視できる程度の差です。

計算量のまとめ

計算量についてまとめると、以下のようになります。

方法 乗除算 加減算
掃き出し法による逆行列の計算 $n^3$ $n(n-1)^2$
行列とベクトルの積 $n^2$ $n(n-1)$
LU分解 $(n-1)(n^2+n+3)/3$ $n(n-1)(2n-1)/6$
前進後退代入 $n^2$ $n(n-1)$
LU分解+前進後退代入による逆行列の計算 $n^3+n-1$ $n(n-1)^2$

LU分解と前進後退代入によって$Ax=b$を解いた場合の計算量では$n^3$に$1/3$がかかっていますから、「逆行列を求める+ベクトルとの積を計算する」の場合に比べて$1/3$程度計算量が減ることがわかります。

NumPyで実験

実際にNumPyで計算時間を比較してみましょう。
以下のようにして行列とベクトルを作ります。

import numpy as np  
A = np.random.rand(1000, 1000)  
b = np.random.rand(1000)  

次に、計算にかかった時間をそれぞれ測ります。

  • 逆行列を求める+ベクトルとの積を計算

    %%timeit  
    inv_x = np.dot(np.linalg.inv(A), b)   

    結果:80.8 ms ± 4.29 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

  • 連立一次方程式を解く

    %%timeit  
    solve_x = np.linalg.solve(A, b)   

    結果:27.7 ms ± 1.21 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

おおよそ3倍くらいの差がつきましたね!

まとめ

連立一次方程式の形に落とし込める場合には、逆行列を求めずに、連立一次方程式として解いてあげましょう!実は、計算量が減ることで数値誤差が増えづらくなり、精度面も連立一次方程式のほうが有利と言われています。色々な面で逆行列を計算するメリットはないのです。

参考書籍:伊理 正夫、藤野 和建.数値計算の常識.

技術ブログをはじめよう Qrunch(クランチ)は、プログラマの技術アプトプットに特化したブログサービスです
駆け出しエンジニアからエキスパートまで全ての方々のアウトプットを歓迎しております!
or 外部アカウントで 登録 / ログイン する
クランチについてもっと詳しく

この記事が掲載されているブログ

@opqrstuvcutの技術とか技術以外のブログ 最近フリーランスになりました

よく一緒に読まれる記事

0件のコメント

ブログ開設 or ログイン してコメントを送ってみよう
目次をみる
技術ブログをはじめよう Qrunch(クランチ)は、プログラマの技術アプトプットに特化したブログサービスです
or 外部アカウントではじめる
10秒で技術ブログが作れます!