生成一个伪随机正定矩阵。

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，参考了“生成正半定矩阵的简单算法”：https://dev59.com/XHRB5IYBdhLWcg3wc280

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。

- Nilesh

---

在numpy中（被scipy所使用的库），数组相乘的符号*代表的是逐个元素相乘，而不是代表代数乘法。如果要进行代数乘法计算，应该使用dot函数。因此，你应该写成B = dot(A, A.transpose())或者B = dot(A, A.T)。 - Warren Weckesser

我会生成一个标准高斯矩阵A，并进行 A * A^t 运算。这会产生一个 Wishart 分布的矩阵，对于这种实验来说条件更好（特征值具有 xhi ^ 2 分布，如果我没记错的话）。 - Alexandre C.

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

- Bitwise

---

太好了！我没有想到那一点，虽然算法本身是稳定的，但病态数据导致了问题。这篇论文描述了一个非常有趣的现象，尽管我需要更多的背景才能完全理解它。谢谢！ 另一个问题是，第二种方法（rand（N，N））是否可以证明生成SPD矩阵？我在想它会不会生成一个秩缺失的矩阵，虽然概率可能非常小。你有什么想法吗？ 我不知道SciPy出了什么问题，对于生成的随机矩阵（1e3），条件数似乎还可以。 - Nilesh

使用第二种方法仍然可能得到病态矩阵并获得负特征值，这实际上非常简单。由于它生成随机矩阵，因此您可以基本上得到任何矩阵（0至1之间），包括病态矩阵，具体而言，您可以得到非常类似于tril(rand(N,N))给您的矩阵（概率很小）。例如：矩阵tril(rand(200,200))+0.0001可能可以由rand(200,200)生成，但将是病态的。 - Bitwise

Answer 2

正如之前所说，三角矩阵的特征值位于对角线上。因此，通过进行

L=tril(rand(n))

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

L=L+n*eye(n)

如果L*L'是正定且良好条件的：

> cond(L*L')

ans =

1.8400

- angainor

Answer 3

在MATLAB中生成一个随机正定矩阵，您的代码应该如下：

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

你确实应该让特征值大于零且 p 为零。

- Malife

---

是的，我忘了在原问题中提到这一点，现在已经更新。谢谢！但我不明白为什么如果L是下三角矩阵（如第一个Matlab代码示例），A就不能被分解因式。 - Nilesh

Answer 4

你问的是下三角情况。让我们看看会发生什么，以及为什么会有问题。这通常是一个好主意，去看一个测试案例。

对于一个简单的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是对角线元素的对数之和。(是的，我知道行列式是奇异性的一个差劣度量，但是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

- user85109