一段简单代码的优化问题

综合讨论区 性能优化
1

各位好，我刚从Matlab转用Julia。在实现一段简单代码的优化时遇到了一些问题，想请大家帮忙看一下。

首先我把Matlab代码“直译”为Julia代码，如下

using LinearAlgebra
using SpecialFunctions
using BenchmarkTools

function main(θ::Array{Float64,1}, ε::Array{Float64,1}, t::Array{Float64,1}, h::Float64)
    α = θ[1]
    b = θ[2]
    σ = zeros(size(ε))
    ω = cumprod(vcat(1, (1 .- (α.+1)./(1:length(ε)-1))))
    
    for i = 1:length(t)
        w_t = ω[1:i]
        σ = σ .- σ[1]
        y_t = σ[i:-1:1]
        σ[i] = h^(-α)*dot(w_t, y_t)
    end
    σ = @. b*(σ + ε*t^(-α) / gamma(1-α))
    return σ
end

α = 0.3
h = 0.0001
t = collect(0:h:2)
θ = [α, 5.0]
ε = 1.0 .* ones(size(t))

@btime main(θ, ε, t, h)`

结果为

3.983 s (115938 allocations: 5.97 GiB)

在浏览了社区内一些关于Julia性能优化的帖子之后，为避免使用全局变量，我将代码写成如下形式

using LinearAlgebra
using SpecialFunctions
using BenchmarkTools

mutable struct Params
    #
    α :: Float64
    b :: Float64
    #
    ω :: Array{Float64,1}
    #
    h :: Float64
    t :: Array{Float64,1}
    ε :: Array{Float64,1}
    σ :: Array{Float64,1}
end

"""
Initial parameters.
"""
function initialparameters(params::Array{Float64,1}, ε::Array{Float64,1}, t::Array{Float64,1}, h::Float64)
    α = params[1]
    b = params[2]
    σ = zeros(size(ε))
    ω = cumprod(vcat(1, (1 .- (α.+1)./(1:length(ε)-1))))
    Params(α, b, ω, h, t, ε, σ)
end

"""
Kernel function for computing σ(t = tn)
"""
function computesigma!(p::Params, i::Int64)
    w_t = p.ω[1:i]
    p.σ = p.σ .- p.σ[1]
    y_t = p.σ[i:-1:1]
    p.σ[i] = p.h^(-p.α)*dot(w_t, y_t)
end

"""
Main function
"""
function main(θ::Array{Float64,1}, ε::Array{Float64,1}, t::Array{Float64,1}, h::Float64)
    p = initialparameters(θ, ε, t, h)
    for i in 1:length(p.t)
        computesigma!(p, i)
    end
    p.σ = @. p.b*(p.σ + p.ε*p.t^(-p.α) / gamma(1-p.α))
    p.σ
end

α = 0.3
h = 0.0001
t = collect(0:h:2)
θ = [α, 5.0]
ε = 1.0 .* ones(size(t))

@btime main(θ, ε, t, h)

结果为

3.960 s (115937 allocations: 5.97 GiB)

与之前的相比没有区别。

请大家帮忙看下第二段代码的写法有什么问题？如果有其他优化建议，更不胜感激。

2

有一些内存是可以重复使用的，比如说：

- p.σ = @. p.b*(p.σ + p.ε*p.t^(-p.α) / gamma(1-p.α))
+ @. p.σ = p.b*(p.σ + p.ε*p.t^(-p.α) / gamma(1-p.α))
3

对数组切片用@view好些吧，但是我这里不知道为啥段错误了……

4 
julia> function main(θ::Array{Float64,1}, ε::Array{Float64,1}, t::Array{Float64,1}, h::Float64)
           α = θ[1]
           b = θ[2]
           σ = zeros(size(ε))
           ω = cumprod(vcat(1, (1 .- (α.+1)./(1:length(ε)-1))))
           
           @inbounds for i = 1:length(t)
               w_t = @view ω[1:i]
               σ .-= σ[1]
               y_t = @view σ[i:-1:1]
               σ[i] = h^(-α)*dot(w_t, y_t)
           end
           @. σ = b*(σ + ε*t^(-α) / gamma(1-α))
           return σ
       end
 main (generic function with 1 method)

julia> @btime main(θ, ε, t, h);
  492.309 ms (25 allocations: 626.14 KiB)
5

多谢！但是我试着跑了几次你改的代码，有几次可以跑出结果，其他时候还是会导致Julia崩溃。