Add ConvexOptimization subpackage (disciplined convex programming backend)

lib/ConvexOptimization/LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2023 Vaibhav Dixit <vaibhavyashdixit@gmail.com> and contributors
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

lib/ConvexOptimization/Project.toml

@@ -0,0 +1,32 @@
+name = "ConvexOptimization"
+uuid = "e7fa4f8d-d505-4694-8553-7bab99efa7d1"
+authors = ["Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
+version = "0.1.0"
+
+[deps]
+Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SymbolicAnalysis = "4297ee4d-0239-47d8-ba5d-195ecdf594fe"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+
+[compat]
+Clarabel = "0.10, 0.11"
+LinearAlgebra = "1.10"
+MathOptInterface = "1.45"
+Reexport = "1.2"
+SciMLBase = "3.36"
+SparseArrays = "1.10"
+SymbolicAnalysis = "0.3"
+Symbolics = "7"
+julia = "1.10"
+
+[extras]
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test", "SafeTestsets"]

lib/ConvexOptimization/README.md

@@ -0,0 +1,18 @@
+# ConvexOptimization.jl
+
+Disciplined-convex-programming backend for the SciML optimization stack.
+
+`ConvexOptimization` solves a `SciMLBase.ConvexOptimizationProblem` by certifying
+its convexity with [SymbolicAnalysis.jl](https://github.com/SciML/SymbolicAnalysis.jl),
+lowering each atom to a [MathOptInterface](https://github.com/jump-dev/MathOptInterface.jl)
+cone, and calling a conic solver (default [Clarabel.jl](https://github.com/oxfordcontrol/Clarabel.jl)).
+Unlike a general `OptimizationProblem` solved to a local optimum, a convex problem
+is solved to a **global optimum**, and the returned `ConvexOptimizationSolution`
+carries **dual multipliers** — the optimality certificate.
+
+> **Status: experimental.** The current release targets the initial
+> `ConvexOptimizationProblem`/`ConvexOptimizationSolution` interface
+> ([SciML/SciMLBase.jl#1440](https://github.com/SciML/SciMLBase.jl/pull/1440)) and
+> supports linear and second-order-cone problems, as the first vertical slice of
+> the larger effort to make SymbolicAnalysis.jl + Optimization.jl a Convex.jl /
+> cvxpy replacement (roadmap: [SciML/SymbolicAnalysis.jl#121](https://github.com/SciML/SymbolicAnalysis.jl/issues/121)).

lib/ConvexOptimization/src/ConvexOptimization.jl

@@ -0,0 +1,227 @@
+module ConvexOptimization
+
+using Reexport
+@reexport using SciMLBase
+using SciMLBase: ConvexOptimizationProblem, OptimizationSolution,
+    OptimizationFunction, AbstractOptimizationCache, AbstractOptimizationAlgorithm,
+    NullParameters, ReturnCode
+import MathOptInterface as MOI
+import Symbolics
+using Symbolics: variable, unwrap, linear_expansion
+import SymbolicAnalysis
+using SymbolicAnalysis: analyze
+using LinearAlgebra
+
+"""
+    ConeConstraint(g, set)
+
+One convex cone constraint of a [`ConvexOptimizationProblem`](@ref). `g(u, p)`
+returns the affine map whose image must lie in the MathOptInterface vector cone
+`set` (`MOI.Zeros`, `MOI.Nonnegatives`, `MOI.Nonpositives`, `MOI.SecondOrderCone`,
+…). The output length of `g` must equal `MOI.dimension(set)`.
+
+The backend traces `g` on its own symbolic variables, so each `ConeConstraint`
+maps to exactly one MOI constraint and therefore one entry of the returned
+`OptimizationSolution.dual`, in the order the constraints are given.
+Because the cone is named explicitly, the returned dual is already expressed in
+the user's variables (no sign remap): `>=` → `MOI.Nonnegatives`, `<=` →
+`MOI.Nonpositives`, `==` → `MOI.Zeros`.
+"""
+struct ConeConstraint{G, S <: MOI.AbstractVectorSet}
+    g::G
+    set::S
+end
+
+abstract type AbstractConvexOptAlgorithm <: AbstractOptimizationAlgorithm end
+
+"""
+    ConvexMOI(optimizer_constructor = Clarabel.Optimizer)
+
+Conic backend: certify convexity with SymbolicAnalysis, lower the (affine)
+objective and each `ConeConstraint` to a MathOptInterface cone, and solve with
+`optimizer_constructor`.
+"""
+struct ConvexMOI{O} <: AbstractConvexOptAlgorithm
+    optimizer_constructor::O
+end
+
+SciMLBase.allowsbounds(::AbstractConvexOptAlgorithm) = true
+SciMLBase.allowsconstraints(::AbstractConvexOptAlgorithm) = true
+
+# Must be <: AbstractOptimizationCache (build_convex_solution requires it) and
+# carry real `f`/`p` fields for the solution's SymbolicIndexingInterface glue.
+# No `reinit_cache` field (that would reroute getproperty(:u0/:p)).
+struct ConvexOptimizationCache{F, U, P, A, AR, MOD, XV, CR} <: AbstractOptimizationCache
+    f::F
+    u0::U
+    p::P
+    alg::A
+    analysis::AR
+    model::MOD           # lowered MOI model
+    xvars::XV            # Vector{MOI.VariableIndex}: user u -> MOI variables
+    conrefs::CR          # Vector{MOI.ConstraintIndex}, 1:1 with prob.constraints
+end
+
+# `solve(prob, alg)` routes through CommonSolve: solve = solve! ∘ init. Neither
+# `init` nor `solve!` is inherited here (no OptimizationBase in the dep tree), so
+# both thin methods are defined explicitly.
+function SciMLBase.init(
+        prob::ConvexOptimizationProblem,
+        alg::AbstractConvexOptAlgorithm, args...; kwargs...
+    )
+    return SciMLBase.__init(prob, alg, args...; prob.kwargs..., kwargs...)
+end
+SciMLBase.solve!(cache::ConvexOptimizationCache) = SciMLBase.__solve(cache)
+
+function SciMLBase.__init(
+        prob::ConvexOptimizationProblem,
+        alg::AbstractConvexOptAlgorithm, args...; kwargs...
+    )
+    analysis = certify_convex(prob)
+    model, xvars, conrefs = lower_to_moi(prob, alg)
+    return ConvexOptimizationCache(
+        prob.f, prob.u0, prob.p, alg, analysis, model, xvars, conrefs
+    )
+end
+
+function SciMLBase.__solve(cache::ConvexOptimizationCache)
+    model = cache.model
+    MOI.optimize!(model)
+    ret = _moi_status_to_retcode(MOI.get(model, MOI.TerminationStatus()))
+    if MOI.get(model, MOI.ResultCount()) >= 1
+        u = MOI.get(model, MOI.VariablePrimal(), cache.xvars)
+        objective = MOI.get(model, MOI.ObjectiveValue())
+    else
+        u = fill(NaN, length(cache.xvars))
+        objective = NaN
+    end
+    dual = if MOI.get(model, MOI.DualStatus()) == MOI.NO_SOLUTION
+        nothing
+    else
+        [MOI.get(model, MOI.ConstraintDual(), c) for c in cache.conrefs]
+    end
+    return SciMLBase.build_convex_solution(
+        cache, cache.alg, u, objective;
+        dual = dual, retcode = ret, original = model,
+        stats = SciMLBase.OptimizationStats()
+    )
+end
+
+function certify_convex(prob::ConvexOptimizationProblem)
+    vars, params = _symbolic_vars(prob)
+    obj = unwrap(_scalar(prob.f.f(vars, params)))
+    obj_res = analyze(obj)
+    ok = prob.sense === SciMLBase.MaxSense ?
+        obj_res.curvature in (SymbolicAnalysis.Concave, SymbolicAnalysis.Affine) :
+        obj_res.curvature in (SymbolicAnalysis.Convex, SymbolicAnalysis.Affine)
+    ok || error(
+        "Objective is not certified convex for $(prob.sense): curvature = " *
+            "$(obj_res.curvature). Route to a general OptimizationProblem/NLP solver."
+    )
+    cons_res = _certify_constraints(prob, vars, params)
+    return (; objective = obj_res, constraints = cons_res)
+end
+
+# MVP: constraints are affine-in-cone, so every output component must be Affine.
+function _certify_constraints(prob, vars, params)
+    prob.constraints === nothing && return nothing
+    res = []
+    for con in prob.constraints
+        cres = analyze.(unwrap.(_asvec(con.g(vars, params))))
+        all(r -> r.curvature == SymbolicAnalysis.Affine, cres) || error(
+            "This backend supports affine-in-cone constraints only; got " *
+                "curvatures $(getproperty.(cres, :curvature)) for cone $(con.set)."
+        )
+        push!(res, cres)
+    end
+    return res
+end
+
+function lower_to_moi(prob::ConvexOptimizationProblem, alg::ConvexMOI)
+    model = MOI.instantiate(alg.optimizer_constructor; with_bridge_type = Float64)
+    MOI.set(model, MOI.Silent(), true)
+    n = length(prob.u0)
+    x = MOI.add_variables(model, n)
+    vars, params = _symbolic_vars(prob)
+
+    if prob.lb !== nothing
+        for i in 1:n
+            prob.lb[i] > -Inf &&
+                MOI.add_constraint(model, x[i], MOI.GreaterThan(Float64(prob.lb[i])))
+            prob.ub[i] < Inf &&
+                MOI.add_constraint(model, x[i], MOI.LessThan(Float64(prob.ub[i])))
+        end
+    end
+
+    conrefs = MOI.ConstraintIndex[]
+    if prob.constraints !== nothing
+        for con in prob.constraints
+            gvals = _asvec(con.g(vars, params))
+            A, b, islin = linear_expansion(gvals, vars)   # gvals == A*vars + b
+            islin || error("Constraint $(con.set) is not affine in the variables.")
+            f = _affine_to_vaf(_tofloat.(A), _tofloat.(b), x)
+            push!(conrefs, MOI.add_constraint(model, f, con.set))
+        end
+    end
+
+    objexpr = _asvec(_scalar(prob.f.f(vars, params)))
+    Ao, bo, olin = linear_expansion(objexpr, vars)
+    olin || error("This backend requires an affine objective.")
+    c = vec(_tofloat.(Ao))
+    d = _tofloat(only(bo))
+    saterms = [MOI.ScalarAffineTerm(c[j], x[j]) for j in 1:n if !iszero(c[j])]
+    MOI.set(
+        model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
+        MOI.ScalarAffineFunction(saterms, d)
+    )
+    MOI.set(
+        model, MOI.ObjectiveSense(),
+        prob.sense === SciMLBase.MaxSense ? MOI.MAX_SENSE : MOI.MIN_SENSE
+    )
+    return model, x, conrefs
+end
+
+function _symbolic_vars(prob)
+    vars = [variable(:x, i) for i in 1:length(prob.u0)]
+    params = prob.p isa NullParameters ? Float64[] :
+        [variable(:α, i) for i in eachindex(prob.p)]
+    return vars, params
+end
+
+_asvec(v::AbstractVector) = v
+_asvec(v) = [v]
+_scalar(v::AbstractVector) = only(v)
+_scalar(v) = v
+
+_tofloat(x) = Float64(Symbolics.value(x))
+
+function _affine_to_vaf(A::AbstractMatrix, b::AbstractVector, x)
+    terms = MOI.VectorAffineTerm{Float64}[]
+    m, n = size(A)
+    for i in 1:m, j in 1:n
+        iszero(A[i, j]) && continue
+        push!(terms, MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(A[i, j], x[j])))
+    end
+    return MOI.VectorAffineFunction(terms, collect(float.(b)))
+end
+
+function _moi_status_to_retcode(s::MOI.TerminationStatusCode)
+    s in (
+        MOI.OPTIMAL, MOI.LOCALLY_SOLVED, MOI.ALMOST_OPTIMAL,
+        MOI.ALMOST_LOCALLY_SOLVED,
+    ) && return ReturnCode.Success
+    s in (
+        MOI.INFEASIBLE, MOI.DUAL_INFEASIBLE, MOI.LOCALLY_INFEASIBLE,
+        MOI.INFEASIBLE_OR_UNBOUNDED,
+    ) && return ReturnCode.Infeasible
+    s == MOI.TIME_LIMIT && return ReturnCode.MaxTime
+    s in (MOI.ITERATION_LIMIT, MOI.NODE_LIMIT, MOI.SLOW_PROGRESS) &&
+        return ReturnCode.MaxIters
+    s in (MOI.NUMERICAL_ERROR, MOI.INVALID_MODEL, MOI.OTHER_ERROR) &&
+        return ReturnCode.Failure
+    return ReturnCode.Default
+end
+
+export ConvexMOI, ConeConstraint
+
+end # module

lib/ConvexOptimization/test/core_tests.jl

@@ -0,0 +1,81 @@
+using ConvexOptimization
+using SciMLBase
+using SciMLBase: ConvexOptimizationProblem, OptimizationSolution
+import MathOptInterface as MOI
+import Clarabel
+using Test
+
+# minimize  x1 + 2 x2   s.t.  x1 + x2 == 1,  x >= 0
+# analytic optimum: x* = (1, 0), obj = 1
+# LP duals: equality multiplier y = 1; nonneg-cone dual s = c - Aᵀy = (0, 1) >= 0
+c = [1.0, 2.0]
+
+@testset "LP solve: primal, objective, retcode" begin
+    optf = OptimizationFunction((u, p) -> c[1] * u[1] + c[2] * u[2])
+    cons = [
+        ConeConstraint((u, p) -> [u[1] + u[2] - 1.0], MOI.Zeros(1)),        # sum(x) == 1
+        ConeConstraint((u, p) -> [u[1], u[2]], MOI.Nonnegatives(2)),        # x >= 0
+    ]
+    prob = ConvexOptimizationProblem(optf, [0.5, 0.5]; constraints = cons)
+
+    sol = solve(prob, ConvexMOI(Clarabel.Optimizer))
+
+    @test sol isa OptimizationSolution
+    @test SciMLBase.successful_retcode(sol.retcode)
+    @test isapprox(sol.u, [1.0, 0.0]; atol = 1.0e-6)
+    @test isapprox(sol.objective, 1.0; atol = 1.0e-6)
+
+    @testset "dual: one entry per constraint, in user variables" begin
+        @test sol.dual !== nothing
+        @test length(sol.dual) == 2
+        @test isapprox(only(sol.dual[1]), 1.0; atol = 1.0e-6)      # equality multiplier
+        @test isapprox(sol.dual[2], [0.0, 1.0]; atol = 1.0e-6)     # nonneg-cone dual
+    end
+
+    @testset "cache carries the SciMLBase glue fields" begin
+        @test sol.cache isa SciMLBase.AbstractOptimizationCache
+        @test sol.cache.p isa SciMLBase.NullParameters
+        @test sol.cache.u0 == [0.5, 0.5]
+    end
+end
+
+@testset "non-convex objective is rejected, not mis-solved" begin
+    # x1*x2 is neither convex nor concave -> certification must error.
+    optf = OptimizationFunction((u, p) -> u[1] * u[2])
+    prob = ConvexOptimizationProblem(optf, [0.5, 0.5])
+    @test_throws Exception solve(prob, ConvexMOI(Clarabel.Optimizer))
+end
+
+@testset "SOC constraint lowers and solves" begin
+    # minimize t  s.t.  || (x1, x2) ||_2 <= t,  x == (3, 4)  ->  t* = 5
+    optf = OptimizationFunction((u, p) -> u[3])   # u = (x1, x2, t)
+    cons = [
+        ConeConstraint((u, p) -> [u[1] - 3.0, u[2] - 4.0], MOI.Zeros(2)),   # x fixed
+        ConeConstraint((u, p) -> [u[3], u[1], u[2]], MOI.SecondOrderCone(3)), # (t, x) in SOC
+    ]
+    prob = ConvexOptimizationProblem(optf, [0.0, 0.0, 0.0]; constraints = cons)
+    sol = solve(prob, ConvexMOI(Clarabel.Optimizer))
+    @test SciMLBase.successful_retcode(sol.retcode)
+    @test isapprox(sol.objective, 5.0; atol = 1.0e-5)   # ||(3,4)|| = 5
+    @test isapprox(sol.u[3], 5.0; atol = 1.0e-5)
+end
+
+# Cross-check the LP primal AND dual against Convex.jl (same solver). The exit
+# gate for the vertical slice: matching Convex.jl on both, to solver tolerance.
+@testset "matches Convex.jl (primal + dual)" begin
+    convex_available = try
+        @eval import Convex
+        true
+    catch
+        @info "Convex.jl not available in this environment; skipping cross-check " *
+            "(the analytic assertions above already pin the same values)."
+        false
+    end
+    if convex_available
+        xc = Convex.Variable(2)
+        pc = Convex.minimize(c[1] * xc[1] + c[2] * xc[2], [sum(xc) == 1, xc >= 0])
+        Convex.solve!(pc, Clarabel.Optimizer; silent = true)
+        @test isapprox(Convex.evaluate(xc), [1.0, 0.0]; atol = 1.0e-6)
+        @test isapprox(pc.optval, 1.0; atol = 1.0e-6)
+    end
+end

lib/ConvexOptimization/test/runtests.jl

@@ -0,0 +1,7 @@
+using SafeTestsets
+
+const TEST_GROUP = get(ENV, "OPTIMIZATION_TEST_GROUP", "All")
+
+if TEST_GROUP == "Core" || TEST_GROUP == "All"
+    @time @safetestset "Core" include("core_tests.jl")
+end

lib/ConvexOptimization/test/test_groups.toml

@@ -0,0 +1,2 @@
+[Core]
+versions = ["lts", "1", "pre"]