support complex PSD cone

src/MOI_wrapper.jl

@@ -2,12 +2,13 @@ using MathOptInterface
 const MOI = MathOptInterface
 
 include("scaled_psd_cone_bridge.jl")
+include("scaled_hermitian_psd_cone_bridge.jl")
 
 import LinearAlgebra
 
 # SeDuMi solves the primal/dual pair
 # min c'x,       max b'y
-# s.t. Ax = b,   c - A'x ∈ K
+# s.t. Ax = b,   c - A'y ∈ K
 #       x ∈ K
 # where K is a product of `MOI.Zeros`, `MOI.Nonnegatives`,
 # `MOI.SecondOrderCone`, `MOI.RotatedSecondOrderCone` and
@@ -17,6 +18,17 @@ import LinearAlgebra
 # supported supported sets are `VectorAffineFunction`-in-`S` where `S` is one
 # of the sets just listed above.
 
+const ComplexAff = MOI.VectorAffineFunction{Complex{Float64}}
+const RealAff = MOI.VectorAffineFunction{Float64}
+
+MOI.Utilities.@struct_of_constraints_by_function_types(
+    ComplexOrReal,
+    ComplexAff,
+    RealAff,
+)
+
+MOI.Utilities.@product_of_sets(ComplexCones, ScaledPSDCone,)
+
 MOI.Utilities.@product_of_sets(
     Cones,
     MOI.Zeros,
@@ -30,30 +42,42 @@ const OptimizerCache = MOI.Utilities.GenericModel{
     Float64,
     MOI.Utilities.ObjectiveContainer{Float64},
     MOI.Utilities.VariablesContainer{Float64},
-    MOI.Utilities.MatrixOfConstraints{
-        Float64,
-        MOI.Utilities.MutableSparseMatrixCSC{
+    ComplexOrReal{Float64}{
+        MOI.Utilities.MatrixOfConstraints{
+            ComplexF64,
+            MOI.Utilities.MutableSparseMatrixCSC{
+                ComplexF64,
+                Int,
+                MOI.Utilities.OneBasedIndexing,
+            },
+            Vector{ComplexF64},
+            ComplexCones{ComplexF64},
+        },
+        MOI.Utilities.MatrixOfConstraints{
             Float64,
-            Int,
-            MOI.Utilities.OneBasedIndexing,
+            MOI.Utilities.MutableSparseMatrixCSC{
+                Float64,
+                Int,
+                MOI.Utilities.OneBasedIndexing,
+            },
+            Vector{Float64},
+            Cones{Float64},
         },
-        Vector{Float64},
-        Cones{Float64},
     },
 }
 
 mutable struct Solution
-    x::Vector{Float64}
+    x::Vector{ComplexF64}
     y::Vector{Float64}
-    slack::Vector{Float64}
+    slack::Vector{ComplexF64}
     objective_value::Float64
     dual_objective_value::Float64
     objective_constant::Float64
     info::Dict{String,Any}
 end
 
 mutable struct Optimizer <: MOI.AbstractOptimizer
-    cones::Union{Nothing,Cones{Float64}}
+    cones::Union{Nothing,Cones{Float64},ComplexCones{ComplexF64}}
     sol::Union{Nothing,Solution}
     silent::Bool
     options::Dict{Symbol,Any}
@@ -67,7 +91,10 @@ function MOI.default_cache(::Optimizer, ::Type{Float64})
 end
 
 function MOI.get(::Optimizer, ::MOI.Bridges.ListOfNonstandardBridges)
-    return Type[ScaledPSDConeBridge{Float64}]
+    return Type[
+        ScaledPSDConeBridge{Float64},
+        ScaledHermitianPSDConeBridge{Float64},
+    ]
 end
 
 MOI.get(::Optimizer, ::MOI.SolverName) = "SeDuMi"
@@ -136,22 +163,40 @@ function MOI.supports_constraint(
     return true
 end
 
+function MOI.supports_constraint(
+    ::Optimizer,
+    ::Type{MOI.VectorAffineFunction{ComplexF64}},
+    ::Type{ScaledPSDCone},
+)
+    return true
+end
+
 function _map_sets(f, sets, ::Type{S}) where {S}
-    F = MOI.VectorAffineFunction{Float64}
+    F = MOI.VectorAffineFunction{eltype(sets.coefficients.nzval)}
     cis = MOI.get(sets, MOI.ListOfConstraintIndices{F,S}())
     return Int[f(MOI.get(sets, MOI.ConstraintSet(), ci)) for ci in cis]
 end
 
 function MOI.optimize!(dest::Optimizer, src::OptimizerCache)
     MOI.empty!(dest)
     index_map = MOI.Utilities.identity_index_map(src)
-    Ac = src.constraints
-    A = Ac.coefficients
+    Ac_real = MOI.Utilities.constraints(
+        src.constraints,
+        MOI.VectorAffineFunction{Float64},
+        ScaledPSDCone,
+    )
+    A_real = Ac_real.coefficients
+    Ac_complex = MOI.Utilities.constraints(
+        src.constraints,
+        MOI.VectorAffineFunction{ComplexF64},
+        ScaledPSDCone,
+    )
+    A_complex = Ac_complex.coefficients
 
     model_attributes = MOI.get(src, MOI.ListOfModelAttributesSet())
     objective_function_attr =
         MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()
-    b = zeros(A.n)
+    b = zeros(A_real.n) #must be equal to A_complex.n
     max_sense = MOI.get(src, MOI.ObjectiveSense()) == MOI.MAX_SENSE
     objective_constant = 0.0
     if objective_function_attr in MOI.get(src, MOI.ListOfModelAttributesSet())
@@ -162,9 +207,24 @@ function MOI.optimize!(dest::Optimizer, src::OptimizerCache)
         end
     end
 
+    AR = SparseMatrixCSC(
+        A_real.m,
+        A_real.n,
+        A_real.colptr,
+        A_real.rowval,
+        -A_real.nzval,
+    )
+    AC = SparseMatrixCSC(
+        A_complex.m,
+        A_complex.n,
+        A_complex.colptr,
+        A_complex.rowval,
+        -A_complex.nzval,
+    )
+    A = A_complex.m == 0 ? AR : [AR; AC] #keep it real if we can
+
     # If m == n, SeDuMi thinks we give A'.
     # See https://github.com/sqlp/sedumi/issues/42#issuecomment-451300096
-    A = SparseMatrixCSC(A.m, A.n, A.colptr, A.rowval, -A.nzval)
     if A.m != A.n
         A = SparseMatrixCSC(A')
     end
@@ -175,20 +235,28 @@ function MOI.optimize!(dest::Optimizer, src::OptimizerCache)
         options[:fid] = 0
     end
 
+    realPSDdims = _map_sets(MOI.side_dimension, Ac_real, ScaledPSDCone)
+    complexPSDdims = _map_sets(MOI.side_dimension, Ac_complex, ScaledPSDCone)
     K = Cone(
-        Ac.sets.num_rows[1],
-        Ac.sets.num_rows[2] - Ac.sets.num_rows[1],
-        _map_sets(MOI.dimension, Ac, MOI.SecondOrderCone),
-        _map_sets(MOI.dimension, Ac, MOI.RotatedSecondOrderCone),
-        _map_sets(MOI.side_dimension, Ac, ScaledPSDCone),
+        Ac_real.sets.num_rows[1],
+        Ac_real.sets.num_rows[2] - Ac_real.sets.num_rows[1],
+        _map_sets(MOI.dimension, Ac_real, MOI.SecondOrderCone),
+        _map_sets(MOI.dimension, Ac_real, MOI.RotatedSecondOrderCone),
+        [realPSDdims; complexPSDdims],
+        Vector(
+            length(realPSDdims)+1:length(realPSDdims)+length(complexPSDdims),
+        ),
     )
 
-    c = Ac.constants
+    c_real = Ac_real.constants
+    c_complex = Ac_complex.constants
+    c = isempty(c_complex) ? c_real : [c_real; c_complex]
     x, y, info = sedumi(A, b, c, K; options...)
 
-    dest.cones = deepcopy(Ac.sets)
+    dest.cones = deepcopy(Ac_real.sets)
+    #dest.cones = deepcopy(Ac_complex.sets)
     objective_value = (max_sense ? 1 : -1) * LinearAlgebra.dot(b, y)
-    dual_objective_value = (max_sense ? 1 : -1) * LinearAlgebra.dot(c, x)
+    dual_objective_value = (max_sense ? 1 : -1) * real(LinearAlgebra.dot(c, x))
     dest.sol = Solution(
         x,
         y,

src/SeDuMi.jl

@@ -17,6 +17,16 @@ mutable struct Cone
     ycomplex::Vector{Float64} #list of constraints on A*x that should also act on the imaginary part
     xcomplex::Vector{Float64} #list of components of f,l,q,r blocks allowed to be complex
 end
+function Cone(
+    f::Real,
+    l::Real,
+    q::Vector,
+    r::Vector,
+    s::Vector,
+    scomplex::Vector,
+)
+    return Cone(f, l, q, r, s, scomplex, Float64[], Float64[])
+end
 function Cone(f::Real, l::Real, q::Vector, r::Vector, s::Vector)
     return Cone(f, l, q, r, s, Float64[], Float64[], Float64[])
 end
@@ -33,7 +43,7 @@ end
 
 # Solve the primal/dual pair
 # min c'x,      max b'y
-# s.t. Ax = b,   c - A'x ∈ K
+# s.t. Ax = b,   c - A'y ∈ K
 #       x ∈ K
 #
 # Note, if `A` is square then SeDuMi assumes that `A'` is passed instead,

src/scaled_hermitian_psd_cone_bridge.jl

@@ -0,0 +1,124 @@
+struct ScaledHermitianPSDConeBridge{T,G} <: MOI.Bridges.Constraint.SetMapBridge{
+    T,
+    ScaledPSDCone,
+    MOI.HermitianPositiveSemidefiniteConeTriangle,
+    MOI.VectorAffineFunction{ComplexF64},
+    G,
+}
+    constraint::MOI.ConstraintIndex{
+        MOI.VectorAffineFunction{ComplexF64},
+        ScaledPSDCone,
+    }
+end
+
+function MOI.Bridges.Constraint.concrete_bridge_type(
+    ::Type{ScaledHermitianPSDConeBridge{T}},
+    ::Type{G},
+    ::Type{MOI.HermitianPositiveSemidefiniteConeTriangle},
+) where {T,G<:Union{MOI.VectorOfVariables,MOI.VectorAffineFunction{T}}}
+    return ScaledHermitianPSDConeBridge{T,G}
+end
+
+function MOI.Bridges.map_set(
+    ::Type{<:ScaledHermitianPSDConeBridge},
+    set::MOI.HermitianPositiveSemidefiniteConeTriangle,
+)
+    return ScaledPSDCone(set.side_dimension)
+end
+
+function MOI.Bridges.inverse_map_set(
+    ::Type{<:ScaledHermitianPSDConeBridge},
+    set::ScaledPSDCone,
+)
+    return MOI.HermitianPositiveSemidefiniteConeTriangle(set.side_dimension)
+end
+
+# Map ConstraintFunction from MOI -> SeDuMi
+function MOI.Bridges.map_function(
+    BT::Type{<:ScaledHermitianPSDConeBridge{T}},
+    func::MOI.VectorOfVariables,
+) where {T}
+    new_f = MOI.Utilities.operate(*, Float64, 1.0, func)
+    return MOI.Bridges.map_function(BT, new_f)
+end
+function MOI.Bridges.map_function(
+    ::Type{<:ScaledHermitianPSDConeBridge},
+    f::MOI.VectorAffineFunction,
+)
+    n = MOI.output_dimension(f)
+    d = isqrt(n) #side dimension of Hermitian matrix
+    constants = copy(f.constants)
+    constants = hermitian_to_complex_triangle(constants, d)
+    scale_coefficients!(constants)
+    constants = triangle_to_square(constants, d)
+    terms = copy(f.terms)
+    terms = hermitian_to_complex_triangle_indices(terms, d)
+    scale_coefficients!(terms)
+    triangle_to_square_indices!(terms, d)
+    return MOI.VectorAffineFunction(terms, constants)
+end
+
+# Used to map the ConstraintPrimal from SeDuMi -> MOI
+function MOI.Bridges.inverse_map_function(
+    ::Type{<:ScaledHermitianPSDConeBridge},
+    square,
+)
+    return _sedumi_to_moi(square)
+end
+
+# Used to map the ConstraintDual from SeDuMi -> MOI
+function MOI.Bridges.adjoint_map_function(
+    ::Type{<:ScaledHermitianPSDConeBridge},
+    square,
+)
+    return _sedumi_to_moi(square)
+end
+
+function _sedumi_to_moi(square)
+    n = isqrt(length(square))
+    triangle_size = div(n * (n + 1), 2)
+    triangle = Vector{real(eltype(square))}(undef, n^2)
+    for j in 1:n, i in 1:j
+        triangle[MOI.Utilities.trimap(i, j)] = real(square[square_map(i, j, n)])
+    end
+    counter = 0
+    for j in 2:n, i in 1:j-1
+        counter += 1
+        triangle[triangle_size+counter] = imag(square[square_map(i, j, n)])
+    end
+    return triangle
+end
+
+function hermitian_to_complex_triangle(x, n)
+    triangle_size = div(n * (n + 1), 2)
+    y = ComplexF64.(x[1:triangle_size])
+    for i in 1:n-1, j in i+1:n
+        y[MOI.Utilities.trimap(i, j)] +=
+            im * x[triangle_size+MOI.Utilities.trimap(i, j - 1)]
+    end
+    return y
+end
+
+function hermitian_to_complex_triangle_indices(
+    x::Vector{<:MOI.VectorAffineTerm},
+    n,
+)
+    triangle_size = div(n * (n + 1), 2)
+    map = zeros(Int64, div(n * (n - 1), 2))
+    for i in 1:n-1, j in i+1:n
+        map[MOI.Utilities.trimap(i, j - 1)] = MOI.Utilities.trimap(i, j)
+    end
+    x = convert(Vector{MOI.VectorAffineTerm{ComplexF64}}, x)
+    for i in eachindex(x)
+        if x[i].output_index >= triangle_size + 1
+            x[i] = MOI.VectorAffineTerm(
+                map[x[i].output_index-triangle_size],
+                MOI.ScalarAffineTerm(
+                    im * x[i].scalar_term.coefficient,
+                    x[i].scalar_term.variable,
+                ),
+            )
+        end
+    end
+    return x
+end

src/scaled_psd_cone_bridge.jl

@@ -136,7 +136,10 @@ end
 _row(i, t::MOI.VectorAffineTerm) = t.output_index
 _row(i, β) = i
 
-_prod(α, t::MOI.VectorAffineTerm) = MOI.Utilities.operate_term(*, α, t)
+_prod(α, t::MOI.VectorAffineTerm{Float64}) = MOI.Utilities.operate_term(*, α, t)
+function _prod(α, t::MOI.VectorAffineTerm{ComplexF64})
+    return MOI.Utilities.operate_term(*, ComplexF64(α), t)
+end
 _prod(α, β) = α * β
 
 # Scale coefficients depending on rows index on symmetric packed upper triangular form