Python sage:我如何计算化学计量矩阵的零空间(内核)？

Question

在拼命尝试从Matlab切换到python时,我遇到了以下问题:

在Matlab中,我能够定义一个矩阵,如:

N = [1  0  0  0 -1 -1 -1  0  0  0;% A
     0  1  0  0  1  0  0 -1 -1  0;% B
     0  0  0  0  0  1  0  1  0 -1;% C
     0  0  0  0  0  0  1  0  0 -1;% D
     0  0  0 -1  0  0  0  0  0  1;% E
     0  0 -1  0  0  0  0  0  1  1]% F

然后可以通过以下方式计算合理基础nullspace(内核):

K_nur= null(N,'r')

和标准正交基础如下:

K_nuo= null(N)

这输出如下:

N =

 1     0     0     0    -1    -1    -1     0     0     0
 0     1     0     0     1     0     0    -1    -1     0
 0     0     0     0     0     1     0     1     0    -1
 0     0     0     0     0     0     1     0     0    -1
 0     0     0    -1     0     0     0     0     0     1
 0     0    -1     0     0     0     0     0     1     1


K_nur =

 1    -1     0     2
-1     1     1     0
 0     0     1     1
 0     0     0     1
 1     0     0     0
 0    -1     0     1
 0     0     0     1
 0     1     0     0
 0     0     1     0
 0     0     0     1


K_nuo =

 0.5933    0.1332    0.3070   -0.3218
-0.0930    0.0433    0.2029    0.7120
 0.1415    0.0084    0.5719    0.2220
 0.3589    0.1682   -0.0620    0.1682
-0.1628    0.4518    0.3389   -0.4617
 0.3972   -0.4867    0.0301   -0.0283
 0.3589    0.1682   -0.0620    0.1682
-0.0383    0.6549   -0.0921    0.1965
-0.2174   -0.1598    0.6339    0.0538
 0.3589    0.1682   -0.0620    0.1682

我一直试图在Python SAGE中复制它,但到目前为止,我没有成功.我的代码看起来像这样:

st1= matrix([
[ 1, 0, 0, 0,-1,-1,-1, 0, 0, 0],
[ 0, 1, 0, 0, 1, 0, 0,-1,-1, 0],
[ 0, 0, 0, 0, 0, 1, 0, 1, 0,-1],
[ 0, 0, 0, 0, 0, 0, 1, 0, 0,-1],
[ 0, 0, 0,-1, 0, 0, 0, 0, 0, 1],
[ 0, 0,-1, 0, 0, 0, 0, 0, 1, 1]])

print st1

null2_or= transpose(st1).kernel()
null2_ra= transpose(st1).kernel().basis()

print "nullr2_or"
print null2_or
print "nullr2_ra"
print null2_ra

注意:在阅读了一些关于此的教程之后介绍了转置,并且与SAGE的性质有关,它自动从左侧计算内核(在这种情况下根本不产生任何结果).

现在的问题是:它会给我打印一些东西......但不是正确的.

输出如下:

sage: load stochiometric.py
[ 1  0  0  0 -1 -1 -1  0  0  0]
[ 0  1  0  0  1  0  0 -1 -1  0]
[ 0  0  0  0  0  1  0  1  0 -1]
[ 0  0  0  0  0  0  1  0  0 -1]
[ 0  0  0 -1  0  0  0  0  0  1]
[ 0  0 -1  0  0  0  0  0  1  1]
nullr2_or
Free module of degree 10 and rank 4 over Integer Ring
Echelon basis matrix:
[ 1  0  0  1  0  0  1  1 -1  1]
[ 0  1  0  1  0 -1  1  2 -1  1]
[ 0  0  1 -1  0  1 -1 -2  2 -1]
[ 0  0  0  0  1 -1  0  1  0  0]
nullr2_ra
[
(1, 0, 0, 1, 0, 0, 1, 1, -1, 1),
(0, 1, 0, 1, 0, -1, 1, 2, -1, 1),
(0, 0, 1, -1, 0, 1, -1, -2, 2, -1),
(0, 0, 0, 0, 1, -1, 0, 1, 0, 0)
]

仔细检查后,您可以看到生成的内核矩阵(nullspace)看起来相似,但不一样.

有没有人知道我需要做什么来获得与Matlab相同的结果,如果可能的话,如何获得正交结果(在Matlab中称为K_nuo).

我试图浏览教程,文档等,但到目前为止,没有运气.

Answer 1

使用SAGE内置函数可能有办法做到这一点; 我不确定.

但是,如果基于numpy/python的解决方案可以,那么:

import numpy as np

def null(A, eps=1e-15):
    """   
    http://mail.scipy.org/pipermail/scipy-user/2005-June/004650.html
    """
    u, s, vh = np.linalg.svd(A)
    n = A.shape[1]   # the number of columns of A
    if len(s)<n:
        expanded_s = np.zeros(n, dtype = s.dtype)
        expanded_s[:len(s)] = s
        s = expanded_s
    null_mask = (s <= eps)
    null_space = np.compress(null_mask, vh, axis=0)
    return np.transpose(null_space)

st1 = np.matrix([
    [ 1, 0, 0, 0,-1,-1,-1, 0, 0, 0],
    [ 0, 1, 0, 0, 1, 0, 0,-1,-1, 0],
    [ 0, 0, 0, 0, 0, 1, 0, 1, 0,-1],
    [ 0, 0, 0, 0, 0, 0, 1, 0, 0,-1],
    [ 0, 0, 0,-1, 0, 0, 0, 0, 0, 1],
    [ 0, 0,-1, 0, 0, 0, 0, 0, 1, 1]])

K = null(st1)
print(K)

产生正交零空间:

[[ 0.59330559  0.13320203  0.30701044 -0.32180406]
 [-0.09297005  0.04333798  0.20286425  0.71195719]
 [ 0.14147329  0.00837169  0.5718718   0.22197807]
 [ 0.35886225  0.16816832 -0.06199711  0.16817506]
 [-0.16275558  0.45177747  0.33887617 -0.46165922]
 [ 0.39719892 -0.48674377  0.03013138 -0.0283199 ]
 [ 0.35886225  0.16816832 -0.06199711  0.16817506]
 [-0.03833668  0.65491209 -0.09212849  0.19649496]
 [-0.21738895 -0.15979664  0.63386891  0.05380301]
 [ 0.35886225  0.16816832 -0.06199711  0.16817506]]

这确认了列具有null space属性:

print(np.allclose(st1*K, 0))
# True

这证实了这K是正交的:

print(np.allclose(K.T*K, np.eye(4)))
# True