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)
你得到了矩阵的SVD分解A.s是特征值的向量.您对几乎为零的特征值感兴趣(参见$ A*x =\lambda*x $,其中$\abs(\ lambda)<\ epsilon $),它由逻辑值向量给出null_mask.
然后,从列表中提取vh对应于几乎为零的特征值的特征向量,这正是您正在寻找的:一种跨越零空间的方法.基本上,您提取行然后转置结果,以便获得具有特征向量作为列的矩阵.
它似乎对我有用:
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)
截至去年(2017 年),scipy 现在null_space在scipy.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)