diff --git a/docs/src/background/duality.md b/docs/src/background/duality.md index 1042e2b5eb..f1c56d27f5 100644 --- a/docs/src/background/duality.md +++ b/docs/src/background/duality.md @@ -113,7 +113,9 @@ and similarly, the dual is: The scalar product is different from the canonical one for the sets [`PositiveSemidefiniteConeTriangle`](@ref), [`LogDetConeTriangle`](@ref), -[`RootDetConeTriangle`](@ref). +[`RootDetConeTriangle`](@ref), +[`SetDotProducts`](@ref) and +[`LinearCombinationInSet`](@ref). If the set ``C_i`` of the section [Duality](@ref) is one of these three cones, then the rows of the matrix ``A_i`` corresponding to off-diagonal entries are diff --git a/docs/src/manual/standard_form.md b/docs/src/manual/standard_form.md index 57b8a8e2d2..4596021400 100644 --- a/docs/src/manual/standard_form.md +++ b/docs/src/manual/standard_form.md @@ -82,6 +82,8 @@ The vector-valued set types implemented in MathOptInterface.jl are: | [`RelativeEntropyCone(d)`](@ref MathOptInterface.RelativeEntropyCone) | ``\{ (u, v, w) \in \mathbb{R}^{d} : u \ge \sum_i w_i \log (\frac{w_i}{v_i}), v_i \ge 0, w_i \ge 0 \}`` | | [`HyperRectangle(l, u)`](@ref MathOptInterface.HyperRectangle) | ``\{x \in \bar{\mathbb{R}}^d: x_i \in [l_i, u_i] \forall i=1,\ldots,d\}`` | | [`NormCone(p, d)`](@ref MathOptInterface.NormCone) | ``\{ (t,x) \in \mathbb{R}^{d} : t \ge \left(\sum\limits_i \lvert x_i \rvert^p\right)^{\frac{1}{p}} \}`` | +| [`SetDotProducts(s, v)`](@ref MathOptInterface.SetDotProducts) | The cone `s` with dot products with the fixed vectors `v`. | +| [`LinearCombinationInSet(s, v)`](@ref MathOptInterface.LinearCombinationInSet) | The cone of vector `(y, x)` such that ``\sum_i y_i v_i + x`` belongs to `s`. | ## Matrix cones diff --git a/docs/src/reference/standard_form.md b/docs/src/reference/standard_form.md index 98ff739a33..3c62a6b632 100644 --- a/docs/src/reference/standard_form.md +++ b/docs/src/reference/standard_form.md @@ -153,4 +153,6 @@ LogDetConeTriangle LogDetConeSquare RootDetConeTriangle RootDetConeSquare +SetDotProducts +LinearCombinationInSet ``` diff --git a/src/Bridges/Constraint/Constraint.jl b/src/Bridges/Constraint/Constraint.jl index 6f6e7bec07..10727ff5be 100644 --- a/src/Bridges/Constraint/Constraint.jl +++ b/src/Bridges/Constraint/Constraint.jl @@ -66,6 +66,7 @@ include("bridges/zero_one.jl") include("bridges/sos1_to_milp.jl") include("bridges/sos2_to_milp.jl") include("bridges/indicator_to_milp.jl") +include("bridges/linear_combination.jl") """ add_all_bridges(bridged_model, ::Type{T}) where {T} @@ -153,6 +154,7 @@ function add_all_bridges(bridged_model, ::Type{T}) where {T} MOI.Bridges.add_bridge(bridged_model, SOS1ToMILPBridge{T}) MOI.Bridges.add_bridge(bridged_model, SOS2ToMILPBridge{T}) MOI.Bridges.add_bridge(bridged_model, IndicatorToMILPBridge{T}) + MOI.Bridges.add_bridge(bridged_model, LinearCombinationBridge{T}) return end diff --git a/src/Bridges/Constraint/bridges/linear_combination.jl b/src/Bridges/Constraint/bridges/linear_combination.jl new file mode 100644 index 0000000000..2656cea235 --- /dev/null +++ b/src/Bridges/Constraint/bridges/linear_combination.jl @@ -0,0 +1,70 @@ +# Copyright (c) 2017: Miles Lubin and contributors +# Copyright (c) 2017: Google Inc. +# +# 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 LinearCombinationBridge{T,S,A,V,F,G} <: + SetMapBridge{T,S,MOI.LinearCombinationInSet{S,A,V},F,G} + constraint::MOI.ConstraintIndex{F,S} + set::MOI.LinearCombinationInSet{S,A,V} +end + +function MOI.supports_constraint( + ::Type{<:LinearCombinationBridge}, + ::Type{<:MOI.AbstractVectorFunction}, + ::Type{<:MOI.LinearCombinationInSet}, +) + return true +end + +function concrete_bridge_type( + ::Type{<:LinearCombinationBridge{T}}, + G::Type{<:MOI.AbstractVectorFunction}, + ::Type{MOI.LinearCombinationInSet{S,A,V}}, +) where {T,S,A,V} + U = MOI.Utilities.promote_operation(*, T, MOI.Utilities.scalar_type(G), T) + F = MOI.Utilities.promote_operation(vcat, T, U) + return LinearCombinationBridge{T,S,A,V,F,G} +end + +function _map_function(set::MOI.LinearCombinationInSet, func) + scalars = MOI.Utilities.eachscalar(func) + return MOI.Utilities.vectorize([ + sum(scalars[j] * set.vectors[j][i] for j in eachindex(scalars)) for + i in 1:MOI.dimension(set.set) + ]) +end + +function bridge_constraint( + ::Type{LinearCombinationBridge{T,S,A,V,F,G}}, + model::MOI.ModelLike, + func::G, + set::MOI.LinearCombinationInSet{S,A,V}, +) where {T,S,A,F,G,V} + mapped_func = _map_function(set, func) + constraint = MOI.add_constraint(model, mapped_func, set.set) + return LinearCombinationBridge{T,S,A,V,F,G}(constraint, set) +end + +function MOI.Bridges.map_set( + ::Type{<:LinearCombinationBridge}, + set::MOI.LinearCombinationInSet, +) + return set.set +end + +function MOI.Bridges.inverse_map_set( + bridge::LinearCombinationBridge, + ::MOI.AbstractSet, +) + return bridge.set +end + +function MOI.Bridges.adjoint_map_function(bridge::LinearCombinationBridge, func) + scalars = MOI.Utilities.eachscalar(func) + return MOI.Utilities.vectorize([ + MOI.Utilities.set_dot(vector, scalars, bridge.set.set) for + vector in bridge.set.vectors + ]) +end diff --git a/src/Bridges/Variable/Variable.jl b/src/Bridges/Variable/Variable.jl index 5c905c9146..38a3b7635c 100644 --- a/src/Bridges/Variable/Variable.jl +++ b/src/Bridges/Variable/Variable.jl @@ -22,6 +22,7 @@ include("bridges/vectorize.jl") include("bridges/zeros.jl") include("bridges/hermitian.jl") include("bridges/parameter.jl") +include("bridges/set_dot.jl") """ add_all_bridges(model, ::Type{T}) where {T} @@ -40,6 +41,7 @@ function add_all_bridges(model, ::Type{T}) where {T} MOI.Bridges.add_bridge(model, RSOCtoPSDBridge{T}) MOI.Bridges.add_bridge(model, HermitianToSymmetricPSDBridge{T}) MOI.Bridges.add_bridge(model, ParameterToEqualToBridge{T}) + MOI.Bridges.add_bridge(model, DotProductsBridge{T}) return end diff --git a/src/Bridges/Variable/bridges/set_dot.jl b/src/Bridges/Variable/bridges/set_dot.jl new file mode 100644 index 0000000000..51cfc3664f --- /dev/null +++ b/src/Bridges/Variable/bridges/set_dot.jl @@ -0,0 +1,82 @@ +# Copyright (c) 2017: Miles Lubin and contributors +# Copyright (c) 2017: Google Inc. +# +# 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 DotProductsBridge{T,S,A,V} <: SetMapBridge{T,S,MOI.SetDotProducts{S,A,V}} + variables::Vector{MOI.VariableIndex} + constraint::MOI.ConstraintIndex{MOI.VectorOfVariables,S} + set::MOI.SetDotProducts{S,A,V} +end + +function supports_constrained_variable( + ::Type{<:DotProductsBridge}, + ::Type{<:MOI.SetDotProducts}, +) + return true +end + +function concrete_bridge_type( + ::Type{<:DotProductsBridge{T}}, + ::Type{MOI.SetDotProducts{S,A,V}}, +) where {T,S,A,V} + return DotProductsBridge{T,S,A,V} +end + +function bridge_constrained_variable( + BT::Type{DotProductsBridge{T,S,A,V}}, + model::MOI.ModelLike, + set::MOI.SetDotProducts{S,A,V}, +) where {T,S,A,V} + variables, constraint = + _add_constrained_var(model, MOI.Bridges.inverse_map_set(BT, set)) + return BT(variables, constraint, set) +end + +function MOI.Bridges.map_set(bridge::DotProductsBridge{T,S}, ::S) where {T,S} + return bridge.set +end + +function MOI.Bridges.inverse_map_set( + ::Type{<:DotProductsBridge}, + set::MOI.SetDotProducts, +) + return set.set +end + +function MOI.Bridges.map_function( + bridge::DotProductsBridge{T}, + func, + i::MOI.Bridges.IndexInVector, +) where {T} + scalars = MOI.Utilities.eachscalar(func) + return MOI.Utilities.set_dot( + bridge.set.vectors[i.value], + scalars, + bridge.set.set, + ) +end + +function MOI.Bridges.map_function(bridge::DotProductsBridge{T}, func) where {T} + scalars = MOI.Utilities.eachscalar(func) + return MOI.Utilities.vectorize([ + MOI.Utilities.set_dot(vector, scalars, bridge.set.set) for + vector in bridge.set.vectors + ]) +end + +# This returns `true` by default for `SetMapBridge` +# but it is not supported for this bridge because `inverse_map_function` +# is not implemented +function MOI.supports( + ::MOI.ModelLike, + ::MOI.VariablePrimalStart, + ::Type{<:DotProductsBridge}, +) + return false +end + +function unbridged_map(::DotProductsBridge, ::Vector{MOI.VariableIndex}) + return nothing +end diff --git a/src/Bridges/Variable/set_map.jl b/src/Bridges/Variable/set_map.jl index 6ebaf57538..6743b8fdb3 100644 --- a/src/Bridges/Variable/set_map.jl +++ b/src/Bridges/Variable/set_map.jl @@ -129,7 +129,7 @@ function MOI.get( bridge::SetMapBridge, ) set = MOI.get(model, attr, bridge.constraint) - return MOI.Bridges.map_set(typeof(bridge), set) + return MOI.Bridges.map_set(bridge, set) end function MOI.set( @@ -149,7 +149,7 @@ function MOI.get( bridge::SetMapBridge, ) value = MOI.get(model, attr, bridge.constraint) - return MOI.Bridges.map_function(typeof(bridge), value) + return MOI.Bridges.map_function(bridge, value) end function MOI.get( @@ -171,7 +171,7 @@ function MOI.get( if any(isnothing, value) return nothing end - return MOI.Bridges.map_function(typeof(bridge), value, i) + return MOI.Bridges.map_function(bridge, value, i) end function MOI.supports( @@ -203,7 +203,7 @@ function MOI.Bridges.bridged_function( i::MOI.Bridges.IndexInVector, ) where {T} func = MOI.Bridges.map_function( - typeof(bridge), + bridge, MOI.VectorOfVariables(bridge.variables), i, ) @@ -212,7 +212,7 @@ end function unbridged_map(bridge::SetMapBridge{T}, vi::MOI.VariableIndex) where {T} F = MOI.ScalarAffineFunction{T} - mapped = MOI.Bridges.inverse_map_function(typeof(bridge), vi) + mapped = MOI.Bridges.inverse_map_function(bridge, vi) return Pair{MOI.VariableIndex,F}[bridge.variable=>mapped] end @@ -222,9 +222,10 @@ function unbridged_map( ) where {T} F = MOI.ScalarAffineFunction{T} func = MOI.VectorOfVariables(vis) - funcs = MOI.Bridges.inverse_map_function(typeof(bridge), func) + funcs = MOI.Bridges.inverse_map_function(bridge, func) scalars = MOI.Utilities.eachscalar(funcs) + # FIXME not correct for SetDotProducts, it won't recover the dot product variables return Pair{MOI.VariableIndex,F}[ - bridge.variables[i] => scalars[i] for i in eachindex(vis) + bridge.variables[i] => scalars[i] for i in eachindex(bridge.variables) ] end diff --git a/src/Bridges/set_map.jl b/src/Bridges/set_map.jl index c665c1dd00..11a51a9906 100644 --- a/src/Bridges/set_map.jl +++ b/src/Bridges/set_map.jl @@ -80,6 +80,10 @@ function map_function(::Type{BT}, func, i::IndexInVector) where {BT} return MOI.Utilities.eachscalar(map_function(BT, func))[i.value] end +function map_function(bridge::AbstractBridge, func, i::IndexInVector) + return map_function(typeof(bridge), func, i) +end + """ inverse_map_function(bridge::MOI.Bridges.AbstractBridge, func) inverse_map_function(::Type{BT}, func) where {BT} diff --git a/src/Test/test_conic.jl b/src/Test/test_conic.jl index b2dca3f3da..48f05000a6 100644 --- a/src/Test/test_conic.jl +++ b/src/Test/test_conic.jl @@ -5770,6 +5770,185 @@ function setup_test( return end +""" +The goal is to find the maximum lower bound `γ` for the polynomial `x^2 - 2x`. +Using samples `-1` and `1`, the polynomial `x^2 - 2x - γ` evaluates at `-γ` +and `2 - γ` respectively. +The dot product with the gram matrix is the evaluation of `[1; x] * [1 x]` hence +`[1; -1] * [1 -1]` and `[1; 1] * [1 1]` respectively. + +The polynomial version is: +max γ +s.t. [-γ, 2 - γ] in SetDotProducts( + PSD(2), + [[1; -1] * [1 -1], [1; 1] * [1 1]], +) +Its dual (moment version) is: +min -y[1] - y[2] +s.t. [-γ, 2 - γ] in LinearCombinationInSet( + PSD(2), + [[1; -1] * [1 -1], [1; 1] * [1 1]], +) +""" +function test_conic_PositiveSemidefinite_RankOne_polynomial( + model::MOI.ModelLike, + config::Config{T}, +) where {T} + set = MOI.SetDotProducts( + MOI.PositiveSemidefiniteConeTriangle(2), + MOI.TriangleVectorization.([ + MOI.PositiveSemidefiniteFactorization(T[1, -1]), + MOI.PositiveSemidefiniteFactorization(T[1, 1]), + ]), + ) + @requires MOI.supports_constraint( + model, + MOI.VectorAffineFunction{T}, + typeof(set), + ) + @requires MOI.supports_incremental_interface(model) + @requires MOI.supports(model, MOI.ObjectiveSense()) + @requires MOI.supports(model, MOI.ObjectiveFunction{MOI.VariableIndex}()) + γ = MOI.add_variable(model) + c = MOI.add_constraint( + model, + MOI.Utilities.operate(vcat, T, T(3) - T(1) * γ, T(-1) - T(1) * γ), + set, + ) + MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE) + MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), γ) + if _supports(config, MOI.optimize!) + @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED + MOI.optimize!(model) + @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status + @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT + if _supports(config, MOI.ConstraintDual) + @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT + end + @test ≈(MOI.get(model, MOI.ObjectiveValue()), T(-1), config) + if _supports(config, MOI.DualObjectiveValue) + @test ≈(MOI.get(model, MOI.DualObjectiveValue()), T(-1), config) + end + @test ≈(MOI.get(model, MOI.VariablePrimal(), γ), T(-1), config) + @test ≈(MOI.get(model, MOI.ConstraintPrimal(), c), T[4, 0], config) + if _supports(config, MOI.ConstraintDual) + @test ≈(MOI.get(model, MOI.ConstraintDual(), c), T[0, 1], config) + end + end + return +end + +function setup_test( + ::typeof(test_conic_PositiveSemidefinite_RankOne_polynomial), + model::MOIU.MockOptimizer, + ::Config{T}, +) where {T<:Real} + A = MOI.TriangleVectorization{ + T, + MOI.PositiveSemidefiniteFactorization{T,Vector{T}}, + } + MOIU.set_mock_optimize!( + model, + (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!( + mock, + T[-1], + ( + MOI.VectorAffineFunction{T}, + MOI.SetDotProducts{ + MOI.PositiveSemidefiniteConeTriangle, + A, + Vector{A}, + }, + ) => [T[0, 1]], + ), + ) + return +end + +""" +The moment version of `test_conic_PositiveSemidefinite_RankOne_polynomial` + +We look for a measure `μ = y1 * δ_{-1} + y2 * δ_{1}` where `δ_{c}` is the Dirac +measure centered at `c`. The objective is +`⟨μ, x^2 - 2x⟩ = y1 * ⟨δ_{-1}, x^2 - 2x⟩ + y2 * ⟨δ_{1}, x^2 - 2x⟩ = 3y1 - y2`. +We want `μ` to be a probability measure so `1 = ⟨μ, 1⟩ = y1 + y2`. +""" +function test_conic_PositiveSemidefinite_RankOne_moment( + model::MOI.ModelLike, + config::Config{T}, +) where {T} + set = MOI.LinearCombinationInSet( + MOI.PositiveSemidefiniteConeTriangle(2), + MOI.TriangleVectorization.([ + MOI.PositiveSemidefiniteFactorization(T[1, -1]), + MOI.PositiveSemidefiniteFactorization(T[1, 1]), + ]), + ) + @requires MOI.supports_add_constrained_variables(model, typeof(set)) + @requires MOI.supports_incremental_interface(model) + @requires MOI.supports(model, MOI.ObjectiveSense()) + @requires MOI.supports( + model, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(), + ) + y, cy = MOI.add_constrained_variables(model, set) + c = MOI.add_constraint(model, T(1) * y[1] + T(1) * y[2], MOI.EqualTo(T(1))) + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.set( + model, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(), + T(3) * y[1] - T(1) * y[2], + ) + if _supports(config, MOI.optimize!) + @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED + MOI.optimize!(model) + @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status + @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT + if _supports(config, MOI.ConstraintDual) + @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT + end + @test ≈(MOI.get(model, MOI.ObjectiveValue()), T(-1), config) + if _supports(config, MOI.DualObjectiveValue) + @test ≈(MOI.get(model, MOI.DualObjectiveValue()), T(-1), config) + end + @test ≈(MOI.get(model, MOI.VariablePrimal(), y), T[0, 1], config) + @test ≈(MOI.get(model, MOI.ConstraintPrimal(), c), T(1), config) + if _supports(config, MOI.ConstraintDual) + @test ≈(MOI.get(model, MOI.ConstraintDual(), cy), T[4, 0], config) + @test ≈(MOI.get(model, MOI.ConstraintDual(), c), T(-1), config) + end + end + return +end + +function setup_test( + ::typeof(test_conic_PositiveSemidefinite_RankOne_moment), + model::MOIU.MockOptimizer, + ::Config{T}, +) where {T<:Real} + A = MOI.TriangleVectorization{ + T, + MOI.PositiveSemidefiniteFactorization{T,Vector{T}}, + } + MOIU.set_mock_optimize!( + model, + (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!( + mock, + T[0, 1], + (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}) => [T(-1)], + ( + MOI.VectorOfVariables, + MOI.LinearCombinationInSet{ + MOI.PositiveSemidefiniteConeTriangle, + A, + Vector{A}, + }, + ) => [T[4, 0]], + ), + ) + return +end + """ _test_det_cone_helper_ellipsoid( model::MOI.ModelLike, diff --git a/src/Utilities/functions.jl b/src/Utilities/functions.jl index e3f0d5a172..a9fb462e48 100644 --- a/src/Utilities/functions.jl +++ b/src/Utilities/functions.jl @@ -638,18 +638,24 @@ end A type that allows iterating over the scalar-functions that comprise an `AbstractVectorFunction`. """ -struct ScalarFunctionIterator{F<:MOI.AbstractVectorFunction,C} +struct ScalarFunctionIterator{F<:MOI.AbstractVectorFunction,C,S} <: + AbstractVector{S} f::F # Cache that can be used to store a precomputed datastructure that allows # an efficient implementation of `getindex`. cache::C + function ScalarFunctionIterator(f::MOI.AbstractVectorFunction, cache) + return new{typeof(f),typeof(cache),scalar_type(typeof(f))}(f, cache) + end end function ScalarFunctionIterator(func::MOI.AbstractVectorFunction) return ScalarFunctionIterator(func, scalar_iterator_cache(func)) end -scalar_iterator_cache(func::MOI.AbstractVectorFunction) = nothing +Base.size(s::ScalarFunctionIterator) = (MOI.output_dimension(s.f),) + +scalar_iterator_cache(::MOI.AbstractVectorFunction) = nothing function output_index_iterator(terms::AbstractVector, output_dimension) start = zeros(Int, output_dimension) diff --git a/src/sets.jl b/src/sets.jl index 534106e9de..b85e0bafc1 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -1791,6 +1791,157 @@ function Base.getproperty( end end +""" + SetDotProducts(set::MOI.AbstractSet, vectors::AbstractVector) + +Given a set `set` of dimension `d` and `m` vectors `a_1`, ..., `a_m` given in `vectors`, this is the set: +``\\{ ((\\langle a_1, x \\rangle, ..., \\langle a_m, x \\rangle) \\in \\mathbb{R}^{m} : x \\in \\text{set} \\}.`` +""" +struct SetDotProducts{S,A,V<:AbstractVector{A}} <: AbstractVectorSet + set::S + vectors::V +end + +function Base.:(==)(s1::SetDotProducts, s2::SetDotProducts) + return s1.set == s2.set && s1.vectors == s2.vectors +end + +function Base.copy(s::SetDotProducts) + return SetDotProducts(copy(s.set), copy(s.vectors)) +end + +dimension(s::SetDotProducts) = length(s.vectors) + +function dual_set(s::SetDotProducts) + return LinearCombinationInSet(s.set, s.vectors) +end + +function dual_set_type(::Type{SetDotProducts{S,A,V}}) where {S,A,V} + return LinearCombinationInSet{S,A,V} +end + +""" + LinearCombinationInSet(set::MOI.AbstractSet, matrices::AbstractVector) + +Given a set `set` of dimension `d` and `m` vectors `a_1`, ..., `a_m` given in `vectors`, this is the set: +``\\{ (y \\in \\mathbb{R}^{m} : \\sum_{i=1}^m y_i a_i \\in \\text{set} \\}.`` +""" +struct LinearCombinationInSet{S,A,V<:AbstractVector{A}} <: AbstractVectorSet + set::S + vectors::V +end + +function Base.:(==)(s1::LinearCombinationInSet, s2::LinearCombinationInSet) + return s1.set == s2.set && s1.vectors == s2.vectors +end + +function Base.copy(s::LinearCombinationInSet) + return LinearCombinationInSet(copy(s.set), copy(s.vectors)) +end + +dimension(s::LinearCombinationInSet) = length(s.vectors) + +function dual_set(s::LinearCombinationInSet) + return SetDotProducts(s.side_dimension, s.vectors) +end + +function dual_set_type(::Type{LinearCombinationInSet{S,A,V}}) where {S,A,V} + return SetDotProducts{S,A,V} +end + +abstract type AbstractFactorization{T,F} <: AbstractMatrix{T} end + +function Base.size(m::AbstractFactorization) + n = size(m.factor, 1) + return (n, n) +end + +""" + struct Factorization{ + T, + F<:Union{AbstractVector{T},AbstractMatrix{T}}, + D<:Union{T,AbstractVector{T}}, + } <: AbstractMatrix{T} + factor::F + scaling::D + end + +Matrix corresponding to `factor * Diagonal(diagonal) * factor'`. +If `factor` is a vector and `diagonal` is a scalar, this corresponds to +the matrix `diagonal * factor * factor'`. +If `factor` is a matrix and `diagonal` is a vector, this corresponds to +the matrix `factor * Diagonal(scaling) * factor'`. +""" +struct Factorization{ + T, + F<:Union{AbstractVector{T},AbstractMatrix{T}}, + D<:Union{T,AbstractVector{T}}, +} <: AbstractFactorization{T,F} + factor::F + scaling::D + function Factorization( + factor::AbstractMatrix{T}, + scaling::AbstractVector{T}, + ) where {T} + if length(scaling) != size(factor, 2) + error( + "Length `$(length(scaling))` of diagonal does not match number of columns `$(size(factor, 2))` of factor", + ) + end + return new{T,typeof(factor),typeof(scaling)}(factor, scaling) + end + function Factorization(factor::AbstractVector{T}, scaling::T) where {T} + return new{T,typeof(factor),typeof(scaling)}(factor, scaling) + end +end + +function Base.getindex(m::Factorization, i::Int, j::Int) + return sum( + m.factor[i, k] * m.scaling[k] * m.factor[j, k]' for + k in eachindex(m.scaling) + ) +end + +""" + struct Factorization{ + T, + F<:Union{AbstractVector{T},AbstractMatrix{T}}, + D<:Union{T,AbstractVector{T}}, + } <: AbstractMatrix{T} + factor::F + scaling::D + end + +Matrix corresponding to `factor * Diagonal(diagonal) * factor'`. +If `factor` is a vector and `diagonal` is a scalar, this corresponds to +the matrix `diagonal * factor * factor'`. +If `factor` is a matrix and `diagonal` is a vector, this corresponds to +the matrix `factor * Diagonal(scaling) * factor'`. +""" +struct PositiveSemidefiniteFactorization{ + T, + F<:Union{AbstractVector{T},AbstractMatrix{T}}, +} <: AbstractFactorization{T,F} + factor::F +end + +function Base.getindex(m::PositiveSemidefiniteFactorization, i::Int, j::Int) + return sum(m.factor[i, k] * m.factor[j, k]' for k in axes(m.factor, 2)) +end + +struct TriangleVectorization{T,M<:AbstractMatrix{T}} <: AbstractVector{T} + matrix::M +end + +function Base.size(v::TriangleVectorization) + n = size(v.matrix, 1) + return (Utilities.trimap(n, n),) +end + +function Base.getindex(v::TriangleVectorization, k::Int) + return getindex(v.matrix, Utilities.inverse_trimap(k)...) +end + """ SOS1{T<:Real}(weights::Vector{T}) diff --git a/test/Bridges/Variable/set_dot.jl b/test/Bridges/Variable/set_dot.jl new file mode 100644 index 0000000000..b9a52377a9 --- /dev/null +++ b/test/Bridges/Variable/set_dot.jl @@ -0,0 +1,72 @@ +# Copyright (c) 2017: Miles Lubin and contributors +# Copyright (c) 2017: Google Inc. +# +# 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. + +module TestVariableDotProducts + +using Test + +import MathOptInterface as MOI + +function runtests() + for name in names(@__MODULE__; all = true) + if startswith("$(name)", "test_") + @testset "$(name)" begin + getfield(@__MODULE__, name)() + end + end + end + return +end + +include("../utilities.jl") + +function test_psd() + MOI.Bridges.runtests( + MOI.Bridges.Variable.DotProductsBridge, + model -> begin + x, _ = MOI.add_constrained_variables( + model, + MOI.SetDotProducts( + MOI.PositiveSemidefiniteConeTriangle(2), + MOI.TriangleVectorization.([ + [ + 1 2.0 + 2 3 + ], + [ + 4 5.0 + 5 6 + ], + ]), + ), + ) + MOI.add_constraint(model, 1.0x[1], MOI.EqualTo(0.0)) + MOI.add_constraint(model, 1.0x[2], MOI.LessThan(0.0)) + end, + model -> begin + Q, _ = MOI.add_constrained_variables( + model, + MOI.PositiveSemidefiniteConeTriangle(2), + ) + MOI.add_constraint( + model, + 1.0 * Q[1] + 4.0 * Q[2] + 3.0 * Q[3], + MOI.EqualTo(0.0), + ) + MOI.add_constraint( + model, + 4.0 * Q[1] + 10.0 * Q[2] + 6.0 * Q[3], + MOI.LessThan(0.0), + ) + end; + cannot_unbridge = true, + ) + return +end + +end # module + +TestVariableDotProducts.runtests() diff --git a/test/sets.jl b/test/sets.jl index 28b0c5000a..3e25051c12 100644 --- a/test/sets.jl +++ b/test/sets.jl @@ -8,6 +8,7 @@ module TestSets using Test import MathOptInterface as MOI +import LinearAlgebra include("dummy.jl") @@ -353,6 +354,41 @@ function test_sets_reified() return end +function _test_factorization(A, B) + @test size(A) == size(B) + @test A ≈ B + d = LinearAlgebra.checksquare(A) + n = div(d * (d + 1), 2) + vA = MOI.TriangleVectorization(A) + @test length(vA) == n + @test eachindex(vA) == Base.OneTo(n) + vB = MOI.TriangleVectorization(B) + @test length(vB) == n + @test eachindex(vA) == Base.OneTo(n) + k = 0 + for j in 1:d + for i in 1:j + k += 1 + @test vA[k] == vB[k] + @test vA[k] == A[i, j] + end + end + return +end + +function test_factorizations() + f = [1, 2] + _test_factorization(f * f', MOI.PositiveSemidefiniteFactorization(f)) + _test_factorization(2 * f * f', MOI.Factorization(f, 2)) + F = [1 2; 3 4; 5 6] + d = [7, 8] + _test_factorization(F * F', MOI.PositiveSemidefiniteFactorization(F)) + return _test_factorization( + F * LinearAlgebra.Diagonal(d) * F', + MOI.Factorization(F, d), + ) +end + function runtests() for name in names(@__MODULE__; all = true) if startswith("$name", "test_")