0

我有奇异矩阵A（10 * 10），它是秩不足的（rank = 9），我有向量b，它在A的范围空间中。现在我对Ax = b的一些解决方案感兴趣。具体来说，这里是我的 A

array([[ 0.        ,  0.        ,  0.        ,  0.86826141,  0.        ,
             0.        ,  0.88788426,  0.        ,  0.4089203 ,  0.88134901],
           [ 0.        ,  0.        ,  0.46416372,  0.        ,  0.        ,
             0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
           [ 0.        ,  0.        ,  0.        ,  0.        ,  0.31303966,
             0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
           [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
             0.        ,  0.3155742 ,  0.        ,  0.64059294,  0.        ],
           [ 0.        ,  0.        ,  0.        ,  0.        ,  0.51349938,
             0.        ,  0.        ,  0.        ,  0.53593509,  0.        ],
           [ 0.        ,  0.01252787,  0.        ,  0.6870415 ,  0.        ,
             0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
           [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
             0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
           [ 0.        ,  0.        ,  0.        ,  0.16643105,  0.        ,
             0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
           [ 0.08626592,  0.        ,  0.        ,  0.        ,  0.        ,
             0.        ,  0.        ,  0.        ,  0.        ,  0.66939531],
           [ 0.43694586,  0.        ,  0.        ,  0.        ,  0.        ,
             0.95941661,  0.        ,  0.52936733,  0.79687149,  0.81463887]])

b 是使用 生成的A.dot(np.ones(10))。现在我想使用 lu 分解来解决这个问题，为此我做了以下

lu_fac=scipy.linalg.lu_factor(X)
scipy.linalg.lu_solve(lu_fac,b)

这使

array([ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan])

在这种情况下， lu_factor 似乎也可以正常工作（有时它确实会给出运行时警告，说“对角线数 %d 正好为零。奇异矩阵”）。为了完整起见，这里是从 lu_factor 验证 PLU 的代码与 A 相同：

L=np.tril(lu_fac[0])
np.fill_diagonal(L,1)
U=np.triu(lu_fac[0])
perm=np.arange(10)
ipiv=lu_factor[1]
for i in range(10):
  temp=perm[i]
  perm[i]=perm[ipiv[i]]
  perm[ipiv[i]]=temp
np.allclose(X[perm,:],L.dot(U))

现在我知道我的矩阵是奇异的，并且我的问题有无数种解决方案。但是我对任何解决方案都感兴趣，我只是很困惑为什么 lu 分解失败，它不能将自由变量设置为 0 并按照我们的教导找到一些解决方案吗？还有什么是运行时警告 “对角线数 %d 正好为零。奇异矩阵”。注意我对解决这个问题的 svd/qr 方法不感兴趣，我只是想知道为什么 lu 对于奇异矩阵会失败。非常感谢任何建议。谢谢。

4

2 回答 2

1 
0 / lu_fac[0][9, 9]

返回nan是因为该条目 - U 的最后一个对角线条目为零。所以这nan成为第 9 个变量的值。然后将其代入上面的等式中，自然而然，其余的nan也是 。SciPy 的 LU 代码，或者更确切地说是它包装的 Fortran 代码，不是为秩亏矩阵设计的，因此它不会为无法确定的变量弥补值。

还有什么是运行时警告“对角线数 %d 正好为零。奇异矩阵”。

警告很明确：算法检测到了一个奇异矩阵，这不是预期的。它还告诉您该实现不适用于奇异矩阵。

有向量 b 在 A 的范围空间中

理论上是这样。在实践中，由于浮点运算中固有的错误，人们无法确定秩亏矩阵的范围空间中是否存在任何内容。您可以计算b = A.dot(...)然后尝试求解 Ax=b，由于处理浮点数时引入的错误，因此不会有解决方案。

顺便说一句：您提到每个方阵都存在 PLU 分解，但 SciPy 不一定旨在计算它。例如，

scipy.linalg.lu_factor(np.array([[0, 1], [0, 0]]))

返回一个带有 NaN 的矩阵。在您的情况下，当尝试找到解决方案并遇到因子 U 的零对角元素时，NaN 会稍后出现。

于 2017-10-04T02:44:52.433 回答
0

正如这里提到的，一个矩阵有一个 LU 分解当且仅当rank(A11) + k >= rank([A11 A12]) + rank([A11 A21])。在您的情况下，rank(A11) = 3、k = 5和rank([A11 A12]) + rank([A11 A21]) = 9。所以你的矩阵不满足条件并且没有LU分解。

于 2017-10-03T23:04:23.470 回答