4

コレスキー分解をインプレースで実行しようとしていますが、その機能が実際にscipyに実装されていないか、理解できないことがあります。後者の場合に備えて、ここに投稿します。これが私がしていることの簡単な例です：

import numpy
import scipy.linalg

numpy.random.seed(0)
X = numpy.random.normal(size=(10,4))
V = numpy.dot(X.transpose(),X)
R = V.copy()
scipy.linalg.cholesky(R,False,overwrite_a=True)
print V
print R

Rが上三角行列で上書きされるとどうなるかと思います。ただし、このコードを実行すると、VとRは同じになります（コレスキーはRを上書きしません）。私はoverwrite_aの目的を誤解していますか、それとも他の間違いを犯していますか？これがscipyの単なるバグである場合、Pythonでコレスキー分解をインプレースで実行できる回避策または別の方法はありますか？

4

3 に答える 3

4

もう一度試してください

>>> scipy.__version__
'0.11.0'
>>> np.__version__
'1.6.2'

それは完璧に機能します：

X = np.random.normal(size=(10,4))
V = np.dot(X.transpose(),X)
print V
R = V.copy()
print R
C = scipy.linalg.cholesky(R,False,overwrite_a=True)
print V
print R

出力は次のとおりです。

[[ 11.22274635   5.10611692   0.70263037   3.14603088] # V before
 [  5.10611692   8.94518939  -3.17865941   1.64689675]
 [  0.70263037  -3.17865941   7.35385131  -2.23948391]
 [  3.14603088   1.64689675  -2.23948391   8.25112653]]
[[ 11.22274635   5.10611692   0.70263037   3.14603088] # R before
 [  5.10611692   8.94518939  -3.17865941   1.64689675]
 [  0.70263037  -3.17865941   7.35385131  -2.23948391]
 [  3.14603088   1.64689675  -2.23948391   8.25112653]]
[[ 11.22274635   5.10611692   0.70263037   3.14603088] # V after
 [  5.10611692   8.94518939  -3.17865941   1.64689675]
 [  0.70263037  -3.17865941   7.35385131  -2.23948391]
 [  3.14603088   1.64689675  -2.23948391   8.25112653]]
[[ 3.35003677  1.52419728  0.20973811  0.93910339] # R after
 [ 0.          2.57332704 -1.35946252  0.08375069]
 [ 0.          0.          2.33703292 -0.99382158]
 [ 0.          0.          0.          2.52478036]]
于 2013-01-18T22:56:29.760 に答える
2

あなたが十分に勇敢であるならば、あなたは避けscipyて、低レベルの電話に行くことができますlinalg.lapack_lite.dpotrf

import numpy as np

# prepare test data
A = np.random.normal(size=(10,10))
A = np.dot(A,A.T)
L = np.tril(A)

# actual in-place cholesky
assert L.dtype is np.dtype(np.float64)
assert L.flags['C_CONTIGUOUS']
n, m = L.shape
assert n == m
result = np.linalg.lapack_lite.dpotrf('U', n, L, n, 0)
assert result['info'] is 0

# check if L is the desired L cholesky factor
assert np.allclose(np.dot(L,L.T), A)
assert np.allclose(L, np.linalg.cholesky(A))

lapackのDPOTRF、Fortranの呼び出し規約、メモリレイアウトを理解する必要があります。終了時にチェックすることを忘れないでくださいresult['info'] == 0。それでも、コードは1行だけであることがわかります。すべての健全性チェックを破棄し、linalg.choleskyこれによって行われるコピーもより効率的である可能性があります。

于 2013-01-19T12:02:18.530 に答える
2

完全を期すために、ウォーレンのコメントに基づく次のコードで問題を解決します。

import numpy
import scipy.linalg

numpy.random.seed(0)
X = numpy.random.normal(size=(10,4))
V = numpy.dot(X.transpose(),X)
R = V.copy('F')
print V.flags['C_CONTIGUOUS']
print R.flags['F_CONTIGUOUS']
scipy.linalg.cholesky(R,False,overwrite_a=True)
print V
print R

私が行ったのは、Rの基になるストレージ形式をCオーダーからFortranオーダーに変更することだけです。これによりoverwrite_a、期待どおりに動作します。

行列が実数値の場合は、転置をコレスキーに渡すだけで、コピーを作成せずにこれを行うことができます。複素エルミート行列の転置は一般にその共役転置と等しくないため、行列が複素数値の場合、転置修正は機能しないと思います。ただし、最初にマトリックスがFortran順序（R.flags['F_CONTIGUOUS'] == True）を使用していることを確認すれば、問題は発生しないはずです。ただし、その戦略の難しさについては、Numpyでの*デフォルト*のデータ順序（CとFortran）の設定を参照してください。

于 2013-01-22T22:03:51.800 に答える