5

我想测试我用 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 会发生这种情况。有任何想法吗？

谢谢，
尼莱什。

4

4 回答 4

6

问题不在于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

于 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)) 的分布很容易计算，而且我懒得考虑条件数的分布。）

啊，但是均匀随机变量的负对数是指数变量，在这种情况下是 lambda = 1 的指数。来自区间 (0,1) 的一组 n 个均匀随机数的对数之和将由中心极限定理是高斯的。该总和的平均值将为-n。因此，由这种方案生成的下三角 nxn 矩阵的行列式将是 exp(-n)。当 n 为 200 时，MATLAB 告诉我

exp(-200)
ans =
   1.3839e-87

因此，对于任何可观大小的矩阵，我们可以看到它的条件很差。更糟糕的是，当您形成产品 L*L' 时，它通常是数字单数。相同的参数适用于条件编号。因此，即使是 20x20 矩阵，看到这样的下三角矩阵的条件数也相当大。然后当我们形成矩阵 L*L' 时，条件将按预期平方。

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 回答