add MOI interface for complex PSD cone

Project.toml

@@ -5,6 +5,7 @@ version = "2.1.0"
 
 [deps]
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 SCS_jll = "f4f2fc5b-1d94-523c-97ea-2ab488bedf4b"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -19,6 +20,7 @@ SCSSCS_MKL_jllExt = ["SCS_MKL_jll"]
 
 [compat]
 MathOptInterface = "1.20"
+MutableArithmetics = "1"
 Pkg = "<0.0.1, ^1.6"
 PrecompileTools = "1"
 SCS_GPU_jll = "=3.2.8"

src/MOI_wrapper/MOI_wrapper.jl

@@ -4,13 +4,16 @@
 # in the LICENSE.md file or at https://opensource.org/licenses/MIT.
 
 include("scaled_psd_cone_bridge.jl")
+include("hermitian_complex_psd_cone_bridge.jl")
+include("scaled_complex_psd_cone_bridge.jl")
 
 MOI.Utilities.@product_of_sets(
     _Cones,
     MOI.Zeros,
     MOI.Nonnegatives,
     MOI.SecondOrderCone,
     ScaledPSDCone,
+    ScaledComplexPSDCone,
     MOI.ExponentialCone,
     MOI.DualExponentialCone,
     MOI.PowerCone{T},
@@ -132,7 +135,11 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
 end
 
 function MOI.get(::Optimizer, ::MOI.Bridges.ListOfNonstandardBridges)
-    return Type[ScaledPSDConeBridge{Cdouble}]
+    return Type[
+        ScaledPSDConeBridge{Cdouble},
+        ScaledComplexPSDConeBridge{Cdouble},
+        HermitianComplexPSDConeBridge{Cdouble},
+    ]
 end
 
 MOI.get(::Optimizer, ::MOI.SolverName) = "SCS"
@@ -239,6 +246,7 @@ function MOI.supports_constraint(
             MOI.Nonnegatives,
             MOI.SecondOrderCone,
             ScaledPSDCone,
+            ScaledComplexPSDCone,
             MOI.ExponentialCone,
             MOI.DualExponentialCone,
             MOI.PowerCone{Cdouble},
@@ -390,9 +398,9 @@ function MOI.optimize!(
         Float64[], # # placeholder: bl
         _map_sets(MOI.dimension, T, Ab, MOI.SecondOrderCone),
         _map_sets(MOI.side_dimension, T, Ab, ScaledPSDCone),
-        Int[], # placeholder complex PSD
-        div(Ab.sets.num_rows[5] - Ab.sets.num_rows[4], 3),
+        _map_sets(MOI.side_dimension, T, Ab, ScaledComplexPSDCone),
         div(Ab.sets.num_rows[6] - Ab.sets.num_rows[5], 3),
+        div(Ab.sets.num_rows[7] - Ab.sets.num_rows[6], 3),
         vcat(
             _map_sets(set -> set.exponent, Float64, Ab, MOI.PowerCone{Cdouble}),
             _map_sets(

src/MOI_wrapper/hermitian_complex_psd_cone_bridge.jl

@@ -0,0 +1,192 @@
+# Copyright (c) 2014: SCS.jl contributors
+#
+# Use of this source code is governed by an MIT-style license that can be found
+# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
+
+"""
+    struct ComplexPositiveSemidefiniteConeTriangle <: MOI.AbstractVectorSet
+        side_dimension::Int
+    end
+
+Similar to `HermitianPositiveSemidefiniteConeTriangle` but its
+vectorization interleaves real and imaginary parts of the off-diagonal.
+"""
+struct ComplexPositiveSemidefiniteConeTriangle <: MOI.AbstractVectorSet
+    side_dimension::Int
+end
+
+function MOI.Utilities.set_with_dimension(
+    ::Type{ComplexPositiveSemidefiniteConeTriangle},
+    dim,
+)
+    return ComplexPositiveSemidefiniteConeTriangle(isqrt(dim))
+end
+
+function Base.copy(x::ComplexPositiveSemidefiniteConeTriangle)
+    return ComplexPositiveSemidefiniteConeTriangle(x.side_dimension)
+end
+
+function MOI.side_dimension(x::ComplexPositiveSemidefiniteConeTriangle)
+    return x.side_dimension
+end
+
+function MOI.dimension(x::ComplexPositiveSemidefiniteConeTriangle)
+    return x.side_dimension^2
+end
+
+import MutableArithmetics as MA
+
+function MOI.Utilities.set_dot(
+    x::AbstractVector{S},
+    y::AbstractVector{T},
+    set::ComplexPositiveSemidefiniteConeTriangle,
+) where {S,T}
+    U = MA.promote_operation(MA.add_mul, S, T)
+    result = zero(U)
+    d = set.side_dimension
+    k = 0
+    for j in 1:d
+        for i in 1:j-1
+            k += 1
+            result = MA.add_mul!!(result, 2, x[k], y[k])
+            k += 1
+            result = MA.add_mul!!(result, 2, x[k], y[k])
+        end
+        k += 1
+        result = MA.add_mul!!(result, x[k], y[k])
+    end
+    return result
+end
+
+function MOI.Utilities.dot_coefficients(
+    a::AbstractVector,
+    set::ComplexPositiveSemidefiniteConeTriangle,
+)
+    d = set.side_dimension
+    b = copy(a)
+    k = 0
+    for j in 1:d
+        for i in 1:j-1
+            k += 1
+            b[k] /= 2
+            k += 1
+            b[k] /= 2
+        end
+        k += 1
+    end
+    return b
+end
+
+function MOI.is_set_dot_scaled(::Type{ComplexPositiveSemidefiniteConeTriangle})
+    return true
+end
+
+struct HermitianComplexPSDConeBridge{T,F} <:
+       MOI.Bridges.Constraint.SetMapBridge{
+    T,
+    ComplexPositiveSemidefiniteConeTriangle,
+    MOI.HermitianPositiveSemidefiniteConeTriangle,
+    F,
+    F,
+}
+    constraint::MOI.ConstraintIndex{F,ComplexPositiveSemidefiniteConeTriangle}
+end
+
+MOI.Bridges.bridging_cost(::Type{<:HermitianComplexPSDConeBridge}) = 0.1
+
+function MOI.Bridges.Constraint.concrete_bridge_type(
+    ::Type{HermitianComplexPSDConeBridge{T}},
+    ::Type{F},
+    ::Type{MOI.HermitianPositiveSemidefiniteConeTriangle},
+) where {T,F<:MOI.AbstractVectorFunction}
+    return HermitianComplexPSDConeBridge{T,F}
+end
+
+function MOI.Bridges.map_set(
+    ::Type{<:HermitianComplexPSDConeBridge},
+    set::MOI.HermitianPositiveSemidefiniteConeTriangle,
+)
+    return ComplexPositiveSemidefiniteConeTriangle(MOI.side_dimension(set))
+end
+
+function MOI.Bridges.inverse_map_set(
+    ::Type{<:HermitianComplexPSDConeBridge},
+    set::ComplexPositiveSemidefiniteConeTriangle,
+)
+    return MOI.HermitianPositiveSemidefiniteConeTriangle(set.side_dimension)
+end
+
+function _hermitian_to_complex(func)
+    vals = MOI.Utilities.eachscalar(func)
+    dim = length(vals)
+    side = isqrt(dim)
+    k_re = 1
+    k_im = div(side * (side + 1), 2) + 1
+    l = 1
+    perm = zeros(Int, dim)
+    for i in 1:side
+        for j in 1:(i-1)
+            perm[l] = k_re
+            perm[l+1] = k_im
+            k_re += 1
+            k_im += 1
+            l += 2
+        end
+        perm[l] = k_re
+        k_re += 1
+        l += 1
+    end
+    return vals[perm]
+end
+
+function _complex_to_hermitian(func)
+    vals = MOI.Utilities.eachscalar(func)
+    dim = length(vals)
+    side = isqrt(dim)
+    k_re = 1
+    k_im = div(side * (side + 1), 2) + 1
+    l = 1
+    perm = zeros(Int, dim)
+    for i in 1:side
+        for j in 1:(i-1)
+            perm[k_re] = l
+            perm[k_im] = l + 1
+            k_re += 1
+            k_im += 1
+            l += 2
+        end
+        perm[k_re] = l
+        k_re += 1
+        l += 1
+    end
+    return vals[perm]
+end
+
+# Map ConstraintFunction from Hermitian -> Complex
+function MOI.Bridges.map_function(::Type{<:HermitianComplexPSDConeBridge}, f)
+    return _hermitian_to_complex(f)
+end
+
+# Used to map the ConstraintPrimal from Complex -> Hermitian
+function MOI.Bridges.inverse_map_function(
+    ::Type{<:HermitianComplexPSDConeBridge},
+    f,
+)
+    return _complex_to_hermitian(f)
+end
+
+# Used to map the ConstraintDual from Complex -> Hermitian
+function MOI.Bridges.adjoint_map_function(
+    ::Type{<:HermitianComplexPSDConeBridge},
+    f,
+)
+    return _complex_to_hermitian(f)
+end
+
+# Used to set ConstraintDualStart
+function MOI.Bridges.inverse_adjoint_map_function(
+    ::Type{<:HermitianComplexPSDConeBridge},
+    f,
+)
+    return _hermitian_to_complex(f)
+end

src/MOI_wrapper/scaled_complex_psd_cone_bridge.jl

@@ -0,0 +1,131 @@
+# Copyright (c) 2014: SCS.jl contributors
+#
+# Use of this source code is governed by an MIT-style license that can be found
+# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
+
+"""
+    struct ScaledComplexPSDCone <: MOI.AbstractVectorSet
+        side_dimension::Int
+    end
+
+Similar to `MOI.Scaled{ComplexPositiveSemidefiniteConeTriangle}` but its
+vectorization is the lower triangular part column-wise.
+"""
+struct ScaledComplexPSDCone <: MOI.AbstractVectorSet
+    side_dimension::Int
+end
+
+function MOI.Utilities.set_with_dimension(::Type{ScaledComplexPSDCone}, dim)
+    return ScaledComplexPSDCone(isqrt(dim))
+end
+
+Base.copy(x::ScaledComplexPSDCone) = ScaledComplexPSDCone(x.side_dimension)
+
+MOI.side_dimension(x::ScaledComplexPSDCone) = x.side_dimension
+
+function MOI.dimension(x::ScaledComplexPSDCone)
+    return x.side_dimension^2
+end
+
+struct ScaledComplexPSDConeBridge{T,F} <: MOI.Bridges.Constraint.SetMapBridge{
+    T,
+    ScaledComplexPSDCone,
+    MOI.Scaled{ComplexPositiveSemidefiniteConeTriangle},
+    F,
+    F,
+}
+    constraint::MOI.ConstraintIndex{F,ScaledComplexPSDCone}
+end
+
+MOI.Bridges.bridging_cost(::Type{<:ScaledComplexPSDConeBridge}) = 0.1
+
+function MOI.Bridges.Constraint.concrete_bridge_type(
+    ::Type{ScaledComplexPSDConeBridge{T}},
+    ::Type{F},
+    ::Type{MOI.Scaled{ComplexPositiveSemidefiniteConeTriangle}},
+) where {T,F<:MOI.AbstractVectorFunction}
+    return ScaledComplexPSDConeBridge{T,F}
+end
+
+function MOI.Bridges.map_set(
+    ::Type{<:ScaledComplexPSDConeBridge},
+    set::MOI.Scaled{ComplexPositiveSemidefiniteConeTriangle},
+)
+    return ScaledComplexPSDCone(MOI.side_dimension(set))
+end
+
+function MOI.Bridges.inverse_map_set(
+    ::Type{<:ScaledComplexPSDConeBridge},
+    set::ScaledComplexPSDCone,
+)
+    return MOI.Scaled(
+        ComplexPositiveSemidefiniteConeTriangle(set.side_dimension),
+    )
+end
+
+function _complex_to_scs(func)
+    vals = MOI.Utilities.eachscalar(func)
+    d = isqrt(length(vals))
+    c = 0
+    perm = zeros(Int, length(vals))
+    for i in 1:d
+        c += 1
+        perm[c] = i^2
+        for j in i+1:d
+            triidx = 2i - 1 + (j - 1)^2
+            c += 1
+            perm[c] = triidx
+            c += 1
+            perm[c] = triidx + 1
+        end
+    end
+    return vals[perm]
+end
+
+function _scs_to_complex(func)
+    vals = MOI.Utilities.eachscalar(func)
+    d = isqrt(length(vals))
+    c = 0
+    perm = zeros(Int, length(vals))
+    for i in 1:d
+        c += 1
+        perm[i^2] = c
+        for j in i+1:d
+            triidx = 2i - 1 + (j - 1)^2
+            c += 1
+            perm[triidx] = c
+            c += 1
+            perm[triidx+1] = c
+        end
+    end
+    return vals[perm]
+end
+
+# Map ConstraintFunction from MOI -> SCS
+function MOI.Bridges.map_function(::Type{<:ScaledComplexPSDConeBridge}, f)
+    return _complex_to_scs(f)
+end
+
+# Used to map the ConstraintPrimal from SCS -> MOI
+function MOI.Bridges.inverse_map_function(
+    ::Type{<:ScaledComplexPSDConeBridge},
+    f,
+)
+    return _scs_to_complex(f)
+end
+
+# Used to map the ConstraintDual from SCS -> MOI
+function MOI.Bridges.adjoint_map_function(
+    ::Type{<:ScaledComplexPSDConeBridge},
+    f,
+)
+    return _scs_to_complex(f)
+end
+
+# Used to set ConstraintDualStart
+function MOI.Bridges.inverse_adjoint_map_function(
+    ::Type{<:ScaledComplexPSDConeBridge},
+    f,
+)
+    return _complex_to_scs(f)
+end

test/MOI_wrapper.jl

@@ -61,7 +61,6 @@ function _test_runtests(linear_solver)
             # Unexpected failures:
             #   TODO(odow): looks like a tolerance issue?
             "test_linear_add_constraints",
-            "test_conic_HermitianPositiveSemidefiniteConeTriangle_2",
             # Expected test failures:
             #   TODO(odow): get not supported for primal/dual starts
             "test_model_ModelFilter_AbstractConstraintAttribute",