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

Non*_*biz 28 numpy matrix linear-algebra scipy svd

我试图找到给定矩阵的零空间(Ax = 0的解空间).我找到了两个例子,但我似乎无法工作.而且,我无法理解他们为实现目标所做的工作,所以我无法调试.我希望有人能够指引我完成这件事.

文档页面(numpy.linalg.svd,和numpy.compress)对我来说是不透明的.我通过创建矩阵C = [A|0],找到减少的行梯形形式并按行求解变量来学习.在这些例子中,我似乎无法遵循它是如何完成的.

感谢您的帮助!

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

A = matrix([
    [2,3,5],
    [-4,2,3]
    ])  
Run Code Online (Sandbox Code Playgroud)

方法(在这里这里找到):

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)
Run Code Online (Sandbox Code Playgroud)

当我尝试它时,我得到一个空矩阵:

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)
>>> 
Run Code Online (Sandbox Code Playgroud)

Idr*_*Idr 18

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]])]
Run Code Online (Sandbox Code Playgroud)

  • 我相信这种方法仅适用于整数,但不适用于浮点数:(。 (3认同)

Wok*_*Wok 9

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

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


Bas*_*ork 6

它似乎对我有用:

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

  • @Nona Urbiz - 它返回一个空矩阵,因为你没有放入一行零,就像Bashwork(和维基百科)在上面做的那样.另外,零空间值返回(`[-0.33,-0.85,0.52]`)被归一化,使得向量的大小为1.维基百科例子不被标准化.如果你只是把`N = NULL(A)`和看看`N/n.max()',你会得到`[-.0625,-1.625,1]`. (7认同)
  • @Bashwork - 我怎么知道以编程方式添加一行零?矩阵必须是正方形吗? (3认同)

spu*_*ick 5

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

实现遵循规范的 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
Run Code Online (Sandbox Code Playgroud)