

現在シミュレーションに必要な逆行列の計算をpythonで行なっています。  
データは ndarray, dtype=np.float64 型の変数 A に格納されており、A はおよそ 400×400
のマトリックスになっています。以下の方法によって逆行列を求めました。

[code]

    import numpy as np
    Ainv = np.linalg.inv(A)
    
[/code]

ここまでは特にエラーなどなく計算できるのですが、この後

[code]

    AAi = np.dot(A,Ainv)
    AiA = np.dot(Ainv,A)
    
[/code]

などもとの行列との逆行列の積を計算すると、本来単位行列になっていて欲しいのに、全く異なる値が返されます（たとえば`AAi[0][0]=6.68,
AiA[0][0]=5.8e+15`など）。

そこであらかじめ

[code]

    A = A[:100,:100]
    
[/code]

などと範囲を区切って同様の計算を行うと、1となっていてほしい項が9.9997e-01,
0になっていて欲しい項は***e-05~***e-11などになっており、桁の精度は別として正しい計算となっています(具体的には要素数が106×106を超えたあたりから一桁目が合わなくなりました。逆に80×80程度だと６桁程度の精度で単位行列となっています)。

行列の要素数の増加によって精度が悪くなり、結果として元の行列と逆行列の積も単位行列ではなくなってしまっていると考えられるのですが、精度良くこれを求める方法はないでしょうか？
多少計算に時間がかかっても構いません。

環境：  
python2.7.9  
numpy1.14.0(lapack,blasは導入済み)


CC BY-SA 3.0
2018-02-07T14:08:06.390
41548
2018-07-20T22:27:06.093
2018-02-09T11:16:11.420
19453
19453
question
python numpy
pythonの逆行列を高精度に計算する方法


データの性質によるのかも知れませんが、次のように一様分布の乱数行列だと、結果はほぼ1,1になりました。

[code]

    import numpy as np
    
    A = np.random.rand(400, 400)
    
    Ainv = np.linalg.inv(A)
    
    AAi = np.dot(A, Ainv)
    AiA = np.dot(Ainv, A)
    
    
    for i in range(400):
        print(AAi[i][i], AiA[i][i])
    
[/code]

結果

[code]

    (0.999999999999984, 0.9999999999998891)
    (1.0000000000000198, 0.9999999999999999)
    (0.9999999999999643, 1.0000000000002853)
    (1.0000000000000118, 1.0000000000000686)
    (1.0000000000001563, 0.9999999999999708)
    (1.000000000000004, 1.0000000000000784)
    以下略
    
[/code]

ただ、inv()では精度が出ないと言う話は聞いたことがあり、その場合solve()を使えとのことです。

[code]

    Ainv = np.linalg.solve(A, np.eye(400))
    
[/code]

これでもほぼ1,1が得られました。なお、np.show_config()は次の通り。実行時ではなくコンパイル時の設定らしいので、参考になるかわかりません。lapack,
blas, openblas, atlasなどとりあえず全部インストールした状態です。

[code]

    lapack_info:
        libraries = ['lapack', 'lapack']
        library_dirs = ['/usr/lib']
        language = f77
    lapack_opt_info:
        libraries = ['lapack', 'lapack', 'blas', 'blas']
        library_dirs = ['/usr/lib']
        language = c
        define_macros = [('NO_ATLAS_INFO', 1), ('HAVE_CBLAS', None)]
    openblas_lapack_info:
      NOT AVAILABLE
    blas_info:
        libraries = ['blas', 'blas']
        library_dirs = ['/usr/lib']
        define_macros = [('HAVE_CBLAS', None)]
        language = c
    atlas_3_10_blas_threads_info:
      NOT AVAILABLE
    atlas_threads_info:
      NOT AVAILABLE
    atlas_3_10_threads_info:
      NOT AVAILABLE
    atlas_blas_info:
      NOT AVAILABLE
    atlas_3_10_blas_info:
      NOT AVAILABLE
    atlas_blas_threads_info:
      NOT AVAILABLE
    openblas_info:
      NOT AVAILABLE
    blas_mkl_info:
      NOT AVAILABLE
    blas_opt_info:
        libraries = ['blas', 'blas']
        library_dirs = ['/usr/lib']
        language = c
        define_macros = [('NO_ATLAS_INFO', 1), ('HAVE_CBLAS', None)]
    atlas_info:
      NOT AVAILABLE
    atlas_3_10_info:
      NOT AVAILABLE
    lapack_mkl_info:
      NOT AVAILABLE
    mkl_info:
      NOT AVAILABLE
    
[/code]


CC BY-SA 3.0
2018-02-09T06:36:13.300
41605
2018-02-09T06:36:13.300
7837
41548
answer


Python
で計算する際に限らず数値計算一般において、逆行列を直接求めると誤差が大きくなることがあります。理論的には行列の[条件数](https://ja.wikipedia.org/wiki/%E6%9D%A1%E4%BB%B6%E6%95%B0)を使った解析ができ、一般に条件数が大きいと逆行列に含まれうる誤差が大きくなります。今回質問者さんが逆行列を求めたい行列の条件数は
[10¹⁸
程度とのことです](https://ja.stackoverflow.com/questions/41548/python%E3%81%AE%E9%80%86%E8%A1%8C%E5%88%97%E3%82%92%E9%AB%98%E7%B2%BE%E5%BA%A6%E3%81%AB%E8%A8%88%E7%AE%97%E3%81%99%E3%82%8B%E6%96%B9%E6%B3%95#comment42682_41548)が、これはとても大きいです。このため
`np.linalg.inv(A)` で求めた逆行列に誤差が多く含まれても不思議ではありません。

そこで、数値計算を行っていて計算式の途中で逆行列が必要になった場合、それを陽には求めず [LU
分解](https://ja.wikipedia.org/wiki/LU%E5%88%86%E8%A7%A3)などを使って間接的に解だけ求めることがよく行われています。`inv()`
ではなく `solve()`
を使った方が精度が出るのは[このためです](https://stackoverflow.com/a/31257909/5989200)。しかし今回はそうではなく[逆行列そのものが欲しい](https://ja.stackoverflow.com/questions/41548/python%E3%81%AE%E9%80%86%E8%A1%8C%E5%88%97%E3%82%92%E9%AB%98%E7%B2%BE%E5%BA%A6%E3%81%AB%E8%A8%88%E7%AE%97%E3%81%99%E3%82%8B%E6%96%B9%E6%B3%95#comment42682_41548)とのことですので、この方法は使えません。(ただし、本当に逆行列そのものが必要なのかは再度検討なさってください。)

そこで、行列の条件数を減らすように式変形できないか考えます。

簡単な例だと、行列の各成分をスケーリングすることで条件数が小さくなることがあります。たとえば[「数値解の品質：条件、安定性、誤差解析」](http://www.nag-j.co.jp/nagcourse/numerical/quality/section4a/section4a.htm)というページから例を頂くと、以下の行列
A は条件数が大きい (つまり、ill-conditioned である) ような行列です。

[![A = \\left\( \\begin{array}{cc} 2 \\times 10^9 & 10^9 \\ 10^{-9} & 2
\\times 10^{-9} \\ \\end{array}
\\right\)](https://i.stack.imgur.com/B5DxL.png)](https://i.stack.imgur.com/B5DxL.png)

しかし適当な対角行列を掛けることで、行列の条件数を小さくすることができます。

[![B = \\left\( \\begin{array}{cc} 10^{-9} & 0 \\ 0 & 10^9 \\ \\end{array}
\\right\) \\left\( \\begin{array}{cc} 2 \\times 10^9 & 10^9 \\ 10^{-9} & 2
\\times 10^{-9} \\ \\end{array} \\right\) = \\left\( \\begin{array}{cc} 2 & 1
\\ 1 & 2 \\ \\end{array}
\\right\)](https://i.stack.imgur.com/2gxGc.png)](https://i.stack.imgur.com/2gxGc.png)

これがスケーリングです。問題の設定によっては、事前にこのような計算を行っておくことで逆行列を精度良く求められます。(とはいえ LAPACK
は自動的にこのようなことをしているそうなので [要出典]、どこまで変わるかは分かりませんが……。)
また、私は詳しくありませんがスケーリングと関連する話題として preconditioning というものがあるらしく、英語版 Wikipedia の
["Preconditioner"](https://en.wikipedia.org/wiki/Preconditioner)
にはまとまった解説が載っているようです。

更に、元の行列を特異値分解した上で逆行列を求めると上手くいくこともあるようです。これは姉妹サイトの Computational Science Stack
Exchange
で私が質問して教えてもらいました。詳細は[こちらをご覧ください](https://scicomp.stackexchange.com/q/29102/27229)。

以上の手法たちは常に使えるというわけではないですが、問題の性質によっては上手くいくかもしれないので、お試しください。

## 参考

(文中でリンクしたページたちです)

  * [条件数](https://ja.wikipedia.org/wiki/%E6%9D%A1%E4%BB%B6%E6%95%B0) -- Wikipedia
  * [Why does numpy.linalg.solve() offer more precise matrix inversions than numpy.linalg.inv()?](https://stackoverflow.com/q/31256252/5989200) -- Stack Overflow
  * [数値解の品質：条件、安定性、誤差解析](http://www.nag-j.co.jp/nagcourse/numerical/quality/section4a/section4a.htm) -- Numerical Algorithms Group
  * [How to directly compute the inverse of an ill-conditioned dense matrix](https://scicomp.stackexchange.com/q/29102/27229) -- Computational Science Stack Exchange


CC BY-SA 4.0
2018-07-20T15:52:58.213
46795
2018-07-20T22:27:06.093
2018-07-20T22:27:06.093
19110
19110
41548
answer