Python sage: How do I computer a nullspace (kernel) for a stoichiometric matrix?

Question

In a desperate attempt to switch from Matlab to python, I am encountering the following problem:

In Matlab, I am able to define a matrix like:

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

The rational basis nullspace (kernel) can then be calculated by:

K_nur= null(N,'r')

And the orthonormal basis like:

K_nuo= null(N)

This outputs the following:

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

I have been trying to replicate this in Python SAGE, but so far, I have had no success. My code looks like this:

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

Note: The transpose was introduced after reading through some tutorials on this and has to do with the nature of SAGE automatically computing the kernel from the left (which in this case yields no result at all).

The problem with this now is: It DOES print me something... But not the right thing.

The output is as follows:

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

Upon closer inspection, you can see that the resulting kernel matrix (nullspace) looks similar, but is not the same.

Does anyone know what I need to do to get the same result as in Matlab and, if possible, how to obtain the orthonormal result (in Matlab called K_nuo).

I have tried to look through the tutorials, documentation etc., but so far, no luck.

Answer 1

Something like this should work:

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

The gram_schmidt method gives a pair of matrices. Type M.gram_schmidt? to see the documentation.

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]

The matrix st1 has integer entries, so Sage treats it as a matrix of integers, and tries to do as much as possible with integer arithmetic, and failing that, rational arithmetic. Because of this, Gram-Schmidt orthonormalization will fail, since it involves taking square roots. This is why the method change_ring(RDF) is there: RDF stands for Real Double Field. You could instead just change one entry of st1 from 1 to 1.0, and then it will treat st1 as a matrix over RDF from the start and you won't need to do this change_ring anywhere.

Answer 2

There might be a way do this with SAGE builtin functions; I'm not sure.

However, if a numpy/python-based solution will do, then:

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)

yields the orthonormal null space:

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

this confirms the columns have the null space property:

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

and this confirms that K is orthonormal:

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