首页 文章

python中稀疏矩阵的伪逆

提问于
浏览 次
6

我正在使用来自神经成像的数据，并且由于数据量很大，我想为我的代码使用稀疏矩阵（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的伪逆？

谢谢

Update 1st June: 再次感谢您的帮助 . 我可以't really show you anything about my data, because the code in python give me some problems. However, in order to understand how I could choose a good k I'试图在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的大小相比应该是什么？我已经花了250（超过1000）并且结果不令人满意（等待时间是可以接受的，但不是很短） . 现在我也可以将结果与已知解决方案进行比较，但是如何选择k呢？我还附上了我获得的250个单值的图，并将它们的正方形归一化 . 我不知道如何更好地在matlab中做一个screeplot . 我现在正在进行更大的k，以确定是否会突然显示该值会小得多 .

再次感谢，珍妮弗

The image shows the 250 computed. I don't know exactly how to create a scree plot in Matlab.

squared normalized single values

python linear-algebra sparse-matrix least-squares matrix-inverse

2 回答

  • 4

    您可以更多地了解scipy.sparse.linalg中提供的替代方案 .

    无论如何，请注意，稀疏矩阵的伪逆最有可能是（非常）密集的，因此在求解稀疏线性系统时，它并不是一个非常富有成效的途径（通常） .

    您可能希望稍微更详细地描述您的特定问题（ dot(A, x)= b+ e ） . 至少指定：

    • 'typical'大小 A

    • 'typical' A 中非零条目的百分比

    • 最小二乘意味着 norm(e) 被最小化，但请指出您的主要兴趣是 x_hat 还是 b_hat ，其中 e= b- b_hat 和 b_hat= dot(A, x_hat)

    Update ：如果您对 A 的等级有了一些了解（并且它远小于列数），您可以尝试使用total least squares方法 . 这是一个简单的实现，其中 k 是要使用的第一个奇异值和向量的数量（即'effective' rank） .

    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]
    回复于
  • 6

    无论我的评论的答案如何，我认为你可以使用Moore-Penrose SVD representation轻松地完成这个 . 使用scipy.sparse.linalg.svds找到SVD，用pseudoinverse替换Sigma，然后乘以V * Sigma_pi * U'找到原始矩阵的伪逆 .

    回复于

相关问题