36

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

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

感谢您的任何帮助!

这是我的示例矩阵,与维基百科示例相同:

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

6 回答 6

24

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]])]
于 2014-11-10T02:09:03.040 回答
11

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

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

于 2011-05-04T20:05:52.020 回答
11

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

实现遵循规范的 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
于 2018-07-19T14:40:21.513 回答
7

它对我来说似乎工作正常:

A = matrix([[2,3,5],[-4,2,3],[0,0,0]])
A * null(A)
>>> [[  4.02455846e-16]
>>>  [  1.94289029e-16]
>>>  [  0.00000000e+00]]
于 2011-05-04T20:29:13.337 回答
3

一种更快但数值稳定性较差的方法是使用rank-revealing QR 分解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]]
于 2015-11-09T22:15:18.227 回答
2

你的方法几乎是正确的。问题是函数 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.]]
于 2014-10-23T19:25:41.213 回答