Chordal Decomposition of SDP

[samuelsonric]

Dear JuMP developers,

I would like to try implementing a chordal decomposition layer for MathOptInterface. If there is still interest in this, I would appreciate some help understanding how to go about it.

I imagine that such a layer would work something like this. Given a semidefinite optimization problem...

1. Convert the problem to a normal form using bridges.

minimize   C  ⋅ X
subject to Aᵢ ⋅ X = bᵢ, i = 1, …, m
           X ∈ 𝕊

2. Perform a chordal decomposition of the transformed problem.
3. Construct a new problem with the decomposed constraints.

I have already implemented 2. Is there an easy way to perform 1 and 3?

[samuelsonric]

This is a duplicate of the issue I made here.

[odow]

@blegat has probably thought about this more than me.

What information do you need as input and output? Can you support Ax + b in S?

I would probably do something like https://github.com/jump-dev/MultiObjectiveAlgorithms.jl, except that instead of intercepting the objective, we intercept the add_constraint(model, f, ::PositiveSemidefiniteConeTriangle) constraints.

The API would be something like:

using JuMP, SCS, ChordalDecomposition
model = Model(() -> ChordalDecomposition.Optimizer(SCS.Optimizer))

[odow]

I'm not sure how much progress @tjdiamandis made in https://github.com/tjdiamandis/Chordal.jl

[odow]

There's also a question whether to expose this as a layer, or whether to expose this as a set of utility functions for people to decompose a matrix at the modeling level.

[tjdiamandis]

As far as I recall, I just implemented utilities in that package: the chordal decomposition and subsequent matrix completion after the solve.

[samuelsonric]

@blegat has probably thought about this more than me.

What information do you need as input and output? Can you support Ax + b in S?

I would probably do something like https://github.com/jump-dev/MultiObjectiveAlgorithms.jl, except that instead of intercepting the objective, we intercept the add_constraint(model, f, ::PositiveSemidefiniteConeTriangle) constraints.

The API would be something like:

using JuMP, SCS, ChordalDecomposition
model = Model(() -> ChordalDecomposition.Optimizer(SCS.Optimizer))

To perform a decomposition, you need to have both the objective and the constraints, since you are decomposing their "aggregate sparsity structure".

https://arxiv.org/abs/1911.05615

[blegat]

To get the full model information at once, you can follow what Dualization.jl is doing. So implement copy_to like
https://github.com/jump-dev/Dualization.jl/blob/42cb27f0122dcc06439e37d38e6034f5143ba0d8/src/MOI_wrapper.jl#L209
We call this the non-incremental interface, see
https://jump.dev/JuMP.jl/stable/moi/tutorials/implementing/#The-non-incremental-interface

[samuelsonric]

I'm not sure how much progress @tjdiamandis made in https://github.com/tjdiamandis/Chordal.jl

This is the package that I intend to use; it does more-or-less the same thing as Chordal.jl.

[samuelsonric]

I'll reopen this when I make some progress.

[samuelsonric]

@odow I was wrong: this can be handled at the constraint level. One is rewriting a constraint

$$S = \left\{ x \in \mathbb{R}^m \mid Ax + b \in \mathbb{S}_+^n \right\}$$

as the equivalent constraint

$$S = \left\{ x \in \mathbb{R}^m \mid Ax + b \in T_1 \mathbb{S}_+^{n_i} + \dots + T_k \mathbb{S}_+^{n_k} \right\}$$

where $$n_i \leq n$$ for all $$1 \leq i \leq k$$. Should this be implemented as a bridge?

[odow]

The main issue is that bridges can't have a loop where they take in a constraint and produce a new one. So we have two options:

1. add a new set, like SparsePositiveSemidefiniteConeTriangle, and bridge that to PositiveSemidefiniteConeTriangle. The user then needs to opt-in by writing JuMP.@constraint(model, x in SparsePositiveSemidefiniteConeTriangle()) instead of @constraint(model, x in PSDCone()).
2. go the ChordalDecomposition.Optimizer route, where we intercept every PSD constraint, potentially decompose, and then pass to an inner optimizer.

If there are multiple options or algorithms that could be tuned or swapped out in the chordal decomposition, we should go with 2. We should probably even go with 2 regardless, and once the kinks are ironed out we could consider a set/bridge.

But the bar for introducing yet-another-psd-set is very high and gets higher all the time...

[joaquimg]

We should certainly go with 2

[odow]

A basic implementation to start with would be:

module MathOptChordalDecomposition

import MathOptInterface as MOI

const CI = MOI.ConstraintIndex{
    MOI.VectorAffineFunction{Float64},
    MOI.PositiveSemidefiniteConeTriangle,
}

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::MOI.AbstractOptimizer
    outer_to_inner::Dict{CI,Vector{CI}}

    function Optimizer(optimizer_factory)
        return new(
            MOI.instantiate(
                optimizer_factory;
                with_cache_type = Float64,
                with_bridge_type = Float64,
            ),
            Dict{CI,Vector{CI}}(),
        )
    end
end

function MOI.empty!(model::Optimizer)
    MOI.empty!(model.inner)
    return
end

function MOI.is_empty(model::Optimizer)
    return MOI.is_empty(model.inner)
end

MOI.supports_incremental_interface(::Optimizer) = true

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
    return MOI.Utilities.default_copy_to(dest, src)
end

MOI.optimize!(model::Optimizer) = MOI.optimize!(model.inner)

const _ATTRIBUTES = Union{
    MOI.AbstractConstraintAttribute,
    MOI.AbstractModelAttribute,
    MOI.AbstractOptimizerAttribute,
    MOI.AbstractVariableAttribute,
}

function MOI.set(model::Optimizer, attr::_ATTRIBUTES, args...)
    MOI.set(model.inner, attr, args...)
    return
end

function MOI.get(model::Optimizer, attr::_ATTRIBUTES, args...)
    if MOI.is_set_by_optimize(attr)
        msg = "MOCD does not support querying this attribute."
        throw(MOI.GetAttributeNotAllowed(attr, msg))
    end
    return MOI.get(model.inner, attr, args...)
end

function MOI.get(model::Optimizer, attr::_ATTRIBUTES, arg::Vector{T}) where {T}
    if MOI.is_set_by_optimize(attr)
        msg = "MOCD does not support querying this attribute."
        throw(MOI.GetAttributeNotAllowed(attr, msg))
    end
    return MOI.get.(model, attr, arg)
end

### AbstractOptimizerAttribute

function MOI.supports(model::Optimizer, arg::MOI.AbstractOptimizerAttribute)
    return MOI.supports(model.inner, arg)
end

function MOI.set(model::Optimizer, attr::MOI.AbstractOptimizerAttribute, value)
    MOI.set(model.inner, attr, value)
    return
end

function MOI.get(model::Optimizer, attr::MOI.AbstractOptimizerAttribute)
    return MOI.get(model.inner, attr)
end

### AbstractModelAttribute

function MOI.supports(model::Optimizer, arg::MOI.AbstractModelAttribute)
    return MOI.supports(model.inner, arg)
end

### AbstractVariableAttribute

function MOI.is_valid(model::Optimizer, x::MOI.VariableIndex)
    return MOI.is_valid(model.inner, x)
end

MOI.add_variable(model::Optimizer) = MOI.add_variable(model.inner)

function MOI.delete(model::Optimizer, x::MOI.VariableIndex)
    MOI.delete(model.inner, x)
    return
end

function MOI.supports(
    model::Optimizer,
    arg::MOI.AbstractVariableAttribute,
    ::Type{MOI.VariableIndex},
)
    return MOI.supports(model.inner, arg, MOI.VariableIndex)
end

function MOI.set(
    model::Optimizer,
    attr::MOI.AbstractVariableAttribute,
    indices::Vector{<:MOI.VariableIndex},
    args::Vector{T},
) where {T}
    MOI.set.(model, attr, indices, args)
    return
end

### AbstractConstraintAttribute

function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex)
    return MOI.is_valid(model.inner, ci)
end

function MOI.supports(
    model::Optimizer,
    arg::MOI.AbstractConstraintAttribute,
    ::Type{MOI.ConstraintIndex{F,S}},
) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
    return MOI.supports(model.inner, arg, MOI.ConstraintIndex{F,S})
end

function MOI.set(
    model::Optimizer,
    attr::MOI.AbstractConstraintAttribute,
    indices::Vector{<:MOI.ConstraintIndex},
    args::Vector{T},
) where {T}
    MOI.set.(model, attr, indices, args)
    return
end

function MOI.get(model::Optimizer, ::Type{MOI.VariableIndex}, args...)
    return MOI.get(model.inner, MOI.VariableIndex, args...)
end

function MOI.get(model::Optimizer, T::Type{<:MOI.ConstraintIndex}, args...)
    return MOI.get(model.inner, T, args...)
end

function MOI.delete(model::Optimizer, ci::MOI.ConstraintIndex)
    MOI.delete(model.inner, ci)
    return
end

function MOI.supports_constraint(
    model::Optimizer,
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    return MOI.supports_constraint(model.inner, F, S)
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.AbstractFunction,
    s::MOI.AbstractSet,
)
    return MOI.add_constraint(model.inner, f, s)
end

# Decomposition

decompose(f, s) = [(f, s)]

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VectorAffineFunction{T},
    s::MOI.PositiveSemidefiniteConeTriangle,
) where {T}
    cis = CI[]
    for (fi, si) in decompose(f, s)
        push!(cis, MOI.add_constraint(model.inner, fi, si))
    end
    model.outer_to_inner[first(cis)] = cis
    return first(cis)
end

function MOI.get(
    model::Optimizer,
    ::MOI.NumberOfConstraints{F,S},
) where {
    F<:MOI.VectorAffineFunction{Float64},
    S<:MOI.PositiveSemidefiniteConeTriangle,
}
    return length(model.outer_to_inner)
end

function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{F,S},
) where {
    F<:MOI.VectorAffineFunction{Float64},
    S<:MOI.PositiveSemidefiniteConeTriangle,
}
    return sort!(collect(keys(model.outer_to_inner)); by = c -> c.value)
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintPrimal,
    ::MOI.ConstraintIndex{F,S},
) where {
    F<:MOI.VectorAffineFunction{Float64},
    S<:MOI.PositiveSemidefiniteConeTriangle,
}
    return # ???
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintDual,
    ::MOI.ConstraintIndex{F,S},
) where {
    F<:MOI.VectorAffineFunction{Float64},
    S<:MOI.PositiveSemidefiniteConeTriangle,
}
    return # ???
end

end  # module

using JuMP
import .MathOptChordalDecomposition
import LinearAlgebra
import SCS
model = Model(() -> MathOptChordalDecomposition.Optimizer(SCS.Optimizer))
@variable(model, X[1:4, 1:4], Symmetric)
A = [1 1 0 0; 1 1 0 0; 0 0 1 1; 0 0 1 1]
@constraint(model, LinearAlgebra.Symmetric(A .* X) in PSDCone())
optimize!(model)