gpt4 book ai didi

julia - Julia 中高效的逐元素矩阵运算

转载 作者:行者123 更新时间:2023-12-04 20:33:22 26 4
gpt4 key购买 nike

我需要执行(复杂)矩阵的离散卷积,并在 Julia 中定义了以下函数:

function convolve(M::Array{Complex{Float64},2}, K::Array{Float64,2}, p::Int)
(n,m) = size(M)
res = zeros(Complex{Float64},n)
for k=1:p
for l=1:n
res[l] += M[l,k]*K[l,end-p+k]
end
end
return res
end

我像这样使用它:
M=complex(rand(2000,2000))
K=rand(2000,2000)
@time convolve(M,K,2000,0)

现在,这比我用 res += M[:,k].*K[:,end-p+k] 替换内循环的矢量化版本要快(大约 3 倍)。 . (我认为这是由于临时数组的大量内存分配,我可以忍受)。

但是矢量化的 MATLAB 代码运行速度大约快 5 倍:
function res = convolve(M, K, p)
n = size(M,1);
res = zeros(n,1);
for k=1:p
res = res + M(:,k).*K(:,end-p+k);
end
end

我做错了什么,我如何让 Julia 像 MATLAB 一样快地执行这种逐元素乘法?是索引问题吗?

注意:我已经检查过 @code_warntype没有类型优柔寡断的有趣业务(没有 AnyUnion 等),但问题可能更微妙。宏 @code_llvm产生了惊人的长输出,但我不是专家,所以我很难看到发生了什么。

最佳答案

以下版本在我的机器上速度更快:

function convolve2(M::Array{Complex{Float64},2}, K::Array{Float64,2}, p::Int)
(n,m) = size(M)
res = zeros(Complex{Float64},n)
offset = size(K,2)-p
(p>m || offset<0) && error("parameter p ($p) out of bounds")
@inbounds for k=1:p
@inbounds @simd for l=1:n
res[l] += M[l,k]*K[l,offset+k]
end
end
return res
end

请注意 @simd目前在许多 CPU 中使用向量指令的加法。

编辑:OP 代码中的性能下降似乎源于使用 endK 的索引中在热循环线上。重新定义 Base.trailingsize@inline使 LLVM 内联 end (在我的机器上)并使两个版本的运行速度大致相同。使用的代码:
import Base: trailingsize
@inline function Base.trailingsize(A, n)
s = 1
for i=n:ndims(A)
s *= size(A,i)
end
return s
end

请参阅对该问题的评论 #19389关于这一点。

关于julia - Julia 中高效的逐元素矩阵运算,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/40744675/

26 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com