8

私は神経画像からのデータを扱っていますが、大量のデータがあるため、コードに疎行列 (scipy.sparse.lil_matrix または csr_matrix) を使用したいと考えています。

特に、最小二乗問題を解くには、行列の疑似逆行列を計算する必要があります。メソッド sparse.lsqr を見つけましたが、あまり効率的ではありません。Moore-Penrose の疑似逆行列を計算する方法はありますか (正規行列の pinv に対応)。

私の行列 A のサイズは約 600'000x2000 で、行列のすべての行に 0 から最大 4 つのゼロ以外の値があります。マトリックス A のサイズは、ボクセル x 繊維束 (白質繊維トラクト) によって与えられ、ボクセル内で交差する最大 4 つのトラクトが予想されます。ほとんどの白質ボクセルでは、少なくとも 1 つのトラクトがあると予想されますが、ラインの約 20% はゼロになる可能性があると言えます。

ベクトル b はスパースであってはなりません。実際には、b には各ボクセルの測定値が含まれており、これは一般にゼロではありません。

エラーを最小限に抑える必要がありますが、ベクトル x にはいくつかの条件もあります。より小さな行列でモデルを試したので、これらの条件を満たすためにシステムを制約する必要はありませんでした (一般に 0

それは何かの助けになりますか？A の疑似逆を避ける方法はありますか?

ありがとう

6 月 1 日更新: ご協力ありがとうございます。私のデータについては何もお見せできません。なぜなら、Python のコードがいくつかの問題を引き起こしているからです。ただし、適切な k を選択する方法を理解するために、Matlab でテスト関数を作成しようとしました。

コードは次のとおりです。

F=zeros(100000,1000);

for k=1:150000
    p=rand(1);
    a=0;
    b=0;
    while a<=0 || b<=0
    a=random('Binomial',100000,p);
    b=random('Binomial',1000,p);
    end
    F(a,b)=rand(1);
end

solution=repmat([0.5,0.5,0.8,0.7,0.9,0.4,0.7,0.7,0.9,0.6],1,100);
size(solution)
solution=solution';
measure=F*solution;
%check=pinvF*measure;
k=250;
F=sparse(F);
[U,S,V]=svds(F,k);
s=svds(F,k);
plot(s)
max(max(U*S*V'-F))
for s=1:k
    if S(s,s)~=0
        S(s,s)=1/S(s,s);
    end
end

inv=V*S'*U';
inv*measure
max(inv*measure-solution)

F のサイズと比較して、k はどうあるべきか、何か考えはありますか? 私は 250 (1000 以上) を取得しましたが、結果は満足のいくものではありません (待ち時間は許容範囲内ですが、短くはありません)。また、結果を既知の解と比較できるようになりましたが、一般的に k を選択するにはどうすればよいでしょうか? また、取得した 250 個の単一値のプロットと、それらの二乗を正規化したものも添付しました。matlabでスクリプロットをより適切に行う方法が正確にはわかりません。私は今、より大きな k を使用して、突然値がはるかに小さくなるかどうかを確認しています。

ありがとう、ジェニファー

画像は、計算された 250 を示しています。 Matlab でスクリー プロットを作成する方法が正確にはわかりません。 二乗正規化単一値

4

2 に答える 2

9

scipy.sparse.linalgで提供されている代替案についてさらに学ぶことができます。

とにかく、スパース行列の疑似逆行列は（非常に）密なものである可能性が高いことに注意してください。したがって、スパース線形システムを解く場合、（一般的に）実際に実りある方法ではありません。

あなたはあなたの特定の問題をもう少し詳細に説明したいかもしれません（dot(A, x)= b+ e）。少なくとも以下を指定します。

  • 「典型的な」サイズのA
  • の非ゼロエントリの「典型的な」パーセンテージA
  • 最小二乗法はそれnorm(e)​​が最小化されていることを意味しますが、あなたの主な関心が上にあるx_hatか上b_hatにあるかe= b- b_hatを示してください。b_hat= dot(A, x_hat)

更新：ランクA（および列数よりはるかに小さい）についてある程度の知識がある場合は、合計最小二乗法を試すことができます。これは単純な実装です。ここで、kは使用する最初の特異値とベクトルの数です（つまり、「有効な」ランク）。

from scipy.sparse import hstack
from scipy.sparse.linalg import svds

def tls(A, b, k= 6):
    """A tls solution of Ax= b, for sparse A."""
    u, s, v= svds(hstack([A, b]), k)
    return v[-1, :-1]/ -v[-1, -1]
于 2011-05-04T08:05:25.810 に答える
7

私のコメントに対する答えに関係なく、ムーア・ペンローズSVD表現を使用してこれをかなり簡単に達成できると思います。SVDをscipy.sparse.linalg.svdsで検索し、Sigmaをその疑似逆行列で置き換えてから、V * Sigma_pi * U'を乗算して、元の行列の疑似逆行列を検索します。

于 2011-05-04T08:06:17.633 に答える