add automated extended formulations for separable spectral cone examples

examples/Examples.jl

@@ -8,7 +8,7 @@ import DataStructures: OrderedDict
 include("common.jl")
 include("common_native.jl")
 include("common_JuMP.jl")
-
+include("extendedforms_JuMP.jl")
 include("setup.jl")
 
 const model_types = [

examples/classicalquantum/JuMP.jl

@@ -7,6 +7,7 @@ and listing 1 in "Efficient optimization of the quantum relative entropy" by H.
 struct ClassicalQuantum{T <: Real} <: ExampleInstanceJuMP{T}
     n::Int
     complex::Bool
+    use_standard_cones::Bool
 end
 
 function build(inst::ClassicalQuantum{T}) where {T <: Float64}
@@ -28,15 +29,14 @@ function build(inst::ClassicalQuantum{T}) where {T <: Float64}
     JuMP.@variable(model, qe_epi)
     JuMP.@objective(model, Max, -qe_epi + dot(prob, Hs))
 
-    cone = Hypatia.EpiPerSepSpectralCone{T}(Hypatia.Cones.NegEntropySSF(),
-        Cones.MatrixCSqr{T, R}, n)
-    entr_sum_vec = zeros(JuMP.AffExpr, MOI.dimension(cone) - 2)
+    entr_sum_vec = zeros(JuMP.AffExpr, Cones.svec_length(R, n))
     ρ_vec = zeros(T, length(entr_sum_vec))
     for (ρ, p) in zip(ρs, prob)
         Cones.smat_to_svec!(ρ_vec, ρ, rt2)
         entr_sum_vec += p * ρ_vec
     end
-    JuMP.@constraint(model, vcat(qe_epi, 1, entr_sum_vec) in cone)
+    add_sepspectral(Cones.NegEntropySSF(), Cones.MatrixCSqr{T, R}, n,
+        vcat(qe_epi, 1, entr_sum_vec), model, inst.use_standard_cones)
 
     return model
 end

examples/classicalquantum/JuMP_test.jl

@@ -1,8 +1,9 @@
 
 insts = OrderedDict()
 insts["minimal"] = [
-    ((2, false),),
-    ((2, true),),
+    ((3, false, false),),
+    ((3, false, true),),
+    ((3, true, false),),
     ]
 insts["fast"] = [
     ((10, false),),

examples/covarianceest/JuMP.jl

@@ -12,18 +12,8 @@ where f is a convex spectral function
 
 struct CovarianceEstJuMP{T <: Real} <: ExampleInstanceJuMP{T}
     d::Int
-    ssf::Cones.SepSpectralFun
-end
-
-function CovarianceEstJuMP{T}(d::Int, ssf_name::Symbol) where {T <: Real}
-    ssf_dict = Dict(
-        :InvSSF => Cones.InvSSF(),
-        :NegLogSSF => Cones.NegLogSSF(),
-        :NegEntropySSF => Cones.NegEntropySSF(),
-        :Power12SSF => Cones.Power12SSF(1.5),
-        )
-    ssf = ssf_dict[ssf_name]
-    return CovarianceEstJuMP{T}(d, ssf)
+    ssf::Symbol
+    use_standard_cones::Bool
 end
 
 function build(inst::CovarianceEstJuMP{T}) where {T <: Float64}
@@ -43,10 +33,17 @@ function build(inst::CovarianceEstJuMP{T}) where {T <: Float64}
     Cones.smat_to_svec!(p_vec, one(T) * p, sqrt(T(2)))
 
     # convex objective
-    f_cone = Hypatia.EpiPerSepSpectralCone{T}(inst.ssf, Cones.MatrixCSqr{T, T}, d)
     JuMP.@variable(model, epi)
     JuMP.@objective(model, Min, epi)
-    JuMP.@constraint(model, vcat(epi, 1, p_vec) in f_cone)
+
+    ssf_dict = Dict(
+        :InvSSF => Cones.InvSSF(),
+        :NegLogSSF => Cones.NegLogSSF(),
+        :NegEntropySSF => Cones.NegEntropySSF(),
+        :Power12SSF => Cones.Power12SSF(1.5),
+        )
+    add_sepspectral(ssf_dict[inst.ssf], Cones.MatrixCSqr{T, T}, d,
+        vcat(epi, 1, p_vec), model, inst.use_standard_cones)
 
     # linear prior constraints
     lin_dim = round(Int, sqrt(d - 1))

examples/covarianceest/JuMP_test.jl

@@ -1,8 +1,10 @@
 
 insts = OrderedDict()
 insts["minimal"] = [
-    ((2, :InvSSF),),
-    ((3, :NegEntropySSF),),
+    ((2, :NegLogSSF, false),),
+    ((2, :NegLogSSF, true),),
+    ((3, :Power12SSF, false),),
+    ((3, :Power12SSF, true),),
     ]
 insts["fast"] = [
     ((10, :InvSSF),),

examples/discretemaxlikelihood/JuMP.jl

@@ -6,6 +6,7 @@ prior
 
 struct DiscreteMaxLikelihood{T <: Real} <: ExampleInstanceJuMP{T}
     d::Int
+    use_standard_cones::Bool
 end
 
 function build(inst::DiscreteMaxLikelihood{T}) where {T <: Float64}
@@ -14,17 +15,17 @@ function build(inst::DiscreteMaxLikelihood{T}) where {T <: Float64}
     freq = rand(1:(2 * d), d)
 
     model = JuMP.Model()
-    JuMP.@variables(model, begin
-        p[1:d]
-        prodp
-    end)
+    JuMP.@variable(model, p[1:d])
+    JuMP.@variable(model, hypo)
+    JuMP.@objective(model, Max, hypo)
+
     JuMP.@constraints(model, begin
-        vcat(prodp, p) in Hypatia.HypoPowerMeanCone{T}(freq / sum(freq))
-        vcat(inv(d), inv(d), p) in Hypatia.EpiPerSepSpectralCone{T}(
-            Cones.NegEntropySSF(), Cones.VectorCSqr{T}, d)
+        vcat(hypo, p) in Hypatia.HypoPowerMeanCone{T}(freq / sum(freq))
         sum(p) == 1
     end)
-    JuMP.@objective(model, Max, prodp)
+
+    add_sepspectral(Cones.NegEntropySSF(), Cones.VectorCSqr{T}, d,
+        vcat(inv(d), inv(d), p), model, inst.use_standard_cones)
 
     return model
 end

examples/discretemaxlikelihood/JuMP_test.jl

@@ -1,7 +1,8 @@
 
 insts = OrderedDict()
 insts["minimal"] = [
-    ((2,),),
+    ((2, false),),
+    ((2, true),),
     ]
 insts["fast"] = [
     ((10,),),

examples/entanglementassisted/JuMP.jl

@@ -44,7 +44,7 @@ function build(inst::EntanglementAssisted{T}) where {T <: Float64}
     Q2_vec = Cones.smat_to_svec!(zeros(JuMP.AffExpr, sa), Q2, rt2)
     Q3_vec = Cones.smat_to_svec!(zeros(JuMP.AffExpr, sb), Q3, rt2)
     RE_cone = Hypatia.EpiTrRelEntropyTriCone{T}(1 + 2 * sa)
-    E_cone = Hypatia.EpiPerSepSpectralCone{T}(Hypatia.Cones.NegEntropySSF(),
+    E_cone = Hypatia.EpiPerSepSpectralCone{T}(Cones.NegEntropySSF(),
         Cones.MatrixCSqr{T, T}, nb)
 
     JuMP.@constraints(model, begin

examples/experimentdesign/JuMP.jl

@@ -16,18 +16,8 @@ experiment (let q ≈ p/2), and f is a convex spectral function
 
 struct ExperimentDesignJuMP{T <: Real} <: ExampleInstanceJuMP{T}
     p::Int
-    ssf::Cones.SepSpectralFun
-end
-
-function ExperimentDesignJuMP{T}(p::Int, ssf_name::Symbol) where {T <: Real}
-    ssf_dict = Dict(
-        :InvSSF => Cones.InvSSF(),
-        :NegLogSSF => Cones.NegLogSSF(),
-        :NegEntropySSF => Cones.NegEntropySSF(),
-        :Power12SSF => Cones.Power12SSF(1.5),
-        )
-    ssf = ssf_dict[ssf_name]
-    return ExperimentDesignJuMP{T}(p, ssf)
+    ssf::Symbol
+    use_standard_cones::Bool
 end
 
 function build(inst::ExperimentDesignJuMP{T}) where {T <: Float64}
@@ -50,8 +40,15 @@ function build(inst::ExperimentDesignJuMP{T}) where {T <: Float64}
     Q = V * diagm(x) * V' # information matrix
     Q_vec = zeros(JuMP.AffExpr, vec_dim)
     Cones.smat_to_svec!(Q_vec, Q, sqrt(T(2)))
-    f_cone = Hypatia.EpiPerSepSpectralCone{T}(inst.ssf, Cones.MatrixCSqr{T, T}, q)
-    JuMP.@constraint(model, vcat(epi, 1, Q_vec) in f_cone)
+
+    ssf_dict = Dict(
+        :InvSSF => Cones.InvSSF(),
+        :NegLogSSF => Cones.NegLogSSF(),
+        :NegEntropySSF => Cones.NegEntropySSF(),
+        :Power12SSF => Cones.Power12SSF(1.5),
+        )
+    add_sepspectral(ssf_dict[inst.ssf], Cones.MatrixCSqr{T, T}, q,
+        vcat(epi, 1, Q_vec), model, inst.use_standard_cones)
 
     return model
 end

examples/experimentdesign/JuMP_test.jl

@@ -1,8 +1,10 @@
 
 insts = OrderedDict()
 insts["minimal"] = [
-    ((2, :InvSSF),),
-    ((3, :NegEntropySSF),),
+    ((2, :InvSSF, false),),
+    ((2, :InvSSF, true),),
+    ((3, :NegEntropySSF, false),),
+    ((3, :NegEntropySSF, true),),
     ]
 insts["fast"] = [
     ((10, :InvSSF),),

examples/extendedforms_JuMP.jl

@@ -0,0 +1,153 @@
+#=
+helpers for constructing extended formulations in JuMP
+=#
+
+# add a separable spectral cone constraint, possibly extended
+function add_sepspectral(
+    h::Cones.SepSpectralFun,
+    csqr::Type{<:Cones.ConeOfSquares{T}},
+    d::Int,
+    aff::Vector{JuMP.AffExpr},
+    model::JuMP.Model,
+    use_standard_cones::Bool,
+    ) where {T <: Real}
+    if use_standard_cones
+        extend_sepspectral(h, csqr, d, aff, model)
+    else
+        JuMP.@constraint(model, aff in
+            Hypatia.EpiPerSepSpectralCone{T}(h, csqr, d))
+    end
+    return
+end
+
+# construct a standard cone EF for a separable spectral vector cone constraint
+function extend_sepspectral(
+    h::Cones.SepSpectralFun,
+    ::Type{Cones.VectorCSqr{T}},
+    d::Int,
+    aff::Vector{JuMP.AffExpr},
+    model::JuMP.Model,
+    ) where {T <: Real}
+    @assert d >= 1
+    @assert d == length(aff) - 2
+
+    # perspective nonnegative NOTE could be redundant
+    JuMP.@constraint(model, aff[2] >= 0)
+
+    # linear constraint
+    x = JuMP.@variable(model, [1:d])
+    JuMP.@constraint(model, aff[1] >= sum(x))
+
+    # 3-dim constraints
+    for i in 1:d
+        aff_i = vcat(x[i], aff[2], aff[2 + i])
+        extend_atom_jump(h, aff_i, model)
+    end
+
+    return
+end
+
+# construct a standard cone EF for a separable spectral matrix cone constraint
+# TODO support complex MatrixCSqr
+function extend_sepspectral(
+    h::Cones.SepSpectralFun,
+    ::Type{Cones.MatrixCSqr{T, T}},
+    d::Int,
+    aff::Vector{JuMP.AffExpr},
+    model::JuMP.Model,
+    ) where {T <: Real}
+    @assert d >= 1
+    @assert Cones.svec_length(d) == length(aff) - 2
+
+    # symmetric matrix (use upper triangle) of aff[2:end]
+    W = zeros(JuMP.AffExpr, d, d)
+    Cones.svec_to_smat!(W, aff[3:end], sqrt(2))
+
+    # eigenvalue ordering
+    λ = JuMP.@variable(model, [1:d])
+    JuMP.@constraint(model, [i in 1:(d - 1)], λ[i] >= λ[i + 1])
+    JuMP.@constraint(model, tr(W) == sum(λ))
+
+    # PSD constraints
+    for i in 1:(d - 1)
+        Z_i = JuMP.@variable(model, [1:d, 1:d], PSD)
+        s_i = JuMP.@variable(model)
+        JuMP.@constraint(model, sum(λ[1:i]) - i * s_i - tr(Z_i) >= 0)
+        mat_i = Symmetric(Z_i - W + s_i * Matrix(I, d, d), :U)
+        JuMP.@SDconstraint(model, mat_i >= 0)
+    end
+
+    # vector separable spectral constraint
+    vec_aff = vcat(aff[1], aff[2], λ)
+    extend_sepspectral(h, Cones.VectorCSqr{T}, d, vec_aff, model)
+
+    return
+end
+
+#=
+InvSSF:
+(x > 0, y > 0, z > 0) : x > y * inv(z / y) = y^2 / z
+↔ x * z > y^2
+↔ (x, z, sqrt(2) * y) ∈ JuMP.RotatedSecondOrderCone
+=#
+function extend_atom_jump(
+    ::Cones.InvSSF,
+    aff::Vector{JuMP.AffExpr},
+    model::JuMP.Model,
+    )
+    @assert length(aff) == 3
+    aff_new = vcat(aff[1], aff[3], sqrt(2) * aff[2])
+    JuMP.@constraint(model, aff_new in JuMP.RotatedSecondOrderCone())
+    return
+end
+
+#=
+NegLogSSF:
+(x, y > 0, z > 0) : x > y * -log(z / y)
+↔ y * exp(-x / y) < z
+↔ (-x, y, z) ∈ MOI.ExponentialCone
+=#
+function extend_atom_jump(
+    ::Cones.NegLogSSF,
+    aff::Vector{JuMP.AffExpr},
+    model::JuMP.Model,
+    )
+    @assert length(aff) == 3
+    aff_new = vcat(-aff[1], aff[2], aff[3])
+    JuMP.@constraint(model, aff_new in MOI.ExponentialCone())
+    return
+end
+
+#=
+NegEntropySSF:
+(x, y > 0, z > 0) : x > z * log(z / y) = -z * log(y / z)
+↔ z * exp(-x / z) < y
+↔ (-x, z, y) ∈ MOI.ExponentialCone
+=#
+function extend_atom_jump(
+    ::Cones.NegEntropySSF,
+    aff::Vector{JuMP.AffExpr},
+    model::JuMP.Model,
+    )
+    @assert length(aff) == 3
+    aff_new = vcat(-aff[1], aff[3], aff[2])
+    JuMP.@constraint(model, aff_new in MOI.ExponentialCone())
+    return
+end
+
+#=
+Power12SSF(p) for p ∈ (1, 2]
+(x > 0, y > 0, z > 0) : x > y * (z / y)^p = y^(1-p) * z^p
+↔ x^(1/p) * y^(1-1/p) > |z|, z > 0
+↔ (x, y, z) ∈ MOI.PowerCone(1/p), z ∈ ℝ₊
+=#
+function extend_atom_jump(
+    h::Cones.Power12SSF,
+    aff::Vector{JuMP.AffExpr},
+    model::JuMP.Model,
+    )
+    @assert length(aff) == 3
+    JuMP.@constraint(model, aff[3] >= 0)
+    JuMP.@constraint(model, aff in MOI.PowerCone(inv(h.p)))
+    return
+end