Python中稀疏矩阵的伪逆

Question

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

特别是，我需要计算矩阵的伪逆来解决最小二乘问题。我找到了sparse.lsqr方法，但效率不高。有没有一种方法可以计算Moore-Penrose伪逆（对应于普通矩阵的pinv）。

我的矩阵A的大小约为600,000 x 2000，在矩阵的每一行中，我将有从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)

你有没有想过k应该与F的大小相比如何？我已经尝试了250（超过1000），但结果并不令人满意（等待时间可接受，但不短）。 现在我可以将结果与已知解进行比较，但通常应该如何选择k呢？ 我还附上了250个单值及其标准化平方的绘图。我不知道如何在matlab中更好地绘制screeplot。我现在正在尝试更大的k值，看看值是否会突然变小。

再次感谢， Jennifer

- Jennifer

---

我对sparse.linalg中lsqr的实现一无所知，但似乎很难想象有更快的计算伪逆的方法；如果有的话，lsqr不就直接调用那个方法然后进行乘法运算吗？编辑：或者问题是你需要做许多这样的计算，并且希望使用伪逆来避免重新计算？ - bnaul

Answer 1

您可以在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]

- eat

---

谢谢您的回答，但在这两种情况下，我都不明白如何处理SVDS中k必须小于秩的事实。 - Jennifer

@Jennifer：首先，实际原因决定了k必须相对较小，请注意u和v可能是密集的。 - eat

@Jennifer：其次，有一些理论原因，考虑到sum(s)（完整的svd）解释了[A b]´中的所有变化，因此您应该找到一个k，其中sum(s[:k])将解释[A b]´中合理数量的变化（考虑这个信号），而sum(s[k:])将解释其余部分（考虑这个噪音）。您可能希望更详细地研究svd背后的线性代数，但在实践中，您可以为一些可行的k绘制s[:k]，如果您能够在那里识别出一个变化点，则将其用作您的k，有关更多详细信息，请搜索“scree plot”。 - eat

@Jennifer：你能给我们展示一下针对你能够在合理时间内计算出的这个k值的scree plot吗？谢谢。 - eat

Answer 2

无论对我的评论做出何种回答，我认为您可以使用Moore-Penrose SVD表示法相当容易地完成此操作。使用scipy.sparse.linalg.svds找到SVD，将Sigma替换为其伪逆，然后乘以V*Sigma_pi*U'，以找到原始矩阵的伪逆。

- bnaul

---

@eat 你说得对，我把 V 和 U 搞反了。我认为其余部分是正确的，但我同意你的答案观点，伪逆可能不是解决她问题的最佳方法。 - bnaul