diff --git a/Project.toml b/Project.toml index b2ea29b..fc4f6ee 100644 --- a/Project.toml +++ b/Project.toml @@ -11,7 +11,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] MATLAB = "0.7.3 - 0.9" -MathOptInterface = "1.1" +MathOptInterface = "1.40" julia = "1.6" [extras] diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 7fd79df..931412a 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -12,7 +12,7 @@ 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 @@ -22,6 +22,20 @@ import LinearAlgebra # supported supported sets are `VectorAffineFunction`-in-`S` where `S` is one # of the sets just listed above. +# If these `const` are remove and coopy-pasted in the `struct` below, then we +# get hit by this assertion error: +# https://github.com/jump-dev/MathOptInterface.jl/blob/86bfa1d37d41f95ae8e2b9e6a7436e17ebe14d81/src/Utilities/struct_of_constraints.jl#L254 +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, @@ -35,22 +49,34 @@ 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 @@ -58,12 +84,14 @@ mutable struct Solution end mutable struct Optimizer <: MOI.AbstractOptimizer - cones::Union{Nothing,Cones{Float64}} + cones_real::Union{Nothing,Cones{Float64}} + cones_complex::Union{Nothing,ComplexCones{ComplexF64}} + complex_offset::Int sol::Union{Nothing,Solution} silent::Bool options::Dict{Symbol,Any} function Optimizer() - return new(nothing, nothing, false, Dict{Symbol,Any}()) + return new(nothing, nothing, 0, nothing, false, Dict{Symbol,Any}()) end end @@ -72,16 +100,18 @@ function MOI.default_cache(::Optimizer, ::Type{Float64}) end function MOI.get(::Optimizer, ::MOI.Bridges.ListOfNonstandardBridges) - return Type[ScaledPSDConeBridge{Float64}] + return Type[ScaledPSDConeBridge{Float64}, ScaledPSDConeBridge{ComplexF64}] end MOI.get(::Optimizer, ::MOI.SolverName) = "SeDuMi" function MOI.is_empty(optimizer::Optimizer) - return optimizer.cones === nothing + return optimizer.cones_real === nothing && + optimizer.cones_complex === nothing end function MOI.empty!(optimizer::Optimizer) - optimizer.cones = nothing + optimizer.cones_real = nothing + optimizer.cones_complex = nothing optimizer.sol = nothing return end @@ -141,8 +171,16 @@ 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 @@ -150,13 +188,23 @@ 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()) @@ -167,9 +215,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 @@ -180,20 +243,27 @@ 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], + convert(Vector{Int}, length(realPSDdims) .+ (1: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_real = deepcopy(Ac_real.sets) + dest.cones_complex = deepcopy(Ac_complex.sets) + dest.complex_offset = length(c_real) 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, @@ -293,6 +363,7 @@ function MOI.get(optimizer::Optimizer, attr::MOI.ObjectiveValue) end return value end + function MOI.get(optimizer::Optimizer, attr::MOI.DualObjectiveValue) MOI.check_result_index_bounds(optimizer, attr) value = optimizer.sol.dual_objective_value @@ -335,6 +406,7 @@ function MOI.get( return MOI.NEARLY_FEASIBLE_POINT end end + function MOI.get( optimizer::Optimizer, attr::MOI.VariablePrimal, @@ -343,22 +415,49 @@ function MOI.get( MOI.check_result_index_bounds(optimizer, attr) return optimizer.sol.y[vi.value] end + +function MOI.get( + optimizer::Optimizer, + attr::MOI.ConstraintPrimal, + ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}}, +) + MOI.check_result_index_bounds(optimizer, attr) + return convert( + Vector{Float64}, + optimizer.sol.slack[MOI.Utilities.rows(optimizer.cones_real, ci)], + ) +end + function MOI.get( optimizer::Optimizer, attr::MOI.ConstraintPrimal, - ci::MOI.ConstraintIndex, + ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{ComplexF64}}, ) MOI.check_result_index_bounds(optimizer, attr) - return optimizer.sol.slack[MOI.Utilities.rows(optimizer.cones, ci)] + rows = MOI.Utilities.rows(optimizer.cones_complex, ci) + return optimizer.sol.slack[optimizer.complex_offset .+ rows] +end + +function MOI.get( + optimizer::Optimizer, + attr::MOI.ConstraintDual, + ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}}, +) + MOI.check_result_index_bounds(optimizer, attr) + return convert( + Vector{Float64}, + optimizer.sol.x[MOI.Utilities.rows(optimizer.cones_real, ci)], + ) end function MOI.get( optimizer::Optimizer, attr::MOI.ConstraintDual, - ci::MOI.ConstraintIndex, + ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{ComplexF64}}, ) MOI.check_result_index_bounds(optimizer, attr) - return optimizer.sol.x[MOI.Utilities.rows(optimizer.cones, ci)] + rows = MOI.Utilities.rows(optimizer.cones_complex, ci) + return optimizer.sol.x[optimizer.complex_offset .+ rows] end MOI.get(optimizer::Optimizer, ::MOI.ResultCount) = 1 diff --git a/src/SeDuMi.jl b/src/SeDuMi.jl index 78d19e3..880b8ad 100644 --- a/src/SeDuMi.jl +++ b/src/SeDuMi.jl @@ -22,10 +22,18 @@ 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) - return Cone(f, l, q, r, s, Float64[], Float64[], Float64[]) + +function Cone( + f::Real, + l::Real, + q::Vector = Float64[], + r::Vector = Float64[], + s::Vector = Float64[], + scomplex::Vector = Float64[], +) + return Cone(f, l, q, r, s, scomplex, Float64[], Float64[]) end -Cone(f::Real, l::Real) = Cone(f, l, Float64[], Float64[], Float64[]) + dimension(K::Cone) = K.f + K.l + sum(K.q) + sum(K.r) + sum(K.s .^ 2) to_vec(x::Vector) = x @@ -38,7 +46,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, diff --git a/src/scaled_psd_cone_bridge.jl b/src/scaled_psd_cone_bridge.jl index 6982e2a..24de1fe 100644 --- a/src/scaled_psd_cone_bridge.jl +++ b/src/scaled_psd_cone_bridge.jl @@ -141,7 +141,9 @@ end _row(i, t::MOI.VectorAffineTerm) = t.output_index _row(i, β) = i -_prod(α, t::MOI.VectorAffineTerm) = MOI.Utilities.operate_term(*, α, t) +function _prod(α, t::MOI.VectorAffineTerm{T}) where {T} + return MOI.Utilities.operate_term(*, convert(T, α), t) +end _prod(α, β) = α * β # Scale coefficients depending on rows index on symmetric packed upper triangular form diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl index 4def97a..f2b88b4 100644 --- a/test/MOI_wrapper.jl +++ b/test/MOI_wrapper.jl @@ -6,11 +6,9 @@ module TestSeDuMi using Test -using MathOptInterface +import MathOptInterface as MOI import SeDuMi -const MOI = MathOptInterface - function runtests() for name in names(@__MODULE__; all = true) if startswith("$(name)", "test_") @@ -32,6 +30,51 @@ function test_options() @test MOI.get(optimizer, MOI.RawOptimizerAttribute("fid")) == 0 end +function test_complex_bridges() + model = MOI.instantiate(SeDuMi.Optimizer; with_bridge_type = Float64) + ComplexF = MOI.VectorAffineFunction{ComplexF64} + S = SeDuMi.ScaledPSDCone + @test MOI.supports_constraint(model, ComplexF, S) + S = MOI.PositiveSemidefiniteConeTriangle + @test MOI.supports_constraint(model, ComplexF, S) + S = MOI.HermitianPositiveSemidefiniteConeTriangle + F = MOI.VectorAffineFunction{Float64} + @test MOI.supports_constraint(model, F, S) + @test MOI.Bridges.bridge_type(model, F, S) == + MOI.Bridges.Constraint.HermitianToComplexSymmetricBridge{ + Float64, + ComplexF, + F, + } + return +end + +function test_complex() + model = MOI.instantiate(SeDuMi.Optimizer; with_cache_type = Float64) + x = MOI.add_variable(model) + T = ComplexF64 + c = MOI.add_constraint( + model, + MOI.Utilities.vectorize([ + one(T) + zero(T) * x, + -1.0im * x, + 1.0im * x, + one(T) + zero(T) * x, + ]), + SeDuMi.ScaledPSDCone(2), + ) + MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE) + obj = 1.0 * x + MOI.set(model, MOI.ObjectiveFunction{typeof(obj)}(), obj) + MOI.optimize!(model) + @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL + @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 1 rtol = 1e-5 + @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ [1, -im, im, 1] rtol = + 1e-5 + @test MOI.get(model, MOI.ConstraintDual(), c) ≈ [0.5, 0.5im, -0.5im, 0.5] rtol = + 1e-5 +end + function test_runtests() model = MOI.Utilities.CachingOptimizer( MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), @@ -42,7 +85,7 @@ function test_runtests() # `Variable.ZerosBridge` makes dual needed by some tests fail. MOI.Bridges.remove_bridge( model.optimizer, - MathOptInterface.Bridges.Variable.ZerosBridge{Float64}, + MOI.Bridges.Variable.ZerosBridge{Float64}, ) MOI.set(model, MOI.Silent(), true) MOI.Test.runtests(