在Julia中通过稀疏向量更新密集向量很慢

Question

我正在使用Julia版本0.4.5并且我遇到了以下问题:据我所知,在稀疏矢量和密集矢量之间获取内积应该与通过稀疏矢量更新密集矢量一样快.后者慢得多.

A = sprand(100000,100000,0.01)
w = rand(100000)

@time for i=1:100000
  w += A[:,i]
end
26.304380 seconds (1.30 M allocations: 150.556 GB, 8.16% gc time)

@time for i=1:100000
  A[:,i]'*w
end
0.815443 seconds (921.91 k allocations: 1.540 GB, 5.58% gc time)

我创建了一个我自己的简单稀疏矩阵类型,加法代码与内积相同.

难道我做错了什么？我觉得应该有一个特殊的功能来做操作w + = A [:,i],但我找不到它.

任何帮助表示赞赏.

Answer 1

我在 GitHub 上问了同样的问题，我们得出了以下结论。从 Julia 0.4 开始添加了 SparseVector 类型，并添加了 BLAS 函数 LinAlg.axpy!，该函数x通过稀疏向量乘以y标量就地更新（可能是稠密）向量a，即执行x += a*y高效。然而，在 Julia 0.4 中它没有正确实现。它仅适用于 Julia 0.5

@time for i=1:100000
  LinAlg.axpy!(1,A[:,i],w)
end
1.041587 seconds (799.49 k allocations: 1.530 GB, 8.01% gc time)

然而，该代码仍然不是最优的，因为它创建了 SparseVector A[:,i]。通过以下功能可以获得更快的版本：

function upd!(w,A,i,c)
  rowval = A.rowval
  nzval = A.nzval
  @inbounds for j = nzrange(A,i)
    w[rowval[j]] += c* nzval[j]
  end
  return w
end

@time for i=1:100000
  upd!(w,A,i,1)
end
0.500323 seconds (99.49 k allocations: 1.518 MB)

这正是我需要实现的目标，经过一番研究我们成功实现了，谢谢大家！