Python（NumPy，SciPy），寻找矩阵的零空间

Question

我正在尝试找到给定矩阵的零空间（Ax=0的解空间）。我已经找到了两个例子，但似乎都无法工作。此外，我不明白他们是如何得出这些结果的，因此我无法进行调试。希望有人能够帮助我完成这个过程。

文档页面（numpy.linalg.svd和numpy.compress）对我来说很难理解。我学会了通过创建矩阵C = [A|0]，找到简化行阶梯形式并按行求解变量来完成此操作。但我似乎无法理解这些示例是如何完成的。

感谢所有提供帮助的人！

以下是我的示例矩阵，与wikipedia example相同：

A = matrix([
    [2,3,5],
    [-4,2,3]
    ])  

方法（在这里找到和在这里）：

import scipy
from scipy import linalg, matrix
def null(A, eps=1e-15):
    u, s, vh = scipy.linalg.svd(A)
    null_mask = (s <= eps)
    null_space = scipy.compress(null_mask, vh, axis=0)
    return scipy.transpose(null_space)

当我尝试运行它时，返回一个空矩阵：

Python 2.6.6 (r266:84292, Sep 15 2010, 16:22:56) 
[GCC 4.4.5] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import scipy
>>> from scipy import linalg, matrix
>>> def null(A, eps=1e-15):
...    u, s, vh = scipy.linalg.svd(A)
...    null_mask = (s <= eps)
...    null_space = scipy.compress(null_mask, vh, axis=0)
...    return scipy.transpose(null_space)
... 
>>> A = matrix([
...     [2,3,5],
...     [-4,2,3]
...     ])  
>>> 
>>> null(A)
array([], shape=(3, 0), dtype=float64)
>>> 

- Nona Urbiz

---

您提供的维基百科页面非常好地解释了，在处理浮点数时为什么要使用奇异值分解（SVD）来计算矩阵的零空间（或解）。您所描述的逐行求解变量的方法会放大任何舍入误差等问题。（这也是为什么您几乎永远不应该显式计算矩阵的逆的原因...） - Joe Kington

Answer 1

Sympy使这很简单。

>>> from sympy import Matrix
>>> A = [[2, 3, 5], [-4, 2, 3], [0, 0, 0]]
>>> A = Matrix(A)
>>> A * A.nullspace()[0]
Matrix([
[0],
[0],
[0]])
>>> A.nullspace()
[Matrix([
[-1/16],
[-13/8],
[    1]])]

- Idr

---

正是我正在寻找的解决方案！ - cxxl

我认为这种方法只适用于整数，而不适用于小数 :(. - Adriaan

@Adriaan 我刚刚用浮点数测试了一下，对我来说可以工作。 - Grant Curell

Answer 2

你可以得到矩阵A的SVD分解。 s 是特征值的向量。你对接近零的特征值感兴趣 (参见$A*x=\lambda*x$，其中$\abs(\lambda)<\epsilon$)，这可以由逻辑值向量null_mask给出。

然后，你从列表vh中提取与接近零特征值相应的特征向量，这正是你要寻找的：一种描述零空间的方法。基本上，你提取行，然后转置结果，这样你就可以得到一个以特征向量为列的矩阵。

- Wok

---

非常感谢您抽出时间回复并帮助我。您的答案对我非常有帮助，但最终我接受了Bashworks的答案，因为它让我找到了解决方案。然而，我之所以能够理解解决方案，是因为您的回复。 - Nona Urbiz

不用担心，我以为你的问题是别的什么东西。 - Wok

Answer 3

截至去年（2017）, scipy现在在scipy.linalg模块中有一个内置的null_space方法 (docs)。

这个implementation遵循了规范的SVD分解，如果你有较旧版本的scipy并且需要自己实现，它非常小（见下文）。然而，如果你是最新的，它就可以使用了。

def null_space(A, rcond=None):
    u, s, vh = svd(A, full_matrices=True)
    M, N = u.shape[0], vh.shape[1]
    if rcond is None:
        rcond = numpy.finfo(s.dtype).eps * max(M, N)
    tol = numpy.amax(s) * rcond
    num = numpy.sum(s > tol, dtype=int)
    Q = vh[num:,:].T.conj()
    return Q

- sputnick

Answer 4

对我来说，看起来这个工作正常：

A = matrix([[2,3,5],[-4,2,3],[0,0,0]])
A * null(A)
>>> [[  4.02455846e-16]
>>>  [  1.94289029e-16]
>>>  [  0.00000000e+00]]

- Bashwork

---

我相信我错过了什么，但维基百科建议这些值应该是 [[-.0625]，[-1.625]，[1]]？ - Nona Urbiz

此外，它对我返回一个空矩阵“[]”。可能出了什么问题？ - Nona Urbiz

因为您没有像Bashwork（以及维基百科）中所示那样输入一行零，所以它返回一个空矩阵。此外，返回的零空间值（[-0.33, -0.85, 0.52]）已被归一化，使向量的大小为1。而维基百科的例子则没有经过归一化处理。如果您只是执行 n = null(A) 并查看 n / n.max()，您将得到 [-.0625, -1.625, 1]。 - Joe Kington

如何以编程的方式添加一行零？这需要矩阵是方阵吗？ - Coder

Answer 5

一种更快但数值稳定性较差的方法是使用透露秩的QR分解，例如使用pivoting=True的scipy.linalg.qr函数：

import numpy as np
from scipy.linalg import qr

def qr_null(A, tol=None):
    Q, R, P = qr(A.T, mode='full', pivoting=True)
    tol = np.finfo(R.dtype).eps if tol is None else tol
    rnk = min(A.shape) - np.abs(np.diag(R))[::-1].searchsorted(tol)
    return Q[:, rnk:].conj()

例如：

A = np.array([[ 2, 3, 5],
              [-4, 2, 3],
              [ 0, 0, 0]])
Z = qr_null(A)

print(A.dot(Z))
#[[  4.44089210e-16]
# [  6.66133815e-16]
# [  0.00000000e+00]]

- ali_m

---

我尝试了一下 Matrix([[2, 2, 1, 0], [2, -2, -1, 1]]) 并得到了不同的答案。使用 sympy，我正确地得到了：[Matrix([ [ 0], [-1/2], [ 1], [ 0]]), Matrix([ [-1/4], [ 1/4], [ 0], [ 1]])]但是使用上述方法，我得到了：array([[0.00000000e+00, 4.16333634e-17], [1.73472348e-16, 0.00000000e+00]])我错过了什么吗？（我在线性代数方面还可以，但我不知道这种方法是什么 - 我只是尝试它） - Grant Curell

Answer 6

您的方法几乎正确。问题在于，函数scipy.linalg.svd返回的s的形状为(K,)，其中K=min(M,N)。因此，在您的示例中，s仅具有两个条目（前两个奇异向量的奇异值）。对您的null函数进行以下更正，即可使其适用于任何大小的矩阵。

import scipy
import numpy as np
from scipy import linalg, matrix
def null(A, eps=1e-12):
...    u, s, vh = scipy.linalg.svd(A)
...    padding = max(0,np.shape(A)[1]-np.shape(s)[0])
...    null_mask = np.concatenate(((s <= eps), np.ones((padding,),dtype=bool)),axis=0)
...    null_space = scipy.compress(null_mask, vh, axis=0)
...    return scipy.transpose(null_space)
A = matrix([[2,3,5],[-4,2,3]])
print A*null(A)
>>>[[  4.44089210e-16]
>>> [  6.66133815e-16]]
A = matrix([[1,0,1,0],[0,1,0,0],[0,0,0,0],[0,0,0,0]])
print null(A)
>>>[[ 0.         -0.70710678]
>>> [ 0.          0.        ]
>>> [ 0.          0.70710678]
>>> [ 1.          0.        ]]
print A*null(A)
>>>[[ 0.  0.]
>>> [ 0.  0.]
>>> [ 0.  0.]
>>> [ 0.  0.]]

- Thomas Wentworth

---

我在工作中一直在使用这段代码，但是我发现了一个问题。eps值为1e-15似乎太小了。特别是考虑矩阵A = np.ones(13,2)。这段代码将报告此矩阵具有秩0的零空间。这是由于scipy.linalg.svd函数报告第二个奇异值高于1e-15。我不太了解该函数背后的算法，但我建议使用eps=1e-12（对于非常大的矩阵可能更低），除非有更多知识的人可以参与讨论。（在无限精度下，第二个奇异值应为0）。 - Thomas Wentworth