生成一个pseduo-random正定矩阵

Question

我想测试一下我用C++编写的简单的Cholesky代码.所以我生成一个随机的低三角形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.我猜测随机性会产生秩不足的矩阵.我对吗？

更新:

我忘了提到以下Matlab代码似乎正如下面Malife所指出的那样:

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

差异是L在第一个代码中是低三角形而不是第二个代码.为什么要这么重要？

在阅读了用于生成正半定矩阵的简单算法后,我还尝试了以下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上.有任何想法吗？

谢谢,
Nilesh.

Answer 1

问题不在于胆怯分解.问题在于随机矩阵L. rand(N,N)比条件要好得多tril(rand(N,N)).看到这一点,比较cond(rand(N,N))来cond(tril(rand(N,N))).我得到了类似于1e3第一个和1e19第二个的东西,因此第二个矩阵的条件数要高得多,并且计算在数值上将不太稳定.这将导致在病态情况下得到一些小的负特征值 - 看这个使用的特征值eig(),一些小的将是负的.

所以我建议使用rand(N,N)生成数值稳定的随机矩阵.

顺便说一句,如果你对发生这种情况的理论感兴趣,你可以看看这篇论文:

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