Python中稀疏矩阵的伪逆

Question

我正在处理神经影像的数据，由于数据量很大，我想在我的代码中使用稀疏矩阵（比如 scipy.sparse.lil_matrix 或 csr_matrix）。

具体来说，我需要计算我的矩阵的伪逆，以解决一个最小二乘问题。我找到了一种方法 sparse.lsqr，但效率不是很好。有没有什么方法可以计算摩尔-彭若斯伪逆（对应于普通矩阵的 pinv）呢？

我的矩阵 A 的大小大约是 600,000 x 2,000，每一行最多有 4 个非零值。矩阵 A 的大小是由体素（voxel）和纤维束（白质纤维束）决定的，我们预计在一个体素中最多会有 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)

你对 k 应该与 F 的大小相比有什么想法吗？我取了 250（超过 1000），结果不太满意（等待时间可以接受，但不短）。现在我可以将结果与已知解进行比较，但一般来说，如何选择 k 呢？我还附上了我得到的 250 个单值及其归一化平方的图。我不太清楚如何在 Matlab 中更好地制作 scree 图。我现在正在尝试更大的 k，看看值是否会突然变得更小。

再次感谢，
Jennifer

Answer 1

不管我的评论有什么答案，我觉得你可以很简单地用Moore-Penrose的奇异值分解（SVD）来实现这个。首先用scipy.sparse.linalg.svds找到SVD，然后把Sigma换成它的伪逆，最后用V*Sigma_pi*U'来计算你原始矩阵的伪逆。

Answer 2

你可以多了解一下在 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]