Python (NumPy, SciPy)，寻找矩阵的零空间

Posted 2023-03-12

【中文标题】Python (NumPy, SciPy)，寻找矩阵的零空间

【英文标题】：Python (NumPy, SciPy), finding the null space of a matrix

【发布时间】：2011-08-18 20:23:52

【问题描述】：

我试图找到给定矩阵的零空间（Ax=0 的解空间）。我找到了两个例子，但我似乎都无法工作。此外，我不明白他们在做什么才能到达那里，所以我无法调试。我希望有人能够指导我完成此操作。

文档页面（numpy.linalg.svd 和 numpy.compress）对我来说是不透明的。我学会了通过创建矩阵C = [A|0] 来做到这一点，找到减少的行梯形形式并逐行求解变量。在这些示例中，我似乎无法理解它是如何完成的。

感谢您的所有帮助！

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

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

方法（found here和here）：

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)
>>> 

【问题讨论】：

您链接到的***页面实际上很好地解释了为什么在处理浮点值时应该使用 SVD 来计算矩阵的零空间（或求解）。 en.wikipedia.org/wiki/…您描述的方法（逐行求解变量）会放大任何舍入误差等（这与您几乎不应该显式计算矩阵的逆矩阵的原因相同......）

【参考方案1】：

截至去年（2017 年），scipy 现在在 scipy.linalg 模块 (docs) 中内置了 null_space 方法。

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

【参考方案2】：

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]])]

【讨论】：

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

我相信这种方法只适用于整数，但不适用于浮点数:(.

@Adriaan 我刚刚用花车测试过，它对我有用。

【参考方案3】：

一种更快但数值稳定性较差的方法是使用a rank-revealing QR decomposition，例如scipy.linalg.qr 和pivoting=True：

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]]

【讨论】：

我试过了，在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] ]) ```我错过了什么吗？ （我很擅长线性代数，但我不知道这个方法是什么——我只是在尝试）

【参考方案4】：

你的方法几乎是正确的。问题是函数 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.]]

【讨论】：

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

【参考方案5】：

对我来说似乎工作正常：

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

【讨论】：

我确定我遗漏了什么，但 Wikipedia 建议值应该是 [ [-.0625], [-1.625], [1] ]?

此外，它为我返回了一个空矩阵[]。有什么问题？

@Nona Urbiz - 它返回一个空矩阵，因为您没有像上面的 Bashwork（和***）那样放入一行零。此外，返回的零空间值 ([-0.33, -0.85, 0.52]) 被归一化，因此向量的大小为 1。***示例未归一化。如果您只使用n = null(A) 并查看n / n.max()，您将获得[-.0625, -1.625, 1]。

@Bashwork - 我怎么知道以编程方式添加一行零？矩阵必须是正方形吗？

【参考方案6】：

你得到矩阵A 的 SVD 分解。 s 是特征值向量。您对几乎为零的特征值感兴趣（参见 $A*x=\lambda*x$ where $\abs(\lambda)null_mask 给出。

然后，您从列表vh 中提取与几乎为零的特征值对应的特征向量，这正是您要寻找的：一种跨越零空间的方法。基本上，您提取行，然后转置结果，以便得到一个以特征向量为列的矩阵。

【讨论】：

非常感谢您抽出宝贵时间回复并帮助我。您的回答对我很有帮助，但我最终接受了 Bashworks 的回答，它让我找到了解决方案。不过，我能够理解解决方案的唯一原因是您的回复。

别担心，我以为你的问题是别的。