在 Julia 中加速递归系统

Question

我需要加速这个递归系统。有人可以帮忙吗？

我注释了所有变量并且只使用了 Float64 一致的函数。我还尝试通过 Memoize.jl 使用 memoization（因为系统是递归的），但我必须在优化算法中计算 U、V 和 Z 数千次，并且我的计算机很快就会耗尽内存。

const Tage = 60.0
const initT = 1975.0
const Ttime = 1990.0
const K = 6.0
const β = 0.95
const s = 0.5

θ0 = [
    0.0011
 -0.0045
 -0.0075
 -0.0013
  3.60   
  0.2587
 -0.0026  
 -0.0528
  0.0060
  1.3772
  0.2932
 -0.0030
  0.0021
  0.0044
 10.2593    
  0.7498
  1.1765
]./10000

function wagem(t::Float64, agem::Float64, edum::Float64, θ::Vector{Float64})
    θ[5] + θ[6]*agem + θ[7]*agem^2 + θ[8]*edum + θ[9]*edum^2
end

function wagef(t::Float64, agef::Float64, eduf::Float64, θ::Vector{Float64})
    θ[10] + θ[11]*agef + θ[12]*agef^2 + θ[13]*eduf + θ[14]*eduf^2
end

function ζ(agem::Float64, edum::Float64, agef::Float64, eduf::Float64, θ::Vector{Float64})
    θ[1]*agem*agef + θ[2]*agem*eduf + θ[3]*edum*agef + θ[4]*edum*eduf
end

function z(t::Float64, agem::Float64, edum::Float64, agef::Float64, eduf::Float64, θ::Vector{Float64})
    ζ(agem, edum, agef, eduf, θ) + wagem(t, agem, edum, θ) + wagef(t, agef, eduf, θ)
end

function dc(θ::Vector{Float64})
    ret::Float64 = 0.0
    ret = θ[15]
    return ret
end

function Z(t::Float64, agem::Float64, edum::Float64, agef::Float64, eduf::Float64, θ::Vector{Float64})
    ret::Float64 = 0.0
    if agem >= Tage || agef >= Tage # terminal age conditions
        ret = 0.0
    elseif t >= Ttime # terminal time condition
        ret = 1.5
    else
        ret = log(
        exp(z(t, agem, edum, agef, eduf, θ) + β*Z(t+1, agem+1, edum, agef+1, eduf, θ))
        + exp(-dc(θ) + β*(V(t+1, agem+1, edum, θ) + U(t+1, agef+1, eduf, θ)))
        )
    end
    return ret
end

function V(t::Float64, agem::Float64, edum::Float64, θ::Vector{Float64})
    ret::Float64 = 0.0
    if agem >= Tage
        ret = 0.0
    elseif t >= Ttime
        ret = 1.5
    else
        suma::Float64 = 0.0
        for agef in 16.0:Tage-1, eduf in 1.0:K
            suma += exp(s*z(t, agem, edum, agef, eduf, θ) + wagem(t, agem, edum, θ) + β*(s*(Z(t+1, agem+1, edum, agef+1, eduf, θ) - V(t+1, agem+1, edum, θ) - U(t+1, agef+1, eduf, θ)) + V(t+1, agem+1, edum, θ)))
        end
        ret = log(
                exp(wagem(t, agem, edum, θ) + β*(V(t+1, agem+1, edum, θ))) + suma
                )
    end
    return ret
end

function U(t::Float64, agef::Float64, eduf::Float64, θ::Vector{Float64})
    ret::Float64 = 0.0
    if agef >= Tage
        ret = 0.0
    elseif t >= Ttime
        ret = 1.5
    else
        suma::Float64 = 0.0
        for agem in 16.0:Tage-1, edum in 1.0:K
            suma += exp((1-s)*z(t, agem, edum, agef, eduf, θ) + wagef(t, agef, eduf, θ) + β*((1-s)*(Z(t+1, agem+1, edum, agef+1, eduf, θ) - V(t+1, agem+1, edum, θ) - U(t+1, agef+1, eduf, θ)) + U(t+1, agef+1, eduf, θ)))
        end
        ret = log(
        exp(wagef(t, agef, eduf, θ) + β*(U(t+1, agef+1, eduf, θ))) + suma
        )
    end
    return ret
end

我希望能够非常快速地计算 U(1984.0, 17.0, 3.0, θ0) 或 V(1975.0, 16.0, 1.0, θ0) 之类的东西。

任何帮助将不胜感激。

Answer 1

我认为代码就像一个怪物！
之后你已经有 3 个递归函数 U(), V(), Z()，U() 调用 V() 反之亦然，然后 V() 调用 Z() 本身调用 U() ......
试试重写它并尽可能地防止递归调用，这将结束一个更有效的解决方案。 检查这个-> recursion versus iteration
之后，您必须在 for 循环内对 U() 和 V() 进行不必要的重复调用，而这些调用必须在外部完成。

第一步：临时变量

function U(t::Float64, agef::Float64, eduf::Float64, θ::Vector{Float64}) 
    ret::Float64 = 0.0
    if agef >= Tage
        ret = 0.0
    elseif t >= Ttime
        ret = 1.5
    else
        suma::Float64 = 0.0
        U1 = U(t+1, agef+1, eduf, θ) # nothing to do with following loop, so I take it outside
        wagef1 = wagef(t, agef, eduf, θ) # same as above comment
        for agem in 16.0:Tage-1, edum in 1.0:K
            suma += exp((1-s)*z(t, agem, edum, agef, eduf, θ) + wagef1 + β*((1-s)*(Z(t+1, agem+1, edum, agef+1, eduf, θ, ui) - V(t+1, agem+1, edum, θ) - U1) + U1))
        end
        ret = log(
                  exp(wagef1 + β*(U1)) + suma
                  )
    end

    return ret
end

julia> @time U(1984.0, 17.0, 3.0, t0) # a test with const Tage = 19.0
  0.869675 seconds (102.77 k allocations: 2.148 MB, 0.57% gc time)
3.785563216949393

U() and V() 中的上述编辑使我的运行速度比原始代码快 28 倍。

第二步：使用正确的类型

接下来是 Float64 数据类型的 mall-usage，IMO 只有 θ,β,s 需要为 Float64 类型，所有其他变量都应声明为整数。

julia> @time U(1984, 17, 3, t0) # a test with const Tage = 19
  0.608982 seconds (125.77 k allocations: 2.530 MB, 0.89% gc time)
3.785563216949393

现在速度提高了 30%。 与之前的代码相比

第三步 Memoize.jl

although already some improvements have been happened for bigger problems they are not enough, but with using the right data types, in the previous step now `Memoize.jl` will be usable:  


using Memoize
const Tage = 60%Int
const initT = 1975%Int
const Ttime = 1990%Int
const K = 6%Int
const β = 0.95
const s = 0.5

t0 = [
    0.0011
 -0.0045
 -0.0075
 -0.0013
  3.60   
  0.2587
 -0.0026  
 -0.0528
  0.0060
  1.3772
  0.2932
 -0.0030
  0.0021
  0.0044
 10.2593    
  0.7498
  1.1765
]./10000

function wagem(t::Int, agem::Int, edum::Int, θ::Vector{Float64})
    θ[5] + θ[6]*agem + θ[7]*agem^2 + θ[8]*edum + θ[9]*edum^2
end

function wagef(t::Int, agef::Int, eduf::Int, θ::Vector{Float64})
    θ[10] + θ[11]*agef + θ[12]*agef^2 + θ[13]*eduf + θ[14]*eduf^2
end

function ζ(agem::Int, edum::Int, agef::Int, eduf::Int, θ::Vector{Float64})
    θ[1]*agem*agef + θ[2]*agem*eduf + θ[3]*edum*agef + θ[4]*edum*eduf
end

function z(t::Int, agem::Int, edum::Int, agef::Int, eduf::Int, θ::Vector{Float64})
    ζ(agem, edum, agef, eduf, θ) + wagem(t, agem, edum, θ) + wagef(t, agef, eduf, θ)
end

function dc(θ::Vector{Float64})
    ret::Float64 = 0.0
    ret = θ[15]
    return ret
end

function Z(t::Int, agem::Int, edum::Int, agef::Int, eduf::Int, θ::Vector{Float64})
    ret::Float64 = 0.0
    if agem >= Tage || agef >= Tage # terminal age conditions
        ret = 0.0
    elseif t >= Ttime # terminal time condition
        ret = 1.5
    else
        ret = log(
        exp(z(t, agem, edum, agef, eduf, θ) + β*Z(t+1, agem+1, edum, agef+1, eduf, θ))
        + exp(-dc(θ) + β*(V(t+1, agem+1, edum, θ) + U(t+1, agef+1, eduf, θ)))
        )
    end
    return ret
end

@memoize function V(t::Int, agem::Int, edum::Int, θ::Vector{Float64})
    ret::Float64 = 0.0
    agef::Int = 0
    eduf::Int = 0
    if agem >= Tage
        ret = 0.0
    elseif t >= Ttime
        ret = 1.5
    else
        suma::Float64 = 0.0
        V1 = V(t+1, agem+1, edum, θ)
        wagem1 = wagem(t, agem, edum, θ)
        for agef in 16:Tage-1, eduf in 1:K
            suma += exp(s*z(t, agem, edum, agef, eduf, θ) + wagem1 + β*(s*(Z(t+1, agem+1, edum, agef+1, eduf, θ) - V1 - U(t+1, agef+1, eduf, θ)) + V1))
        end
        ret = log(
                exp(wagem1 + β*(V1)) + suma
                )
    end
    return ret
end

@memoize  function U(t::Int, agef::Int, eduf::Int, θ::Vector{Float64})
    ret::Float64 = 0.0
    agem::Int = 0
    edum::Int = 0
    if agef >= Tage
        ret = 0.0
    elseif t >= Ttime
        ret = 1.5
    else
        suma::Float64 = 0.0
        U1 = U(t+1, agef+1, eduf, θ)
        wagef1 = wagef(t, agef, eduf, θ)
        for agem in 16:Tage-1, edum in 1:K
            suma += exp((1-s)*z(t, agem, edum, agef, eduf, θ) + wagef1 + β*((1-s)*(Z(t+1, agem+1, edum, agef+1, eduf, θ) - V(t+1, agem+1, edum, θ) - U1) + U1))
        end
        ret = log(
                  exp(wagef1 + β*(U1)) + suma
                  )
    end  
    return ret
end

现在Tage = 60：

julia> @time U(1984, 17, 3, t0) # a test with const Tage = 60
  3.789108 seconds (27.08 M allocations: 494.135 MB, 15.12% gc time)
16.99266161490023

我也测试了Int16 版本：

julia> @time U(1984%Int16, 17%Int16, 3%Int16, t0)
  3.596871 seconds (27.05 M allocations: 413.808 MB, 11.93% gc time)
16.99266161490023

与Int32 版本相比，它的运行速度提高了 5%，内存分配减少了 20%