c++ - 生成伪随机正定矩阵

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 在第一个代码中是下三角的，而在第二个代码中不是。为什么这很重要？

在阅读 A simple algorithm for generating positive-semidefinite matrices 后，我还用 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 会发生这种情况。有什么想法吗？

谢谢，
尼勒什。

Answer 1

问题不在于 cholesky 分解。问题出在随机矩阵 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