Python sage：如何计算化学计量矩阵的零空间（内核）？答案

【问题标题】：Python sage: How do I computer a nullspace (kernel) for a stoichiometric matrix?Python sage：如何计算化学计量矩阵的零空间（内核）？
【发布时间】：2012-11-01 10:21:37
【问题描述】：

在从 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

然后可以通过以下方式计算有理基础零空间（内核）：

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

仔细检查后，您可以看到生成的内核矩阵（零空间）看起来相似，但并不相同。

有谁知道我需要做什么才能获得与在 Matlab 中相同的结果，如果可能的话，如何获得正交结果（在 Matlab 中称为 K_nuo）。

我尝试过浏览教程、文档等，但到目前为止，还没有成功。

【问题讨论】：

标签： python matlab matrix sage

【解决方案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]]

这确认列具有空空间属性：

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

这证实了K 是正交的：

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

【讨论】：

哇！非常感谢您的出色回答。我很想给一个以上的赞成票，但不幸的是，我不能......这真的解释了一切。

【解决方案2】：

这样的事情应该可以工作：

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]])
sage: K = st1.right_kernel(); K
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]
sage: M = K.basis_matrix()

gram_schmidt 方法给出了一对矩阵。输入 M.gram_schmidt? 以查看文档。

sage: M.gram_schmidt() # rows are orthogonal, not orthonormal
(
[     1      0      0      1      0      0      1      1     -1      1]
[    -1      1      0      0      0     -1      0      1      0      0]
[  5/12    3/4      1    1/6      0    1/4    1/6  -1/12    5/6    1/6]
[ 12/31 -25/62   4/31  -9/62      1 -29/62  -9/62  10/31  17/62  -9/62],

[    1     0     0     0]
[    1     1     0     0]
[ -7/6  -3/4     1     0]
[  1/6   1/2 -4/31     1]
)
sage: M.gram_schmidt()[0] # rows are orthogonal, not orthonormal
[     1      0      0      1      0      0      1      1     -1      1]
[    -1      1      0      0      0     -1      0      1      0      0]
[  5/12    3/4      1    1/6      0    1/4    1/6  -1/12    5/6    1/6]
[ 12/31 -25/62   4/31  -9/62      1 -29/62  -9/62  10/31  17/62  -9/62]
sage: M.change_ring(RDF).gram_schmidt()[0] # orthonormal
[  0.408248290464              0.0              0.0   0.408248290464              0.0              0.0   0.408248290464   0.408248290464  -0.408248290464   0.408248290464]
[            -0.5              0.5              0.0              0.0              0.0             -0.5              0.0              0.5              0.0              0.0]
[  0.259237923683   0.466628262629   0.622171016838   0.103695169473              0.0    0.15554275421   0.103695169473 -0.0518475847365   0.518475847365   0.103695169473]
[  0.289303646409   -0.30135796501  0.0964345488031  -0.108488867403   0.747367753224  -0.349575239411  -0.108488867403   0.241086372008   0.204923416206  -0.108488867403]

矩阵st1 具有整数项，因此 Sage 将其视为整数矩阵，并尝试尽可能多地使用整数算术，如果失败，则使用有理算术。因此，Gram-Schmidt 正交归一化将失败，因为它涉及取平方根。这就是change_ring(RDF) 方法存在的原因：RDF 代表 Real Double Field。您可以改为将st1 的一个条目从1 更改为1.0，然后它将从一开始就将st1 视为RDF 上的矩阵，您无需在任何地方执行change_ring。

【讨论】：

Doggone 它，约翰，我只是在中间 ;-) 我将在下面添加一个澄清。
非常感谢。这些解释对我帮助很大。特别是关于有理数和整数算术的警告。我会尝试更频繁地使用它。
您也可以考虑将 Sage 问题发布到 ask.sagemath.org。我认为阅读该网站的 Sage 用户和开发人员比这个网站要多。
嗨，约翰。感谢您的补充评论。我以后会这样做的。

【解决方案3】：

为了扩展 John 的出色答案，我认为您对于同一个向量空间只有两个不同的基数。注意他使用right_kernel。

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]])
sage: st2 = matrix([[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]])
sage: st2 = st2.transpose()
sage: st2
[ 1 -1  0  0  1  0  0  0  0  0]
[-1  1  0  0  0 -1  0  1  0  0]
[ 0  1  1  0  0  0  0  0  1  0]
[ 2  0  1  1  0  1  1  0  0  1]
sage: st1.right_kernel()
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]
sage: st2.row_space()
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]

你的空间是一样的，只是在 Sage 和 Matlab 中的基础不同。

【讨论】：

感谢您进一步向我扩展和解释这一点。不得不承认我对矩阵魔法不是太精通，所以没有考虑这种可能性。