4

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

我试图浏览教程、文档等，但到目前为止，还没有运气。

4

3 回答 3

5

使用 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
于 2012-11-01T11:23:49.277 回答
3

为了扩展约翰的出色答案，我认为您对于同一个向量空间只有两个不同的基数。注意他使用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 中的基础不同。

于 2012-11-01T16:22:38.060 回答
3

像这样的东西应该工作：

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代表真正的双场。您可以改为将st1from的一个条目更改1为1.0，然后它将st1从一开始就将其视为 RDF 上的矩阵，您无需在change_ring任何地方执行此操作。

于 2012-11-01T15:58:28.757 回答