Julia加速递归系统

我需要加快这个递归系统。谁能帮忙？Julia加速递归系统

我注释了所有变量，并且只使用了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）非常快。

任何帮助，将不胜感激。

来源

2015-11-12 amrods

它看起来像一个无限的，你会得到正确的结果吗？ –

你也可以看看https://github.com/one-more-minute/Lazy中的'@ lazy'宏。jl – SalchiPapa

我觉得代码就像一个怪物！
你已经有3个递归函数U(), V(), Z()后，U()电话V()，反之亦然，那么V()电话Z()本身调用U() ......
尝试重写它并防止递归调用尽可能多的，这会结束于更有效的解决方案。 检查这个 - >recursion versus iteration
之后，你必须U()和V()内for循环，必须在外部进行不必要的重复呼叫。

第一步：临时变量

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数据类型的商场使用，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

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

来源

2015-11-13 05:35:17

到目前为止，我一直认为迭代比递归更好，但现在我怀疑了，因此在最近的编辑中，我删除了一个建议使用迭代而不是递归方法的语句。任何人都可以针对上述问题编写更高效的纯迭代方法吗？ –

非常感谢@reza。在它的第一个版本中，我有几个变量作为整数，我将它们全部切换到Float64，因为我读到了Julia在不需要在类型之间进行转换时速度更快，但我不是计算机人员和我用于过去没有明确的类型（我正在谈论Mathematica和其他统计软件包）。我将检查您对我的代码所做的编辑。只是为了记录，系统必须有一个解决方案，因为我有t和年龄的终端条件。这只是一个反向递归系统。 – amrods

这样的问题可以并行吗？ – amrods

Julia加速递归系统

回答

相关问题