より効率的な方法は、それが対称で正の半定理であることを知っている行列を反転する

Question

私はnumpyとpythonで反転共分散行列です。共分散行列は対称で正の半定理である。

私は、対称陽性半限定マトリックスに最適化されたアルゴリズムがnumpy.linalg.inv()より速く存在するかどうか疑問に思っていました（もちろん、その実装にはPythonから簡単にアクセスできます）。私はnumpy.linalgで何かを見つけたり、ウェブを検索したりすることはできませんでした。

EDIT：

@yixuanによって観察されたように、半正定値行列は、一般的には厳密に可逆ではありません。私は私の場合、正の明確な行列を得たことを確認したので、正の確定性のために働く答えを受け入れました。とにかく、LAPACKの低レベルのルーチンでは、DSY*ルーチンがsimmetric/hermitian行列のために最適化されていますが、scipyには見つからないようです（多分インストールされたバージョンの問題です）。

Answer 1

APIはまだ存在しませんが、低レベルのLAPACK ?POTRIルーチンファミリを使用することができます。

Docstring:  
inv_a,info = dpotri(c,[lower,overwrite_c]) 

Wrapper for ``dpotri``. 

Parameters 
---------- 
c : input rank-2 array('d') with bounds (n,n) 

Other Parameters 
---------------- 
overwrite_c : input int, optional 
    Default: 0 
lower : input int, optional 
    Default: 0 

Returns 
------- 
inv_a : rank-2 array('d') with bounds (n,n) and c storage 
info : int 
Call signature: sp.linalg.lapack.dpotri(*args, **kwargs) 

最も重要なのはinfo出力され、次のように

sp.linalg.lapack.dpotriのドキュメント文字列です。それがゼロであれば、正の確定性にかかわらず、式を成功裏に解くことを意味する。これは低レベル呼び出しであるため、他のチェックは実行されません。

>>>> M = np.random.rand(10,10) 
>>>> M = M + M.T 
>>>> # Make it pos def 
>>>> M += (1.5*np.abs(np.min(np.linalg.eigvals(M))) + 1) * np.eye(10) 
>>>> zz , _ = sp.linalg.lapack.dpotrf(M, False, False) 
>>>> inv_M , info = sp.linalg.lapack.dpotri(zz) 
>>>> # lapack only returns the upper or lower triangular part 
>>>> inv_M = np.triu(inv_M) + np.triu(inv_M, k=1).T 

また、あなたがスピード

>>>> %timeit sp.linalg.lapack.dpotrf(M) 
The slowest run took 17.86 times longer than the fastest. This could mean that an intermediate result is being cached. 
1000000 loops, best of 3: 1.15 µs per loop 

>>>> %timeit sp.linalg.lapack.dpotri(M) 
The slowest run took 24.09 times longer than the fastest. This could mean that an intermediate result is being cached. 
100000 loops, best of 3: 2.08 µs per loop 

>>>> III = np.eye(10) 

>>>> %timeit sp.linalg.solve(M,III, sym_pos=1,overwrite_b=1) 
10000 loops, best of 3: 40.6 µs per loop 

を比較するのであれば、あなたはかなり無視できない速利益を得ます。複素数で作業している場合は、代わりにzpotriを使用する必要があります。

質問はあなたが逆を必要とするかどうかです。 solve(B,A)がそれに適しているため、B⁻¹ * Aを計算する必要がある場合はおそらくありません。

Answer 2

私の知るところでは、対称行列の標準行列逆関数はありません。一般に、ソルバのスピードアップを得るには、スパースネスなどの制約が必要です。しかし、scipy.linalgを見ると、エルミート（対称）行列に最適化されたいくつかの固有値ルーチンがあることがわかります。例えば

、私はランダム200x200の密行列を生成し、私が得る固有値解決しません：だから何逆に利点はなく

from scipy.linalg import inv,pinvh,eig,eigh 
B = np.rand(200,200) 
B = B+B.T 
%timeit inv(B) 
1000 loops, best of 3: 915 µs per loop 

%timeit pinvh(B) 
100 loops, best of 3: 6.93 ms per loop 

を：

%timeit eig(B) 
10 loops, best of 3: 39.1 ms per loop 

%timeit eigh(B) 
100 loops, best of 3: 4.9 ms per loop 

固有値のクールな8倍速の高速化。

あなたの行列が疎な場合は、いくつかのソルバを持つscipy.sparse.linalgをチェックする必要があります.bicgやcgのようないくつかはエルミート行列を必要とするため、高速かもしれません。しかし、行列が疎である場合にのみ妥当であり、特定の解ベクトルbを解くだけであり、行列構造によっては実際には高速ではないかもしれません。あなたは本当に見つけるためにそれをベンチマークする必要があります。

私は同様の質問をC++ solversについて尋ねましたが、実際にはアプリケーションに依存しており、問題のための最良のソルバーを選ぶ必要があります。

Answer 3

あなたは共分散行列は、それ以外のマトリックスが反転ではない、むしろ半 -definiteより正定（P.D.）であることを確認する必要がありますまず第一に。

p.d.行列反転は、通常、LU分解によって行われるが、p。一般に計算コストを削減するコレスキー分解を用いることができる。 scipyのダウンロードで

、linalg.solve()機能は以下...行列がp.dであると仮定し、パラメータsym_posを持っている簡単な例です：

import numpy as np 
from scipy import linalg 
import time 

def pd_inv(a): 
    n = a.shape[0] 
    I = np.identity(n) 
    return linalg.solve(a, I, sym_pos = True, overwrite_b = True) 

n = 50 
A = np.random.rand(n, n) 
B = A.dot(A.T) 

start = time.clock() 
for i in xrange(0, 10000): 
    res = linalg.inv(B) 

end = time.clock() 
print end - start 
## 1.324752 

start = time.clock() 
for i in xrange(0, 10000): 
    res = pd_inv(B) 

end = time.clock() 
print end - start 
## 1.109778 

私のマシンでは、pd_inv()は小さな行列（〜100×100）のためのいくつかの利点があります。大きな行列の場合はほとんど違いがなく、時にはpd_inv()がさらに遅くなります。