python - 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

合理的基底ヌルスペース（零空間）は、次の方法で計算できます。

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

私はこれをPythonSAGEで複製しようと試みてきましたが、これまでのところ成功していません。私のコードは次のようになります。

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と同じ結果を取得するために私が何をする必要があるか、そして可能であれば、正規直交結果を取得する方法（K_nuoと呼ばれるMatlab）を知っている人はいますか？

チュートリアルやドキュメントなどを調べてみましたが、今のところ運がありません。

score 5 · Accepted Answer

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スペースプロパティがあることが確認されます。

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

Kこれは、それが正規直交であることを確認します。

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

score 3 · Accepted Answer

このようなものが機能するはずです：

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はそれを整数の行列として扱い、整数演算で可能な限り多くのことを実行しようとしますが、それに失敗すると、有理演算になります。このため、グラムシュミット正規化は平方根を取る必要があるため失敗します。これがメソッドchange_ring(RDF)が存在する理由です。RealDoubleFieldのRDF略です。代わりに、st1からの1つのエントリを変更するだけで、最初からRDF上のマトリックスとして扱われ、どこでもこれを行う必要はありません。11.0st1change_ring

score 3 · Accepted Answer

ジョンの素晴らしい答えを拡張するために、同じベクトル空間に対して2つの異なるベースがあると思います。彼の使用に注意してください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の拠点が異なります。

python - Python sage：化学量論的行列のヌルスペース（零空間）を計算するにはどうすればよいですか？

3 に答える 3

Related

Reference