了解如何使用 lu 分解求解奇异方阵？

Question

我有秩亏的奇异矩阵 A(10*10)(秩=9)，我有向量 b，它在 A 的 space 范围内。现在我对 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 在这种情况下似乎工作正常（有时它确实会给出运行时间警告说 "Diagonal number %d is exactly zero. Singular matrix"）。为了完整起见，这里是用于验证来自 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 并按照我们所教的那样找到一些解决方案吗？运行时间警告 "Diagonal number %d is exactly zero. Singular matrix" 又是怎么回事？请注意，我对 svd/qr 解决此问题的方法不感兴趣，我只是想知道为什么 lu 对奇异矩阵失败。非常感谢任何建议。谢谢。

Answer 1

如前所述 here，当且仅当 rank(A11) + k >= rank([A11 A12]) + rank([A11 A21]) 时矩阵具有 LU 因式分解。在您的情况下，rank(A11) = 3、k = 5、和 rank([A11 A12]) + rank([A11 A21]) = 9。所以你的矩阵不满足条件，没有LU分解。

Answer 2

0 / lu_fac[0][9, 9]

returns nan 因为该条目 - U 的最后一个对角线条目为零。所以这个nan就变成了第9个变量的值。然后将它代入上面的方程式，自然而然地，剩下的结果也为 nan。 SciPy的LU代码，或者更确切地说是它包装的Fortran代码，不是为秩亏矩阵设计的，所以它不会为无法确定的变量补值。

Also what is the deal with the run time warning "Diagonal number %d is exactly zero. Singular matrix".

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

have vector b which is in range space of A

这是理论上的。实际上，由于浮点运算中固有的错误，人们无法确定秩亏矩阵 space 范围内的任何内容。您可以计算 b = A.dot(...)，然后尝试求解 Ax=b，但由于操作浮点数时引入的错误，不会有解。

顺便说一下：您提到每个方阵都存在 PLU 因式分解，但 SciPy 不一定是为计算它而设计的。例如，

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

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

了解如何使用 lu 分解求解奇异方阵？

Gettings nans on using lu factorization for solving singular square matrix?

python

linear-algebra

scipy

matrix-factorization