julia 中内建函数与嵌套循环的效率

Introduction

前面做宽动态成像(前面的文章)时遇到了运行瓶颈，主要是在计算量上，要计算的式子如下：

[egin{eqnarray} I_{out}(x,y)=sum_{i=1}^{n}sum_{j=1}^{m}B_{ij}(x,y)U(x-rx_{ij},y-ry_{ij})I_{d_{ij}}(x,y) end{eqnarray} ]

而Bij(x,y)的计算也有求和运算：

[egin{eqnarray} B_{ij}(x,y)=frac{e^{- left( frac{ (x-rx_{ij})^2 }{2sigma^2_x} + frac{ (y-ry_{ij})^2 }{2sigma^2_y} ight) }} {sum_{p=1}^{m}sum_{q=1}^{n}e^{- left( frac{ (x-rx_{pq})^2 }{2sigma^2_x} + frac{ (y-ry_{pq})^2 }{2sigma^2_y} ight) } } end{eqnarray} ]

做实验的时候设置m=20,n=15,H,W分别为597,397，结果一循环就回不来了。。。（对julia的循环优化也是充满信心）

后面试了下用空间换效率的办法：增加一维存储m×n个的矩阵，然后用内置求和，发现提升明显。。。

...
BM=zeros(H,W,M*N)*0.
X=collect(1:W)'.-zeros(H)
Y=collect(1:H).-zeros(W)'

for pn=1:N
    for pm=1:M
        BM[:,:,pn+(pm-1)*N]=1./exp( (X-MNCenMatrix[pm,pn,1]).^2/(2*sigmaX^2) + (Y-MNCenMatrix[pm,pn,2]).^2/(2*sigmaY^2) ) .*
                                     convert(Array{Int64,2},epsilon.>map(max,abs(X-MNCenMatrix[pm,pn,1]),abs(Y-MNCenMatrix[pm,pn,2])))
    end
end
BM=BM./sum(BM,3)

for pn=1:N
    for pm=1:M
        retI[:]=retI[:]+(BM[:,:,pn+(pm-1)*N].*ILM[:,:,pn+(pm-1)*N])[:]
    end
end
...

Followup

issue 1: 与内建函数的差距

然后就想看看内建函数和嵌套函数差了多远：

H,W=900,800
C=300
m=rand(H,W,C)
n=rand(H,W)

function matrixSum!(m::Array{Float64,3},ret::Array{Float64,2})
H,W,C=size(m)
for pw=1:W
    for ph=1:H
        for pc=1:C
            ret[ph,pw]=ret[ph,pw]+m[ph,pw,pc]
        end
    end
end
end

@time n1=sum(m,3)
@time matrixSum!(m,n)

# i3 CPU 530 @ 2.93GHz × 4 
0.687958 seconds (19 allocations: 5.494 MB)
2.001088 seconds (6.04 k allocations: 246.260 KB)

对这结果甚是无语。。。分配了5.5MB也比你跑得快，是分配次数多限制了提升？难道这样的嵌套写法还可以大幅优化？（我一直以为这就是内存最简写法了）

issue 2: 内存访问顺序

突然想到访问次序有些关键，于是调整了下：

H,W=900,800
C=300
m=rand(H,W,C)
n=rand(H,W)
function matrixSum!(m::Array{Float64,3},ret::Array{Float64,2})
H,W,C=size(m)

for pc=1:C
    for pw=1:W
        for ph=1:H
            ret[ph,pw]=ret[ph,pw]+m[ph,pw,pc]
        end
    end
end
end

@time n1=sum(m,3)
@time matrixSum!(m,n)

# 首次运行
0.879469 seconds (193.28 k allocations: 13.598 MB)
0.787520 seconds (6.47 k allocations: 269.166 KB)

# 再次运行
0.609809 seconds (19 allocations: 5.494 MB)
0.788194 seconds (6.04 k allocations: 246.260 KB)

就内建函数来看，内存分配的次数并没有太多的决定时间消耗;
就内建函数与循环嵌套对比来看，刚才是访问顺序的问题。
那么内存为什么分配了那么多次呢？

issue 3: 内存分配次数

突然想到julia的动态编译机制，于是拿出来单独运行了次：

julia> @time matrixSum!(m,n)
  0.765424 seconds (4 allocations: 160 bytes)

看来是这样的。

Conclusion

julia 的优化确实做得不错，几乎达到了极致;
访问顺序很关键;
分配次数并不一定会是瓶颈。