在 Julia 中用稀疏向量更新密集向量很慢
Updating a dense vector by a sparse vector in Julia is slow
我使用的是 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]这个操作我感觉应该有专门的函数做,但是没找到
感谢任何帮助。
假设您要计算 w += c * A[:, i]
,有一种简单的方法可以对其进行向量化:
>>> A = sprand(100000, 100000, 0.01)
>>> c = rand(100000)
>>> r1 = zeros(100000)
>>> @time for i = 1:100000
>>> r1 += A[:, i] * c[i]
>>> end
29.997412 seconds (1.90 M allocations: 152.077 GB, 12.73% gc time)
>>> @time r2 = sum(A .* c', 2);
1.191850 seconds (50 allocations: 1.493 GB, 0.14% gc time)
>>> all(r1 == r2)
true
首先,创建一个向量 c
要与之相乘的常量。然后按 c
的值对 A
的列进行多重播放(A .* c'
,它在内部进行广播)。最后,减少 A
的列(sum(.., 2)
部分)。
我在 GitHub 上问了同样的问题,我们得出了以下结论。从 Julia 0.4 开始添加了 SparseVector 类型,并添加了 BLAS 函数 LinAlg.axpy!,它通过稀疏向量 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)
这正是我需要实现的目标,经过一些研究我们设法实现了目标,谢谢大家!
我使用的是 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]这个操作我感觉应该有专门的函数做,但是没找到
感谢任何帮助。
假设您要计算 w += c * A[:, i]
,有一种简单的方法可以对其进行向量化:
>>> A = sprand(100000, 100000, 0.01)
>>> c = rand(100000)
>>> r1 = zeros(100000)
>>> @time for i = 1:100000
>>> r1 += A[:, i] * c[i]
>>> end
29.997412 seconds (1.90 M allocations: 152.077 GB, 12.73% gc time)
>>> @time r2 = sum(A .* c', 2);
1.191850 seconds (50 allocations: 1.493 GB, 0.14% gc time)
>>> all(r1 == r2)
true
首先,创建一个向量 c
要与之相乘的常量。然后按 c
的值对 A
的列进行多重播放(A .* c'
,它在内部进行广播)。最后,减少 A
的列(sum(.., 2)
部分)。
我在 GitHub 上问了同样的问题,我们得出了以下结论。从 Julia 0.4 开始添加了 SparseVector 类型,并添加了 BLAS 函数 LinAlg.axpy!,它通过稀疏向量 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)
这正是我需要实现的目标,经过一些研究我们设法实现了目标,谢谢大家!