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)

你知道 k 与 F 的大小相比应该是多少吗？我已经拿了250（超过1000），结果并不理想（等待时间可以接受，但不短）。现在我也可以将结果与已知解决方案进行比较，但一般来说如何选择 k 呢？我还附上了我得到的 250 个单个值的图，并将它们的平方归一化。我不知道如何更好地在 matlab 中做一个 screeplot。我现在正在处理更大的 k，看看值是否会突然变小。

再次感谢，詹妮弗

该图像显示了计算的 250。 我不确切知道如何在 Matlab 中创建碎石图。 平方归一化单个值

4

2 回答 2

9

您可以更多地研究scipy.sparse.linalg提供的替代方案。

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

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

  • “典型”大小A
  • 中非零条目的“典型”百分比A
  • 最小二乘法意味着norm(e)最小化，但请说明您的主要兴趣是 onx_hat还是 on 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

不管我的评论的答案如何，我认为您可以使用Moore-Penrose SVD 表示相当轻松地完成此任务。使用 scipy.sparse.linalg.svds 找到 SVD，将 Sigma 替换为其伪逆，然后将 V*Sigma_pi*U' 相乘以找到原始矩阵的伪逆。

于 2011-05-04T08:06:17.633 回答