python - A 不可逆，但 A * 1.00001 是

Question

我有一个奇怪的问题，我认为它来自数字准确性(但我不确定)。

我有一个稀疏矩阵(数据在最后)

x
Out[25]: 
<11250x11250 sparse matrix of type '<class 'numpy.float64'>'
    with 44547 stored elements in Compressed Sparse Column format>

和一个向量

z = np.zeros((x.shape[0]))
z[idx] = 0.1

我正在努力解决

from scipy.sparse.linalg import spsolve
F = spsolve(x, z)

但是我明白了

scipy.sparse.linalg.dsolve.linsolve.MatrixRankWarning: Matrix is exactly singular

我很确定 x 应该 允许反转，我将在后面描述它。然而，指示性的已经是

spsolve(x* 1.00000001, z)
Out[53]: array([-0., -0., -0., ..., -0., -0., -0.])
spsolve(x/ 1.00000001, z)
Out[54]: array([-0., -0., -0., ..., -0., -0., -0.])

如果 x 不可逆，则两者都不应该是整个矩阵的线性缩放。

那么这里发生了什么？这些是 float 问题吗？

---

背景

我有一个矩阵A，它满足与转移矩阵相对应的连续时间生成矩阵的典型属性:

• 对角线为负
• 非对角线是非负的
• 行加起来是一

然后我想计算平稳分布。为此，我创建了一个零向量 z。对于任意索引 idx，我设置 z[idx] = 0.1。我修复了 A[idx] = 0。然后，(未归一化的)平稳分布 g 可以求解为

A.T g = z

在上一题中，A[idx] = 0 的调整后的A 被标记为x。

---

数据

我不确定如何提供x 背后的数据。当矩阵非常大时会出现问题，所以我不能只将内容粘贴到这里。数据生成问题冗长复杂，也不适合提供。

我想出的最好办法是存储和上传(稀疏)npz ，可以通过 scipy.sparse.load_npz 加载。如果有更好的方法，请告诉我。

---

我还绘制了与 A 关联的最大奇异值:

但是，尚不清楚这些奇异值的单峰性在多大程度上可以预测数值稳定性。接下来是我设法求逆的矩阵的最大奇异值，它们看起来同样平坦:

Answer 1

您可以尝试的一项健全性检查是计算 A 矩阵的特征值谱。随机过程的平稳分布本质上是具有特征值单位的转移矩阵的特征向量。然而，为了有一个定义明确的平衡分布，必须有一个唯一的单位特征值，并且所有其他特征值必须具有小于一的(实部)。因此，如果你计算 A 的奇异值分解，你获得的对角矩阵将揭示一些关于平衡状态的稳定性，以及系统衰减到该平衡的多个时间尺度。如果您的随机过程表现良好，则绘制奇异值图以降序排序的数量级应该在 n=0 处有一个孤立的峰值。但是，如果下降更为缓慢，则表明系统可能只会非常缓慢地稳定到平衡状态，并且随着问题维度的增加可能会出现数值病态。

另一件值得检查的事情 Ornstein-Uhlenbeck过程是，当您增加系统的维度时，您可能需要减少转移矩阵的有效时间步长。这让人想起需要减少抛物型微分方程数值解的时间步长，以对应于描述扩散过程的网格的空间长度尺度，如果要避免数值不稳定性的话。