如何减少Julia代码内存分配？优化不动点求解程序

我因Julia的速度优势开始使用它，目前在求解一个不动点问题。现有代码运行速度尚可，但希望进一步优化性能。

算法概述 #

• 初始映射C0将空间(b,y)映射到行动空间c，即C0(b,y)
• 通过公式从C0生成规则Ct
• 借助额外约束更新b（记为bt），得到规则Ct(bt,y)
• 插值转换网格从bt到原始b，更新C0为C1
• 迭代直至C1与C0的距离低于收敛阈值

定义的结构与函数 #

我定义了两个结构分别存储参数和决策规则：

struct Parm
    lC::Array{Float64, 2}    # Lower limit
    uC::Array{Float64, 2}    # Upper limit
    γ::Float64               # CRRA coefficient
    δ::Float64               # factor in the euler
    γ1::Float64              # 
    r1::Float64              # inverse of the gross interest rate
    yb1::Array{Float64, 2}   # y - b(t+1)
    P::Array{Float64, 2}     # Transpose of transition matrix
end

mutable struct Upd1
    pol::Array{Float64,2}    # policy function
    b::Array{Float64, 1}     # exogenous grid for interpolation
    dif::Float64             # updating difference
end

同时实现了核心逻辑函数：

function eulerm(x::Upd1,p::Parm)
    ct = p.δ *(x.pol.^(-p.γ)*p.P).^(-p.γ1); #Euler equation 
    bt = p.r1.*(ct .+ p.yb1);               #Endeogenous grid for bonds
    return ct,bt
end

function interp0!(bt::Array{Float64},ct::Array{Float64},x::Upd1, p::Parm)
    polold = x.pol;
    polnew = similar(x.pol); 
    @inbounds @simd for col in 1:size(bt,2)
        F1 = LinearInterpolation(bt[:,col], ct[:,col],extrapolation_bc=Line());
        polnew[:,col] = F1(x.b);
    end
    polnew[polnew .< p.lC] .= p.lC[polnew .< p.lC];
    polnew[polnew .> p.uC] .= p.uC[polnew .> p.uC];
    dif = maximum(abs.(polnew - polold));
    return polnew,dif
end

function updating!(x::Upd1,p::Parm)
    ct, bt  = eulerm(x,p);       # endogeneous grid
    x.pol, x.dif   = interp0!(bt,ct,x,p);
end

function conver(x::Upd1,p::Parm)
    while x.dif>1e-8
        updating!(x,p);
    end
end

其中eulerm实现算法步骤2和3，updating!完成更新，conver负责迭代收敛。

核心问题与优化尝试 #

核心问题出在interp0!函数的循环中，该部分内存分配最多。尝试直接修改x.pol而非定义polnew会导致结果错误（仅两次迭代就收敛），推测是因为polold与x.pol指向同一对象，更新时会同时修改两者的值。

已对interp0!做如下修改，速度从约1.5秒提升至1.3秒，内存分配也有所减少：

function interp0!(bt::Array{Float64},ct::Array{Float64},x::Upd1, p::Parm)
    polold = x.pol;
    polnew = similar(x.pol); 
    @inbounds for col in 1:size(bt,2)
        F1 = LinearInterpolation(@view(bt[:,col]), @view(ct[:,col]),extrapolation_bc=Line());
        polnew[:,col] = F1(x.b);
    end
    for j in eachindex(polnew)
        polnew[j] < p.lC[j] ? polnew[j] =  p.lC[j] : nothing
        polnew[j] > p.uC[j] ? polnew[j] =  p.uC[j] : nothing
    end  
    dif = maximum(abs.(polnew - polold));
    return polnew,dif
end

遗留疑问 #

1. 将polnew[:,col] = F1(x.b)改为polnew[:,col] .= F1(x.b)会减少内存分配，但速度反而变慢，原因是什么？
2. 使用@time和@btime测试得到差异极大的结果：

   up1 = Upd1(c0,b,1);
   @time conver(up1,pp)
     1.338042 seconds (385.72 k allocations: 1.157 GiB, 3.37% gc time)
   up1 = Upd1(c0,b,1);
   @btime conver(up1,pp)
     4.200 ns (0 allocations: 0 bytes)
   

   是否意味着所有耗时都来自编译阶段的内存分配？

内容的提问来源于stack exchange，提问作者Christian Fernando Velasquez J