Implement Shubert-Piyavskii for multimodal deterministic global optimization

[Luis-Varona]

I'd like to create a PR adding the Shubert-Piyavskii method to Optim.jl's univariate solvers. The Shubert-Piyavskii method optimizes a univariate, Lipschitz continuous function inside a specified interval. Given an accurate Lipschitz constant, it is guaranteed to come within an arbitrary epsilon of the true global minimum. Unlike most other deterministic global optimizations algorithm (e.g., golden-section search), it is not restricted to unimodal functions; moreover, it does not require an initial guess. It deterministically samples points from ever-narrowing subintervals of the search space until the sampled best point is sufficiently close to the lower bound given by the Lipschitz constant.

I feel like SP would complement the existing implementations of Brent's method and golden-section search in src/univariate/solvers; both of these are, unfortunately, quite limited in their effectiveness when applied to multimodal functions.

[pkofod]

Do you have a draft version of the code (maybe outside the optim setup)?

[Luis-Varona]

Yes @pkofod I do! Will send later.

[Luis-Varona]

module ShubertPiyavskii

using DataStructures: BinaryMaxHeap

include("error_messages.jl")

export GlobalMaximizationResult, maximize_shubert

const RealT = Float64
# TODO: Another "standard" option is 1e-5... but SP runs too slow then. Decide?
const DEFAULT_TOL = 0.01

struct GlobalMaximizationResult
    maximizer::RealT
    maximum::RealT
    tol::RealT
end

struct _SearchInterval
    x_lower::RealT
    x_upper::RealT
    y_lower::RealT
    y_upper::RealT
    upper_bound::RealT

    function _SearchInterval(
        x_lower::RealT, x_upper::RealT, y_lower::RealT, y_upper::RealT, lipschitz::RealT
    )
        upper_bound = (lipschitz * (x_upper - x_lower) + y_lower + y_upper) / 2
        return new(x_lower, x_upper, y_lower, y_upper, upper_bound)
    end
end

Base.isless(A::_SearchInterval, B::_SearchInterval) = A.upper_bound < B.upper_bound

# TODO: Add parallelization over subintervals of the initial optimization range
function maximize_shubert(
    objective::Function,
    x_lower::RealT,
    x_upper::RealT,
    lipschitz::RealT;
    tol::RealT=DEFAULT_TOL,
)
    x_lower < x_upper || throw(DomainError([x_lower, x_upper], OPTIM_RANGE_ERR))
    lipschitz >= 0 || throw(DomainError(lipschitz, LIPSCHITZ_ERR))
    tol > 0 || throw(DomainError(tol, TOL_ERR))

    y_lower = objective(x_lower)
    y_upper = objective(x_upper)
    maximizer, maximum = (y_lower >= y_upper) ? (x_lower, y_lower) : (x_upper, y_upper)

    search_intervals = BinaryMaxHeap{_SearchInterval}()
    initial_interval = _SearchInterval(x_lower, x_upper, y_lower, y_upper, lipschitz)
    push!(search_intervals, initial_interval)

    while !isempty(search_intervals)
        search_interval = pop!(search_intervals)

        if search_interval.upper_bound - maximum >= tol
            child1, child2 = _split_search_interval(search_interval, objective, lipschitz)

            for child in (child1, child2)
                if child.upper_bound - maximum >= tol
                    push!(search_intervals, child)
                end
            end

            x_sample = child1.x_upper
            y_sample = child1.y_upper

            if maximum <= y_sample
                maximizer, maximum = (x_sample, y_sample)
            end
        end
    end

    return GlobalMaximizationResult(maximizer, maximum, tol)
end

@inline # TODO: Should benchmark to see if this actually helps
function _split_search_interval(
    search_interval::_SearchInterval, objective::Function, lipschitz::RealT
)
    x_lower = search_interval.x_lower
    x_upper = search_interval.x_upper
    y_lower = search_interval.y_lower
    y_upper = search_interval.y_upper

    x_sample = clamp(
        (x_lower + x_upper + (y_upper - y_lower) / lipschitz) / 2, x_lower, x_upper
    )
    y_sample = objective(x_sample)

    child1 = _SearchInterval(x_lower, x_sample, y_lower, y_sample, lipschitz)
    child2 = _SearchInterval(x_sample, x_upper, y_sample, y_upper, lipschitz)
    return child1, child2
end

end

[Luis-Varona]

Certainly, this implementation isn't as well-written as I'd prefer—it' just some draft code from another project of mine. Nevertheless, @pkofod, I hope it gets the point across.

(Indeed, my need for this in a project ultimately not closely related to these types of optimization algorithms is what motivated me to propose that we integrate Shubert–Piyavskii into Optim.jl.)

[pkofod]

I think the error messages file is missing, but thanks :)

[Luis-Varona]

I think the error messages file is missing, but thanks :)

Yes, @pkofod, it's just draft code from another project that already has an error_messages.jl file. Not related to the logic 🙂

Anyway, what are your thoughts on this? Would you like to integrate this algorithm into Optim.jl?

[Luis-Varona]

Hi, @pkofod, just following up on whether you'd actually like to include this in Optim.jl.