2017-11-10 14 views
0

私はSageでプログラミングするのが初めてです。私は矩形のR * C行列(R行とC列)を持ち、MのランクはRとCの両方よりも(おそらく)小さいです。ターゲットベクトルTが列のサブセットのスパンにあるかどうかをチェックしたいM.私はSageに以下のコードを書いています(私はMとTのやり方がやや面倒なので、コード全体は含めていません)。私はちょうどコードが私が望むことをしているかどうかをチェックしたい。ベクトルがSageの行列の列のサブセットにまたがるかどうか確認してください

簡単に、これは私のコードがしようとしていることです:Mは私の与えられた行列です、私は最初にMの列のスパン(最初のif条件)にあることを確認します。そうであれば、M(C列を持つ)をM列の列(これは最初のwhileループと同じです)を持つ行列M1にトリムします。その後、列の1つずつを削除して残りの列に範囲がTを含むかどうかを確認します(これは2番目のwhileループです)。 2番目のwhileループでは、まずM2から基本的にM1のコピーである列を削除し、この行列M3を呼び出します。 M3に。ベクトルTを増やし、階数が減少するかどうかを調べます。 Tは既にM2の範囲にあるので、ランク([M2 T])はランク(M2)と同じでなければならない。今度は列cを削除してM2をTに増やしてもランクは下がらないので、cはTを生成する必要はないことがわかります。このようにして、Tを生成するために必要な列のみを保持します。

私が試した例の正解ですが、このコードを大規模に変化するエントリを持つ行列で実行しようとしています(たとえば、最大値は20^20、最小値は1です)、通常行列の次元はしたがって、週末に数百のテストケースを実行する予定です。何かが怪しい/間違っているかどうかを教えてもらえると本当に役に立ちます。私は精度のエラーに遭遇するでしょうか?上記のようにすべての値/範囲で動作するようにコードを変更するにはどうすればよいですか?また、私のコードを高速化する方法がある場合(または短い/より良い方法で同じものを書く場合)、私は知りたいと思います。

R = 155 
C= 167 
T = vector(QQ, R) 
M1 = matrix(ZZ, R, C) 

M1 = M 
C1 = C 
i2 = 0 

if rank(M.augment(T)) == rank(M): 
print("The rank of M is") 
print(rank(M)) 
while i2 < C1 : 
    if rank(M1.delete_columns([i2])) == rank(M1) : 
    M1 = M1.delete_columns([i2]) 
    C1 = C1 - 1 
    else : 
    i2 = i2+1 

C2 = M1.ncols() 
print("The number of columns in the trimmed down matrix M1 is") 
print(C2) 


i3 = 0 
M2 = M1 
print("The rank of M1 which is now also the rank of M2 is") 
print(rank(M2)) 
while i3 < C2 : 
    M3 = M2.delete_columns([i3]) 
    if rank(M3.augment(T)) < rank(M2) : 
    M2 = M3 
    C2 = C2 - 1 
    else : 
    i3 = i3 + 1 

print("Rank of matrix M is") 
print(M.rank()) 
+0

私はまた、ランクパッケージがセージでどのように機能するのか興味がある。 MatrixSpace(QQ、R、C)の代わりにMatrixSpace(RR、R、C)を使用しています。なぜこれが起こるのですか?私の目的のために、行列の正確な階数が必要であり、いくつかの許容誤差より大きい特異値の数に基づく数値近似ではありません(これは私がSageを使用するようにしたmatlabの考え方です)。私はシンボリックランク - MatrixSpace(SR、R、C)を使用しようとしましたが、コードは非常に遅くなりました。だから、これらの線の周りの任意のポインタ/ヒント? – Nikhil

+0

私は数値アナリストではありませんが、マトリックスのランク計算はマトリックスのエントリーの小さな変化に非常に敏感であると思われますので、RRのような不正確なフィールドのランク計算は疑わしいです。 QQのような、数値的な不確実性のない分野では、ランクがうまく機能します。 –

答えて

1

私はベクトルTは別の行列Mの列のサブセットから構築いくつかの行列M1のイメージにあったかどうかを決定するためにセージを使用したい場合は、私はこれだろう。

M1 = M.matrix_from_columns([list of indices of the columns to use]) 
T in M1.column_space() 

またはwhileループを使用して毎回M1を変更してください。 (しかし、私はT in M1.column_space()がランクの同等性をテストするよりもうまくいくはずだと思う)

関連する問題