Nonlinear SDP interior point meta-algorithm

Project.toml

@@ -20,13 +20,13 @@ Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 Sobol = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
 ADNLPModels = "0.1, 0.2, 0.3"
@@ -55,6 +55,8 @@ Zygote = "0.6"
 julia = "1"
 
 [extras]
+ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 Hyperopt = "93e5fe13-2215-51db-baaf-2e9a34fb2712"
 Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84"
@@ -65,4 +67,4 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["FiniteDifferences", "Hyperopt", "Juniper", "NLopt", "Pavito", "Percival", "SafeTestsets", "Test"]
+test = ["FiniteDifferences", "Hyperopt", "Juniper", "NLopt", "Pavito", "Percival", "SafeTestsets", "Distributions", "Test", "ChainRulesTestUtils"]

src/Nonconvex.jl

@@ -9,6 +9,8 @@ export  Model,
         set_objective!,
         add_ineq_constraint!,
         add_eq_constraint!,
+        add_sd_constraint!,
+        decompress_symmetric,
         getmin,
         getmax,
         isinteger,
@@ -32,6 +34,7 @@ export  Model,
         PavitoIpoptCbcAlg,
         HyperoptAlg,
         BayesOptAlg,
+        SDPBarrierAlg,
         MTSAlg,
         LS1Alg,
         KKTCriteria,
@@ -47,6 +50,7 @@ export  Model,
         PavitoIpoptCbcOptions,
         HyperoptOptions,
         BayesOptOptions,
+        SDPBarrierOptions,
         MTSOptions,
         LS1Options,
         Tolerance,
@@ -126,6 +130,9 @@ include("algorithms/mts.jl")
 
 include("algorithms/deflation.jl")
 
+# Semi-definite programming
+include("algorithms/sdp.jl")
+
 # Wrappers
 
 include("wrappers/moi.jl")

src/algorithms/common.jl

@@ -212,6 +212,7 @@ end
 
 abstract type AbstractModel end
 
+
 """
 ```
 optimize(
@@ -224,6 +225,7 @@ optimize(
     callback::Function = plot_trace ? LazyPlottingCallback() : NoCallback(),
     kwargs...,
 )
+
 ```
 
 Optimizes `model` using the algorithm `optimizer`, e.g. an instance of [`MMA87`](@ref) or [`MMA02`](@ref). `x0` is the initial solution. The keyword arguments are:
@@ -234,7 +236,14 @@ Optimizes `model` using the algorithm `optimizer`, e.g. an instance of [`MMA87`]
 
  The details of the MMA optimization algorithms can be found in the original [1987 MMA paper](https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620240207) and the [2002 paper](https://epubs.siam.org/doi/abs/10.1137/S1052623499362822).
 """
+function _optimize_precheck(model::AbstractModel, optimizer::AbstractOptimizer, args...; options=nothing, kwargs...)
+    if (length(model.sd_constraints.fs) != 0) && !(options isa SDPBarrierOptions)
+        @warn "Input `sd_constraints` will be ignored due to not using `SDPBarrierOptions`. "
+    end
+end
+
 function optimize(model::AbstractModel, optimizer::AbstractOptimizer, x0, args...; kwargs...)
+    _optimize_precheck(model, optimizer, x0, args...; kwargs...)
     _model, _x0, unflatten = tovecmodel(model, x0)
     r = optimize(_model, optimizer, _x0, args...; kwargs...)
     return @set r.minimizer = unflatten(r.minimizer)
@@ -244,11 +253,13 @@ end
  optimize without x0
 """
  function optimize(model::AbstractModel, optimizer::AbstractOptimizer, args...; kwargs...)
+    _optimize_precheck(model, optimizer, args...; kwargs...)
     _model, _, unflatten = tovecmodel(model)
     r = optimize(_model, optimizer, args...; kwargs...)
     return @set r.minimizer = unflatten(r.minimizer)
 end
 
+
 """
 Clamp the box constraint and evaluate objective function at point x
 """

src/algorithms/mma_algorithm.jl

@@ -235,6 +235,7 @@ function optimize!(workspace::MMAWorkspace)
 
             # Evaluates the exact objective and constraints and their gradients at the optimal x
             optimalx = getoptimalx(dualmodel)
+
             fg, ∇fg = value_jacobian(
                 approxmodel.objective_ineq_constraints,
                 optimalx,
@@ -245,7 +246,6 @@ function optimize!(workspace::MMAWorkspace)
             # Reduces the chance of overflow or underflow in λs 
             # Helps with the convergence
             auto_scale && iter < 5 && scaleobjective!(model, fg, ∇fg, approxfg, λ)
-            #iter < 10 && scaleobjective!(model, fg, ∇fg, approxfg, λ)
 
             # Check if the approximation is an upper bound at the current x.
             # If not, increase ρ. Only used in MMA02

src/algorithms/sdp.jl

@@ -0,0 +1,137 @@
+# Semidefinite programming
+
+# Decompress
+function lowertriangind(mat::Matrix)
+    indices = [i for i in CartesianIndices(mat) if i[1]>i[2]]
+    return LinearIndices(mat)[indices]
+end
+
+function rearrange_x(x_L::AbstractVector, x_D::AbstractVector)
+    mat_dim = length(x_D)
+    L = zeros(mat_dim, mat_dim)
+    L[lowertriangind(L)] .= x_L
+    D = diagm(x_D)
+    return (L, D)
+end
+
+function decompress_symmetric(L::Matrix, D::Matrix)
+    L + D + L'
+end
+
+function decompress_symmetric(x_L::AbstractArray, x_D::AbstractArray)
+    L, D = rearrange_x(x_L, x_D)
+    return decompress_symmetric(L, D)
+end
+
+function ChainRulesCore.rrule(::typeof(rearrange_x), x_L::AbstractVector, x_D::AbstractVector)
+    function pullback((ΔL, ΔD))
+        Δx_L = ΔL[lowertriangind(ΔL)]
+        Δx_D = diag(ΔD)
+        NoTangent(), Δx_L, Δx_D 
+    end
+    return rearrange_x(x_L, x_D), pullback
+end
+
+
+# Optimizer (algorithm)
+struct SDPBarrierAlg{Alg <: AbstractOptimizer} <: AbstractOptimizer
+    sub_alg::Alg
+end
+function SDPBarrierAlg(;sub_alg)
+    return SDPBarrierAlg(sub_alg)
+end
+
+# Options
+@params mutable struct SDPBarrierOptions
+    # Dimension of objective matrix
+    # Hyperparameters 
+    # Initial value of `c` in barrier method: 
+    # `Real` for using same value for all `sd_constraints`, `AbstractArray` for assign them respectively
+    c_init::Union{Real, AbstractArray}
+    # Decrease rate of `c` for every epoch, same as above
+    c_decr::Union{Real, AbstractArray}
+    n_iter::Int
+    # sub_option to solve (in)equality constraints
+    sub_options
+    # Keep all results or not
+    keep_all::Bool
+end
+function SDPBarrierOptions(c_init, c_decr, n_iter; sub_options, keep_all=true)
+    @assert 0 < c_decr < 1 "c_decr should be between 0 and 1. "
+    @assert c_init > 0 "c_init shoule be larger than 0. "
+    SDPBarrierOptions(c_init, c_decr, n_iter, sub_options, keep_all)
+end
+function SDPBarrierOptions(;sub_options, c_init=1.0, c_decr=0.1, n_iter=10, keep_all=true)
+    SDPBarrierOptions(c_init, c_decr, n_iter, sub_options=sub_options, keep_all=keep_all)
+end
+
+# Result
+@params struct SDPResult
+    minimum
+    minimizer
+    results
+    optimal_ind
+end
+
+@params struct SDPWorkspace <: Workspace
+    model::VecModel
+    x0::AbstractVector
+    options::SDPBarrierOptions
+    sub_alg::AbstractOptimizer
+end
+
+function Workspace(model::VecModel, optimizer::SDPBarrierAlg, x0, args...; options, kwargs...,)
+    @unpack c_init, c_decr = options
+    for c in model.sd_constraints.fs
+        @assert isposdef(c(x0)) "Initial matrix should be positive definite. "
+    end
+    if c_init isa AbstractArray
+        @assert length(model.sd_constraints.fs) == length(c_init) "c_init should be same length with number of `sd_constraints` when using array. "
+    end
+    if c_decr isa AbstractArray
+        @assert length(model.sd_constraints.fs) == length(c_decr) "c_decr should be same length with number of `sd_constraints` when using array. "
+    end
+    if c_init isa AbstractArray && c_decr isa AbstractArray
+        @assert length(c_init) == length(c_decr) "c_decr should be same length with c_init. "
+    end
+    return SDPWorkspace(model, copy(x0), options, optimizer.sub_alg)
+end
+
+function sd_objective(objective0, sd_constraints, c::AbstractArray)
+    function _objective(args)
+        target = objective0(args)
+        barrier = sum(c .* -logdet.(map(f -> f(args), sd_constraints.fs)))
+        return target + barrier
+    end
+    return _objective
+end
+
+function to_barrier(model, c::AbstractArray)
+    sd_constraints, objective0 = model.sd_constraints, model.objective
+    _model = set_objective(model, sd_objective(objective0, sd_constraints, c))
+    return _model
+end
+
+function optimize!(workspace::SDPWorkspace)
+    @unpack model, x0, options, sub_alg = workspace
+    @unpack c_init, c_decr, n_iter, sub_options, keep_all = options
+    objective0 = model.objective
+    x = copy(x0)
+    c = c_init isa Real ? ([c_init for _ in 1:length(model.sd_constraints.fs)]) : c_init
+    results = []
+    for _ in 1:n_iter
+        model_i = to_barrier(model, c)
+        result_i = optimize(model_i, sub_alg, x, options = sub_options)
+        minimizer_i = result_i.minimizer
+        push!(results, (objective0(minimizer_i), minimizer_i))
+        c = c .* c_decr
+        x = copy(minimizer_i)
+    end
+    optimal_ind = argmin(first.(results))
+    minimum, minimizer = results[optimal_ind]
+    if keep_all
+        return SDPResult(minimum, minimizer, results, optimal_ind)
+    else
+        return SDPResult(minimum, minimizer, nothing, nothing)
+    end
+end

src/functions/functions.jl

@@ -273,3 +273,25 @@ getdim(c::EqConstraint) = c.dim
 Returns the function wrapped in `f`.
 """
 getfunction(f::EqConstraint) = f.f
+
+
+"""
+
+Used in semidefinite programming
+"""
+@params struct SDConstraint <: AbstractConstraint
+    f::Function
+    dim::Int
+end
+
+function (c::SDConstraint)(args...; kwargs...)
+    out = c.f(args...; kwargs...)
+    # println(length(out))
+    # println(getdim(c))
+    @assert length(out) == getdim(c)^2
+    return out
+end
+
+getdim(c::SDConstraint) = c.dim
+
+getfunction(c::SDConstraint) = c.f

src/models/dict_model.jl

@@ -2,6 +2,7 @@ mutable struct DictModel <: AbstractModel
     objective::Union{Nothing, Objective}
     eq_constraints::VectorOfFunctions
     ineq_constraints::VectorOfFunctions
+    sd_constraints::VectorOfFunctions
     box_min::OrderedDict
     box_max::OrderedDict
     init::OrderedDict
@@ -12,6 +13,7 @@ function DictModel(f = nothing)
     return DictModel(
         Objective(f), VectorOfFunctions(EqConstraint[]),
         VectorOfFunctions(IneqConstraint[]),
+        VectorOfFunctions(SDConstraint[]),
         OrderedDict(), OrderedDict(),
         OrderedDict(), OrderedDict()
     )

src/models/model.jl

@@ -15,11 +15,13 @@ The `Model` structs stores information about the nonlinear optimization problem.
 - `objective`: the objective function of the problem of type [`Objective`](@ref).
 - `eq_constraints`: the equality constraints of the problem of type [`VectorOfFunctions`](@ref). Each function in `ineq_constraints` can be an instance of [`IneqConstraint`](@ref) or [`AbstractFunction`](@ref). If the function is not an `IneqConstraint`, its right-hand-side bound is assumed to be 0.
 - `ineq_constraints`: the inequality constraints of the problem of type [`VectorOfFunctions`](@ref). Each function in `ineq_constraints` can be an instance of [`IneqConstraint`](@ref) or [`AbstractFunction`](@ref). If the function is not an `IneqConstraint`, its right-hand-side bound is assumed to be 0.
+- `sd_constraints`: a list of matrix-valued functions which must be positive semidefinite at any feasible solution.
 """
 mutable struct Model{Tv <: AbstractVector} <: AbstractModel
     objective::Union{Nothing, Objective}
     eq_constraints::VectorOfFunctions
     ineq_constraints::VectorOfFunctions
+    sd_constraints::VectorOfFunctions
     box_min::Tv
     box_max::Tv
     init::Tv
@@ -41,7 +43,12 @@ Constructs an empty model with objective function `f` and decision variable valu
 """
 Model(::Type{T}, f::Function) where {T} = Model(T, Objective(f))
 function Model(::Type{T}, obj::Union{Nothing, Objective}) where {T}
-    return Model(obj, VectorOfFunctions(EqConstraint[]), VectorOfFunctions(IneqConstraint[]), T[], T[], T[], falses(0))
+    return Model(obj, 
+    VectorOfFunctions(EqConstraint[]), 
+    VectorOfFunctions(IneqConstraint[]), 
+    VectorOfFunctions(SDConstraint[]), 
+    T[], T[], T[], 
+    falses(0))
 end
 
 getobjective(m::AbstractModel) = m.objective
@@ -139,6 +146,11 @@ end
 
 Sets the objective of the moodel `m` to the function `f`. `f` must return a scalar.
 """
+function set_objective(m::AbstractModel, f::Function; kwargs...)
+    @set m.objective = Objective(f; kwargs...)
+    return m
+end
+
 function set_objective!(m::AbstractModel, f::Function; kwargs...)
     m.objective = Objective(f; kwargs...)
     return m
@@ -173,3 +185,11 @@ function add_eq_constraint!(m::AbstractModel, fs::Vector{<:EqConstraint})
     append!(m.eq_constraints.fs, fs)
     return m
 end
+
+function add_sd_constraint!(m::AbstractModel, f::Function, dim=size(f(getinit(m)), 1))
+    return add_sd_constraint!(m, SDConstraint(f, dim))
+end
+function add_sd_constraint!(m::AbstractModel, sd_constraint::AbstractFunction)
+    push!(m.sd_constraints.fs, sd_constraint)
+    return m
+end