5

C++で書いた簡単なコレスキーコードをテストしたかったのです。したがって、ランダムな下三角Lを生成し、その転置を乗算してAを生成します。

A = L * Lt;

しかし、私のコードはAを因数分解できません。そこで、Matlabでこれを試しました。

N=200; L=tril(rand(N, N)); A=L*L'; [lc,p]=chol(A,'lower'); p

これはゼロ以外のpを出力します。これは、MatlabもAを因数分解できないことを意味します。ランダム性によってランクが不足している行列が生成されると思います。私は正しいですか？

アップデート：

以下のMalifeが指摘しているように、次のMatlabコードが機能しているように見えることを忘れました。

N=200; L=rand(N, N); A=L*L'; [lc,p]=chol(A,'lower'); p

違いは、Lは最初のコードでは下三角であり、2番目のコードではないことです。なぜそれが重要なのですか？

また、正値半数行列を生成するための簡単なアルゴリズムを読んだ後、scipyで次のことを試しました。

from scipy import random, linalg
A = random.rand(100, 100)
B = A*A.transpose()
linalg.cholesky(B)

しかし、それはエラーになります：

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python2.7/dist-packages/scipy/linalg/decomp_cholesky.py", line 66, in cholesky
    c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=True)
  File "/usr/lib/python2.7/dist-packages/scipy/linalg/decomp_cholesky.py", line 24, in _cholesky
    raise LinAlgError("%d-th leading minor not positive definite" % info)
numpy.linalg.linalg.LinAlgError: 2-th leading minor not positive definite

なぜそれがscipyで起こっているのか分かりません。何か案は？

ありがとう、
ナイルシュ。

4

4 に答える 4

6

問題はコレスキー分解ではありません。問題はランダム行列にありLます。 rand(N,N)よりもはるかに良い条件ですtril(rand(N,N))。これを確認するには、と比較cond(rand(N,N))してくださいcond(tril(rand(N,N)))。1e31番目と2番目のようなものが得1e19られたので、2番目の行列の条件数ははるかに高く、計算は数値的に不安定になります。これにより、悪条件の場合にいくつかの小さな負の固有値が得られます。これを使用して固有値を確認するにはeig()、いくつかの小さな固有値が負になります。

したがってrand(N,N)、数値的に安定したランダム行列を生成するために使用することをお勧めします。

ところで、なぜこれが起こるのかという理論に興味があるなら、この論文を見ることができます：

http://epubs.siam.org/doi/abs/10.1137/S0895479896312869

于 2012-09-07T19:09:38.647 に答える
1

前に述べたように、三角行列の固有値は対角線上にあります。したがって、

L=tril(rand(n))

eig（L）が正の値のみを生成することを確認しました。対角線に十分な大きさの正の数を追加することにより、L *L'の条件数を改善できます。

L=L+n*eye(n)

L * L'は正定値であり、条件が整っています。

> cond(L*L')

ans =

1.8400
于 2012-09-08T10:09:15.937 に答える
0

MATLABでランダムな正定値行列を生成するには、コードは次のようになります。

N=200;
L=rand(N, N); 
A=L*transpose(L); 
[lc,p]=chol(A,'lower');
eig(A)
p

そして、実際には、固有値がゼロより大きく、ゼロでpある必要があります。

于 2012-09-07T17:49:37.553 に答える
0

下三角の場合についてお伺いします。何が起こるのか、なぜ問題があるのか​​を見てみましょう。これは、テストケースを確認するために行うのがよい場合がよくあります。

単純な5x5マトリックスの場合、

L = tril(rand(5))
L =
      0.72194            0            0            0            0
     0.027804      0.78422            0            0            0
      0.26607     0.097189      0.77554            0            0
      0.96157      0.71437      0.98738      0.66828            0
     0.024571     0.046486      0.94515      0.38009     0.087634

eig(L)
ans =
     0.087634
      0.66828
      0.77554
      0.78422
      0.72194

もちろん、三角行列の固有値は対角要素にすぎません。randによって生成される要素は常に0から1の間であるため、平均して約1/2になります。おそらく、Lの行列式の分布を調べることが役立つでしょう。log（det（L））の分布を検討することをお勧めします。行列式は単に対角要素の積になるため、対数は対角要素の対数の合計になります。（はい、行列式は特異性の尺度としては不十分ですが、log（det（L））の分布は簡単に計算でき、条件数の分布について考えるのが面倒です。）

ああ、しかし、一様確率変数の負の対数は指数変量であり、この場合はラムダ= 1の指数です。間隔（0,1）からのn個の一様確率変数のセットの対数の合計は中心極限定理はガウス分布です。その合計の平均は-nになります。したがって、このようなスキームによって生成される下三角nxn行列の行列式は、exp（-n）になります。nが200の場合、MATLABは次のように指示します。

exp(-200)
ans =
   1.3839e-87

したがって、かなりのサイズのマトリックスの場合、コンディショニングが不十分であることがわかります。さらに悪いことに、積L * L'を形成するとき、それは一般に数値的に特異になります。同じ引数が条件数に適用されます。したがって、20x20の行列の場合でも、そのような下三角行列の条件数がかなり大きいことを確認してください。次に、行列L * L'を作成すると、条件は期待どおりに2乗されます。

L = tril(rand(20));

cond(L)
ans =
   1.9066e+07

cond(L*L')
ans =
   3.6325e+14

完全なマトリックスの場合、どれだけ優れているかを確認してください。

A = rand(20);

cond(A)
ans =
       253.74

cond(A*A')
ans =
        64384
于 2012-09-07T23:17:04.967 に答える