生成一个 pseduo-随机正定矩阵

generating a pseduo-random positive definite matrix

本文关键字：随机 pseduo- 一个      更新时间：2023-10-16

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

在阅读了生成正半定矩阵的简单算法后，我还尝试了以下内容：

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会发生这种情况。有什么想法吗？

谢谢
尼莱什。

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

如前所述，三角矩阵的特征值位于对角线上。因此，通过做

L=tril(rand(n))

您确保 eig(L( 仅产生正值。您可以通过在对角线上添加足够大的正数来改善 L*L' 的条件数，例如

L=L+n*eye(n)
L

*L' 是正定的，条件良好：

> cond(L*L')
ans =
1.8400

要在 MATLAB 中生成随机正定矩阵，您的代码应为：

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

而且您确实应该让特征值大于零，p为零。

你问的是下三角形的情况。让我们看看会发生什么，以及为什么会有问题。这通常是一件好事，看看测试用例。

对于简单的 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
相关文章：