QuasiNewtonMethods.jl

Gradient-based numerical optimization.
Author chriselrod
Popularity
1 Star
Updated Last
1 Year Ago
Started In
September 2019

QuasiNewtonMethods

Stable Dev CI CI (Julia nightly) Codecov


This library aims to be fast. It's intended use is for optimizing in statistical logdensity functions, in particular in conjunction with ProbabilityModels.jl and InplaceDHMC.jl (both libraries are still undergoing major development and are not yet usable). The API thus requires using logdensity functions, which ProbabilityModels will automatically define for a given model:

using QuasiNewtonMethods, StrideArrays
using Test

struct Rosenbrock end
function QuasiNewtonMethods.logdensity(::Rosenbrock, θ)
    s = zero(eltype(θ))
    N = length(θ) >> 1
    @inbounds @simd for i  1:N
        s -= 100(θ[i+N] - θ[i]^2)^2 + (1 - θ[i])^2
    end
    if isodd(length(θ))
        @inbounds δ = 1 - θ[end]
        -muladd(δ, δ, -s)
    else
        s
    end
end

function QuasiNewtonMethods.∂logdensity!(∇, ::Rosenbrock, θ) 
    s = zero(eltype(θ))
    N = length(θ) >> 1
    @inbounds @simd for i  1:N
        s -= 100(θ[i+N] - θ[i]^2)^2 + (1 - θ[i])^2
        ∇[i] = 400(θ[i+N] - θ[i]^2)*θ[i] + 2(1 - θ[i])
        ∇[i+N] = 200(θ[i]^2 - θ[i+N])
    end
    if isodd(length(θ))
        @inbounds δ = 1 - θ[end]
        s = -muladd(δ, δ, -s)
        @inbounds ∇[end] = 2δ
    end
    s
end

Example usage, and benchmark vs the equivalent method from Optim.jl:

julia> n = 60 # set size
60

julia> state = QuasiNewtonMethods.BFGSState{n}(undef);

julia> x = @StrideArray randn(StaticInt(n));

julia> @test abs(optimize!(state, Rosenbrock(), x)) < eps()
Test Passed

julia> @show QuasiNewtonMethods.optimum(state) .- 1;
QuasiNewtonMethods.optimum(state) .- 1 = [2.751932015598868e-11, 1.3031797863050087e-12, -1.5009105069907491e-12, 2.6655788687435233e-11, -2.244759933489604e-12, -1.579680830587904e-11, -1.1838985347623066e-10, 1.3630208073323047e-11, 1.982880526441022e-11, 5.3439475067307285e-11, -3.896738487441098e-11, 2.4940494114389367e-11, 2.1896706670077037e-11, -2.1127433136314266e-11, 1.4427570249608834e-11, 2.329803017175891e-11, -3.941846848931618e-12, 3.2440716779547074e-13, -5.52979884105298e-12, 1.6714185591126807e-11, -3.831268635678953e-12, 3.045141916402372e-11, 1.3429257705865894e-12, 1.957722872703016e-11, 9.442890913646806e-12, -4.360312111373332e-11, 2.250799546743565e-11, 1.6193268947972683e-11, -1.954936212911207e-11, -7.409961533255682e-12, 5.451172846449026e-11, 2.3572255258841324e-12, -5.270783809407931e-12, 5.5249804731261065e-11, -4.586775403936372e-12, -3.0561220221159147e-11, -2.37073916053987e-10, 2.8603786006442533e-11, 3.746403187676606e-11, 1.092077539510683e-10, -7.855049943827908e-11, 4.777400697264511e-11, 4.427724853428572e-11, -4.329003822078903e-11, 3.108424628806006e-11, 4.3983705566574827e-11, -6.195599588920686e-12, 3.228528555609955e-13, -8.1592510525752e-12, 3.213473931396038e-11, -7.304934435126142e-12, 5.92046411895808e-11, 3.838485085339016e-12, 4.019495847273902e-11, 1.7372103755519674e-11, -8.670275608579914e-11, 4.6741499559743716e-11, 3.335220988276433e-11, -3.7821412668392895e-11, -1.6192935881065296e-11]

julia> @test all(x -> x  1, QuasiNewtonMethods.optimum(state))
Test Passed

julia> @test maximum(abs, QuasiNewtonMethods.gradient(state)) < 1e-8
Test Passed

julia> @benchmark optimize!($state, Rosenbrock(), $x)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     450.655 μs (0.00% GC)
  median time:      452.243 μs (0.00% GC)
  mean time:        452.309 μs (0.00% GC)
  maximum time:     506.655 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> using Optim, LineSearches

julia> xa = Vector(x);

julia> @benchmark Optim.maximize(x -> QuasiNewtonMethods.logdensity(Rosenbrock(), x), (∇, x) -> QuasiNewtonMethods.∂logdensity!(∇, Rosenbrock(), x), $xa, $(BFGS(linesearch=BackTracking(order=2))))
BenchmarkTools.Trial:
  memory estimate:  291.27 KiB
  allocs estimate:  2483
  --------------
  minimum time:     2.640 ms (0.00% GC)
  median time:      2.657 ms (0.00% GC)
  mean time:        2.673 ms (0.31% GC)
  maximum time:     4.242 ms (34.59% GC)
  --------------
  samples:          1871
  evals/sample:     1

Note that in most problems, evaluating the logdensity function will be the bottleneck, not the speed of the optimization library itself. Thus don't expect a performance improvement like this for real problems. Additionally, QuasiNewtonMethods.jl only provides a backtracking linesearch at the moment. If a different optimization algorithm provides better results, yielding convergence in fewer function evaluations, then again Optim.jl is likely to be faster.

Used By Packages

No packages found.