From b36bb74e7144e2e9d2a9988e87f300210be5491d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 5 Jul 2023 17:46:17 +0200 Subject: [PATCH 1/3] Add Scaled set --- docs/src/changelog.md | 2 +- docs/src/manual/standard_form.md | 2 +- docs/src/reference/standard_form.md | 2 +- .../src/submodules/Bridges/list_of_bridges.md | 4 +- src/Bridges/Constraint/Constraint.jl | 9 +- .../Constraint/bridges/set_dot_scaling.jl | 212 +++++++++++ .../bridges/symmetric_matrix_scaling.jl | 200 ---------- src/FileFormats/MOF/MOF.jl | 2 +- src/FileFormats/MOF/read.jl | 5 +- src/FileFormats/MOF/write.jl | 13 +- src/Test/test_basic_constraint.jl | 4 +- src/Test/test_conic.jl | 6 +- src/Utilities/Utilities.jl | 1 + src/Utilities/matrix_of_constraints.jl | 11 +- src/Utilities/model.jl | 2 +- src/Utilities/results.jl | 140 ------- src/Utilities/set_dot.jl | 352 ++++++++++++++++++ src/Utilities/sets.jl | 39 -- src/sets.jl | 258 ++++++++----- ...c_matrix_scaling.jl => set_dot_scaling.jl} | 16 +- test/Bridges/lazy_bridge_optimizer.jl | 2 +- test/Utilities/matrix_of_constraints.jl | 6 + test/Utilities/set_dot.jl | 53 +++ test/Utilities/sets.jl | 11 +- 24 files changed, 847 insertions(+), 505 deletions(-) create mode 100644 src/Bridges/Constraint/bridges/set_dot_scaling.jl delete mode 100644 src/Bridges/Constraint/bridges/symmetric_matrix_scaling.jl create mode 100644 src/Utilities/set_dot.jl rename test/Bridges/Constraint/{symmetric_matrix_scaling.jl => set_dot_scaling.jl} (84%) create mode 100644 test/Utilities/set_dot.jl diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 249fd42543..3569cdb7e7 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -109,7 +109,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - [`hessian_objective_structure`](@ref) - Added new sets - [`NormCone`](@ref) (#2119) - - [`ScaledPositiveSemidefiniteConeTriangle`](@ref) (#2154) + - `ScaledPositiveSemidefiniteConeTriangle` (#2154) ### Fixed diff --git a/docs/src/manual/standard_form.md b/docs/src/manual/standard_form.md index 431c1fc03b..9458cf3355 100644 --- a/docs/src/manual/standard_form.md +++ b/docs/src/manual/standard_form.md @@ -93,13 +93,13 @@ The matrix-valued set types implemented in MathOptInterface.jl are: | [`RootDetConeSquare(d)`](@ref MathOptInterface.RootDetConeSquare) | ``\{ (t,X) \in \mathbb{R}^{1+d^2} : t \le \det(X)^{1/d}, X \mbox{ is a PSD matrix} \}`` | | [`PositiveSemidefiniteConeTriangle(d)`](@ref MathOptInterface.PositiveSemidefiniteConeTriangle) | ``\{ X \in \mathbb{R}^{d(d+1)/2} : X \mbox{ is the upper triangle of a PSD matrix} \}`` | | [`PositiveSemidefiniteConeSquare(d)`](@ref MathOptInterface.PositiveSemidefiniteConeSquare) | ``\{ X \in \mathbb{R}^{d^2} : X \mbox{ is a PSD matrix} \}`` | -| [`ScaledPositiveSemidefiniteConeTriangle(d)`](@ref MathOptInterface.ScaledPositiveSemidefiniteConeTriangle) | ``\{ X \in \mathbb{R}^{d(d+1)/2} : X \mbox{ is a PSD matrix} \}`` | | [`LogDetConeTriangle(d)`](@ref MathOptInterface.LogDetConeTriangle) | ``\{ (t,u,X) \in \mathbb{R}^{2+d(1+d)/2} : t \le u\log(\det(X/u)), X \mbox{ is the upper triangle of a PSD matrix}, u > 0 \}`` | | [`LogDetConeSquare(d)`](@ref MathOptInterface.LogDetConeSquare) | ``\{ (t,u,X) \in \mathbb{R}^{2+d^2} : t \le u \log(\det(X/u)), X \mbox{ is a PSD matrix}, u > 0 \}`` | | [`NormSpectralCone(r, c)`](@ref MathOptInterface.NormSpectralCone) | ``\{ (t, X) \in \mathbb{R}^{1 + r \times c} : t \ge \sigma_1(X), X \mbox{ is a } r\times c\mbox{ matrix} \}`` | [`NormNuclearCone(r, c)`](@ref MathOptInterface.NormNuclearCone) | ``\{ (t, X) \in \mathbb{R}^{1 + r \times c} : t \ge \sum_i \sigma_i(X), X \mbox{ is a } r\times c\mbox{ matrix} \}`` | | [`HermitianPositiveSemidefiniteConeTriangle(d)`](@ref MathOptInterface.HermitianPositiveSemidefiniteConeTriangle) | The cone of Hermitian positive semidefinite matrices, with `side_dimension` rows and columns. | +| [`Scaled(S)`](@ref MathOptInterface.Scaled) | The set `S` scaled so that [`Utilities.set_dot`](@ref MathOptInterface.Utilities.set_dot) corresponds to `LinearAlgebra.dot` | Some of these cones can take two forms: `XXXConeTriangle` and `XXXConeSquare`. diff --git a/docs/src/reference/standard_form.md b/docs/src/reference/standard_form.md index 9b4e43ab50..98ff739a33 100644 --- a/docs/src/reference/standard_form.md +++ b/docs/src/reference/standard_form.md @@ -100,6 +100,7 @@ SOS2 Indicator Complements HyperRectangle +Scaled ``` ## Constraint programming sets @@ -148,7 +149,6 @@ List of recognized matrix sets. PositiveSemidefiniteConeTriangle PositiveSemidefiniteConeSquare HermitianPositiveSemidefiniteConeTriangle -ScaledPositiveSemidefiniteConeTriangle LogDetConeTriangle LogDetConeSquare RootDetConeTriangle diff --git a/docs/src/submodules/Bridges/list_of_bridges.md b/docs/src/submodules/Bridges/list_of_bridges.md index a83c5a6858..3a7c8e2c4b 100644 --- a/docs/src/submodules/Bridges/list_of_bridges.md +++ b/docs/src/submodules/Bridges/list_of_bridges.md @@ -57,8 +57,8 @@ Bridges.Constraint.NormSpectralBridge Bridges.Constraint.NormNuclearBridge Bridges.Constraint.SquareBridge Bridges.Constraint.HermitianToSymmetricPSDBridge -Bridges.Constraint.SymmetricMatrixScalingBridge -Bridges.Constraint.SymmetricMatrixInverseScalingBridge +Bridges.Constraint.SetDotScalingBridge +Bridges.Constraint.SetDotInverseScalingBridge Bridges.Constraint.RootDetBridge Bridges.Constraint.LogDetBridge Bridges.Constraint.IndicatorActiveOnFalseBridge diff --git a/src/Bridges/Constraint/Constraint.jl b/src/Bridges/Constraint/Constraint.jl index 3e14f12434..576bec6cd3 100644 --- a/src/Bridges/Constraint/Constraint.jl +++ b/src/Bridges/Constraint/Constraint.jl @@ -58,7 +58,7 @@ include("bridges/split_complex_zeros.jl") include("bridges/split_hyperrectangle.jl") include("bridges/hermitian.jl") include("bridges/square.jl") -include("bridges/symmetric_matrix_scaling.jl") +include("bridges/set_dot_scaling.jl") include("bridges/table.jl") include("bridges/vectorize.jl") include("bridges/zero_one.jl") @@ -109,11 +109,8 @@ function add_all_bridges(bridged_model, ::Type{T}) where {T} MOI.Bridges.add_bridge(bridged_model, NormNuclearBridge{T}) MOI.Bridges.add_bridge(bridged_model, HermitianToSymmetricPSDBridge{T}) MOI.Bridges.add_bridge(bridged_model, SquareBridge{T}) - MOI.Bridges.add_bridge(bridged_model, SymmetricMatrixScalingBridge{T}) - MOI.Bridges.add_bridge( - bridged_model, - SymmetricMatrixInverseScalingBridge{T}, - ) + MOI.Bridges.add_bridge(bridged_model, SetDotScalingBridge{T}) + MOI.Bridges.add_bridge(bridged_model, SetDotInverseScalingBridge{T}) MOI.Bridges.add_bridge(bridged_model, LogDetBridge{T}) MOI.Bridges.add_bridge(bridged_model, RootDetBridge{T}) MOI.Bridges.add_bridge(bridged_model, RSOCtoSOCBridge{T}) diff --git a/src/Bridges/Constraint/bridges/set_dot_scaling.jl b/src/Bridges/Constraint/bridges/set_dot_scaling.jl new file mode 100644 index 0000000000..de9ec6680a --- /dev/null +++ b/src/Bridges/Constraint/bridges/set_dot_scaling.jl @@ -0,0 +1,212 @@ +# 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. + +""" + SetDotScalingBridge{T,S,F,G} <: Bridges.Constraint.AbstractBridge + +`SetDotScalingBridge` implements the reformulation from constraints +in `S` to constraints in [`MOI.Scaled{S}`](@ref MOI.Scaled). + +## Source node + +`SetDotScalingBridge` supports: + + * `G` in `S` + +## Target node + +`SetDotScalingBridge` creates: + + * `F` in [`MOI.Scaled{S}`](@ref MOI.Scaled) +""" +struct SetDotScalingBridge{T,S,F,G} <: SetMapBridge{T,MOI.Scaled{S},S,F,G} + constraint::MOI.ConstraintIndex{F,MOI.Scaled{S}} +end + +const SetDotScaling{T,OT<:MOI.ModelLike} = + SingleBridgeOptimizer{SetDotScalingBridge{T},OT} + +function concrete_bridge_type( + ::Type{<:SetDotScalingBridge{T}}, + G::Type{<:MOI.AbstractVectorFunction}, + ::Type{S}, +) where {T,S<:MOI.AbstractVectorSet} + F = MOI.Utilities.promote_operation( + *, + T, + LinearAlgebra.Diagonal{T,Vector{T}}, + G, + ) + return SetDotScalingBridge{T,S,F,G} +end + +function MOI.Bridges.map_set( + ::Type{<:SetDotScalingBridge{T,S}}, + set::S, +) where {T,S} + return MOI.Scaled(set) +end + +function MOI.Bridges.inverse_map_set( + ::Type{<:SetDotScalingBridge{T,S}}, + set::MOI.Scaled{S}, +) where {T,S} + return set.set +end + +_length(f::MOI.AbstractVectorFunction) = MOI.output_dimension(f) +_length(f::AbstractVector) = length(f) + +function _scale(::Type{T}, ::Type{S}, func) where {T,S} + set = MOI.Utilities.set_with_dimension(S, _length(func)) + scale = MOI.Utilities.SetDotScalingVector{T}(set) + return MOI.Utilities.operate(*, T, LinearAlgebra.Diagonal(scale), func) +end + +function _inverse_scale(::Type{T}, ::Type{S}, func) where {T,S} + set = MOI.Utilities.set_with_dimension(S, _length(func)) + scale = MOI.Utilities.SetDotScalingVector{T}(set) + inv_scale = MOI.Utilities.lazy_map(T, inv, scale) + return MOI.Utilities.operate(*, T, LinearAlgebra.Diagonal(inv_scale), func) +end + +function MOI.Bridges.map_function( + ::Type{<:SetDotScalingBridge{T,S}}, + func, +) where {T,S} + return _scale(T, S, func) +end + +function MOI.Bridges.inverse_map_function( + ::Type{<:SetDotScalingBridge{T,S}}, + func, +) where {T,S} + return _inverse_scale(T, S, func) +end + +# Since the map is a diagonal matrix `D`, it is symmetric so one would initially +# expect `adjoint_map_function` to be the same as `map_function`. However, the +# scalar product for the scaled PSD cone is `_2 = x'y` but the scalar +# product for the PSD cone additionally scales the offdiagonal entries by `2` +# hence by `D^2` so `_1 = x'D^2y`. +# So `_2 = _1` hence the adjoint of `D` is its inverse! +function MOI.Bridges.adjoint_map_function( + ::Type{<:SetDotScalingBridge{T,S}}, + func, +) where {T,S} + return _inverse_scale(T, S, func) +end + +function MOI.Bridges.inverse_adjoint_map_function( + ::Type{<:SetDotScalingBridge{T,S}}, + func, +) where {T,S} + return _scale(T, S, func) +end + +""" + SetDotInverseScalingBridge{T,S,F,G} <: Bridges.Constraint.AbstractBridge + +`SetDotInverseScalingBridge` implements the reformulation from constraints +in the `MOI.Scaled{S}` to constraints in the `S`. + +## Source node + +`SetDotInverseScalingBridge` supports: + + * `G` in [`MOI.Scaled{S}`](@ref MOI.Scaled) + +## Target node + +`SetDotInverseScalingBridge` creates: + + * `F` in `S` +""" +struct SetDotInverseScalingBridge{T,S,F,G} <: + SetMapBridge{T,S,MOI.Scaled{S},F,G} + constraint::MOI.ConstraintIndex{F,S} +end + +const SetDotInverseScaling{T,OT<:MOI.ModelLike} = + SingleBridgeOptimizer{SetDotInverseScalingBridge{T},OT} + +function concrete_bridge_type( + ::Type{<:SetDotInverseScalingBridge{T}}, + G::Type{<:MOI.AbstractVectorFunction}, + ::Type{MOI.Scaled{S}}, +) where {T,S<:MOI.AbstractVectorSet} + F = MOI.Utilities.promote_operation( + *, + T, + LinearAlgebra.Diagonal{T,Vector{T}}, + G, + ) + return SetDotInverseScalingBridge{T,S,F,G} +end + +function MOI.Bridges.map_set( + ::Type{<:SetDotInverseScalingBridge{T,S}}, + set::MOI.Scaled{S}, +) where {T,S} + return set.set +end + +function MOI.Bridges.inverse_map_set( + ::Type{<:SetDotInverseScalingBridge{T,S}}, + set::S, +) where {T,S} + return MOI.Scaled(set) +end + +function MOI.Bridges.map_function( + ::Type{<:SetDotInverseScalingBridge{T,S}}, + func, +) where {T,S} + return _inverse_scale(T, S, func) +end + +function MOI.Bridges.inverse_map_function( + ::Type{<:SetDotInverseScalingBridge{T,S}}, + func, +) where {T,S} + return _scale(T, S, func) +end + +function MOI.Bridges.adjoint_map_function( + ::Type{<:SetDotInverseScalingBridge{T,S}}, + func, +) where {T,S} + return _inverse_scale(T, S, func) +end + +function MOI.Bridges.inverse_adjoint_map_function( + ::Type{<:SetDotInverseScalingBridge{T,S}}, + func, +) where {T,S} + return _scale(T, S, func) +end + +# Since the set type is not defined, the default `MOI.supports_constraint` +# for `SetMapBridge` does not work +function MOI.supports_constraint( + ::Type{<:SetDotScalingBridge}, + ::Type{<:MOI.AbstractVectorFunction}, + S::Type{<:MOI.AbstractVectorSet}, +) + return MOI.is_set_dot_scaled(S) +end + +function MOI.supports_constraint( + ::Type{<:SetDotInverseScalingBridge}, + ::Type{<:MOI.AbstractVectorFunction}, + ::Type{MOI.Scaled{S}}, +) where {S<:MOI.AbstractVectorSet} + return true +end + +# TODO remove in MOI v2 +const SymmetricMatrixScalingBridge = SetDotScalingBridge +const SymmetricMatrixInverseScalingBridge = SetDotInverseScalingBridge diff --git a/src/Bridges/Constraint/bridges/symmetric_matrix_scaling.jl b/src/Bridges/Constraint/bridges/symmetric_matrix_scaling.jl deleted file mode 100644 index f46cbc37f6..0000000000 --- a/src/Bridges/Constraint/bridges/symmetric_matrix_scaling.jl +++ /dev/null @@ -1,200 +0,0 @@ -# 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. - -""" - SymmetricMatrixScalingBridge{T,F,G} <: Bridges.Constraint.AbstractBridge - -`SymmetricMatrixScalingBridge` implements the reformulation from constraints -in the `MOI.PositiveSemidefiniteConeTriangle` to constraints -in the `MOI.ScaledPositiveSemidefiniteConeTriangle`. - -## Source node - -`SymmetricMatrixScalingBridge` supports: - - * `G` in [`MOI.PositiveSemidefiniteConeTriangle`](@ref) - -## Target node - -`SymmetricMatrixScalingBridge` creates: - - * `F` in [`MOI.ScaledPositiveSemidefiniteConeTriangle`](@ref) -""" -struct SymmetricMatrixScalingBridge{T,F,G} <: SetMapBridge{ - T, - MOI.ScaledPositiveSemidefiniteConeTriangle, - MOI.PositiveSemidefiniteConeTriangle, - F, - G, -} - constraint::MOI.ConstraintIndex{ - F, - MOI.ScaledPositiveSemidefiniteConeTriangle, - } -end - -const SymmetricMatrixScaling{T,OT<:MOI.ModelLike} = - SingleBridgeOptimizer{SymmetricMatrixScalingBridge{T},OT} - -function concrete_bridge_type( - ::Type{<:SymmetricMatrixScalingBridge{T}}, - G::Type{<:MOI.AbstractVectorFunction}, - ::Type{MOI.PositiveSemidefiniteConeTriangle}, -) where {T} - F = MOI.Utilities.promote_operation( - *, - T, - LinearAlgebra.Diagonal{T,Vector{T}}, - G, - ) - return SymmetricMatrixScalingBridge{T,F,G} -end - -function MOI.Bridges.map_set( - ::Type{<:SymmetricMatrixScalingBridge}, - set::MOI.PositiveSemidefiniteConeTriangle, -) - return MOI.ScaledPositiveSemidefiniteConeTriangle(set.side_dimension) -end - -function MOI.Bridges.inverse_map_set( - ::Type{<:SymmetricMatrixScalingBridge}, - set::MOI.ScaledPositiveSemidefiniteConeTriangle, -) - return MOI.PositiveSemidefiniteConeTriangle(set.side_dimension) -end - -_length(f::MOI.AbstractVectorFunction) = MOI.output_dimension(f) -_length(f::AbstractVector) = length(f) - -function MOI.Bridges.map_function( - ::Type{<:SymmetricMatrixScalingBridge{T}}, - func, -) where {T} - scale = MOI.Utilities.symmetric_matrix_scaling_vector(T, _length(func)) - return MOI.Utilities.operate(*, T, LinearAlgebra.Diagonal(scale), func) -end - -function MOI.Bridges.inverse_map_function( - ::Type{<:SymmetricMatrixScalingBridge{T}}, - func, -) where {T} - N = _length(func) - scale = MOI.Utilities.symmetric_matrix_inverse_scaling_vector(T, N) - return MOI.Utilities.operate(*, T, LinearAlgebra.Diagonal(scale), func) -end - -# Since the map is a diagonal matrix `D`, it is symmetric so one would initially -# expect `adjoint_map_function` to be the same as `map_function`. However, the -# scalar product for the scaled PSD cone is `_2 = x'y` but the scalar -# product for the PSD cone additionally scales the offdiagonal entries by `2` -# hence by `D^2` so `_1 = x'D^2y`. -# So `_2 = _1` hence the adjoint of `D` is its inverse! -function MOI.Bridges.adjoint_map_function( - BT::Type{<:SymmetricMatrixScalingBridge}, - func, -) - return MOI.Bridges.inverse_map_function(BT, func) -end - -function MOI.Bridges.inverse_adjoint_map_function( - BT::Type{<:SymmetricMatrixScalingBridge}, - func, -) - return MOI.Bridges.map_function(BT, func) -end - -""" - SymmetricMatrixInverseScalingBridge{T,F,G} <: Bridges.Constraint.AbstractBridge - -`SymmetricMatrixInverseScalingBridge` implements the reformulation from constraints -in the `MOI.ScaledPositiveSemidefiniteConeTriangle` to constraints -in the `MOI.PositiveSemidefiniteConeTriangle`. - -## Source node - -`SymmetricMatrixInverseScalingBridge` supports: - - * `G` in [`MOI.ScaledPositiveSemidefiniteConeTriangle`](@ref) - -## Target node - -`SymmetricMatrixInverseScalingBridge` creates: - - * `F` in [`MOI.PositiveSemidefiniteConeTriangle`](@ref) -""" -struct SymmetricMatrixInverseScalingBridge{T,F,G} <: SetMapBridge{ - T, - MOI.PositiveSemidefiniteConeTriangle, - MOI.ScaledPositiveSemidefiniteConeTriangle, - F, - G, -} - constraint::MOI.ConstraintIndex{F,MOI.PositiveSemidefiniteConeTriangle} -end - -const SymmetricMatrixInverseScaling{T,OT<:MOI.ModelLike} = - SingleBridgeOptimizer{SymmetricMatrixInverseScalingBridge{T},OT} - -function concrete_bridge_type( - ::Type{<:SymmetricMatrixInverseScalingBridge{T}}, - G::Type{<:MOI.AbstractVectorFunction}, - ::Type{MOI.ScaledPositiveSemidefiniteConeTriangle}, -) where {T} - F = MOI.Utilities.promote_operation( - *, - T, - LinearAlgebra.Diagonal{T,Vector{T}}, - G, - ) - return SymmetricMatrixInverseScalingBridge{T,F,G} -end - -function MOI.Bridges.map_set( - ::Type{<:SymmetricMatrixInverseScalingBridge}, - set::MOI.ScaledPositiveSemidefiniteConeTriangle, -) - return MOI.PositiveSemidefiniteConeTriangle(set.side_dimension) -end - -function MOI.Bridges.inverse_map_set( - ::Type{<:SymmetricMatrixInverseScalingBridge}, - set::MOI.PositiveSemidefiniteConeTriangle, -) - return MOI.ScaledPositiveSemidefiniteConeTriangle(set.side_dimension) -end - -function MOI.Bridges.map_function( - ::Type{<:SymmetricMatrixInverseScalingBridge{T}}, - func, -) where {T} - N = _length(func) - scale = MOI.Utilities.symmetric_matrix_inverse_scaling_vector(T, N) - return MOI.Utilities.operate(*, T, LinearAlgebra.Diagonal(scale), func) -end - -function MOI.Bridges.inverse_map_function( - ::Type{<:SymmetricMatrixInverseScalingBridge{T}}, - func, -) where {T} - N = _length(func) - scale = MOI.Utilities.symmetric_matrix_scaling_vector(T, N) - return MOI.Utilities.operate(*, T, LinearAlgebra.Diagonal(scale), func) -end - -function MOI.Bridges.adjoint_map_function( - BT::Type{<:SymmetricMatrixInverseScalingBridge}, - func, -) - return MOI.Bridges.inverse_map_function(BT, func) -end - -function MOI.Bridges.inverse_adjoint_map_function( - BT::Type{<:SymmetricMatrixInverseScalingBridge}, - func, -) - return MOI.Bridges.map_function(BT, func) -end diff --git a/src/FileFormats/MOF/MOF.jl b/src/FileFormats/MOF/MOF.jl index ac3db58f74..5d050721fd 100644 --- a/src/FileFormats/MOF/MOF.jl +++ b/src/FileFormats/MOF/MOF.jl @@ -64,7 +64,7 @@ MOI.Utilities.@model( MOI.LogDetConeTriangle, MOI.LogDetConeSquare, MOI.PositiveSemidefiniteConeTriangle, - MOI.ScaledPositiveSemidefiniteConeTriangle, + MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle}, MOI.PositiveSemidefiniteConeSquare, MOI.HermitianPositiveSemidefiniteConeTriangle, MOI.AllDifferent, diff --git a/src/FileFormats/MOF/read.jl b/src/FileFormats/MOF/read.jl index be391cf4e5..c5846a5882 100644 --- a/src/FileFormats/MOF/read.jl +++ b/src/FileFormats/MOF/read.jl @@ -382,7 +382,7 @@ end LogDetConeTriangle, LogDetConeSquare, PositiveSemidefiniteConeTriangle, - ScaledPositiveSemidefiniteConeTriangle, + ScaledPositiveSemidefiniteConeTriangle, # Required for v1.4 PositiveSemidefiniteConeSquare, HermitianPositiveSemidefiniteConeTriangle, ExponentialCone, @@ -539,7 +539,8 @@ function set_to_moi( ::Val{:ScaledPositiveSemidefiniteConeTriangle}, object::Object, ) - return MOI.ScaledPositiveSemidefiniteConeTriangle(object["side_dimension"]) + d = object["side_dimension"] + return MOI.Scaled(MOI.PositiveSemidefiniteConeTriangle(d)) end function set_to_moi(::Val{:PositiveSemidefiniteConeSquare}, object::Object) diff --git a/src/FileFormats/MOF/write.jl b/src/FileFormats/MOF/write.jl index 7eca362126..c2c3344285 100644 --- a/src/FileFormats/MOF/write.jl +++ b/src/FileFormats/MOF/write.jl @@ -300,9 +300,6 @@ head_name(::Type{MOI.LogDetConeSquare}) = "LogDetConeSquare" function head_name(::Type{MOI.PositiveSemidefiniteConeTriangle}) return "PositiveSemidefiniteConeTriangle" end -function head_name(::Type{MOI.ScaledPositiveSemidefiniteConeTriangle}) - return "ScaledPositiveSemidefiniteConeTriangle" -end function head_name(::Type{MOI.PositiveSemidefiniteConeSquare}) return "PositiveSemidefiniteConeSquare" end @@ -359,3 +356,13 @@ function moi_to_object( "table" => [set.table[i, :] for i in 1:size(set.table, 1)], ) end + +function moi_to_object( + set::MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle}, + name_map::Dict{MOI.VariableIndex,String}, +) + return OrderedObject( + "type" => "ScaledPositiveSemidefiniteConeTriangle", + "side_dimension" => set.set.side_dimension, + ) +end diff --git a/src/Test/test_basic_constraint.jl b/src/Test/test_basic_constraint.jl index 6f7dae5a41..32fd3d698d 100644 --- a/src/Test/test_basic_constraint.jl +++ b/src/Test/test_basic_constraint.jl @@ -131,8 +131,8 @@ end function _set(::Type{MOI.HermitianPositiveSemidefiniteConeTriangle}) return MOI.HermitianPositiveSemidefiniteConeTriangle(3) end -function _set(::Type{MOI.ScaledPositiveSemidefiniteConeTriangle}) - return MOI.ScaledPositiveSemidefiniteConeTriangle(3) +function _set(::Type{MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle}}) + return MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle}(3) end _set(::Type{MOI.LogDetConeTriangle}) = MOI.LogDetConeTriangle(3) _set(::Type{MOI.LogDetConeSquare}) = MOI.LogDetConeSquare(3) diff --git a/src/Test/test_conic.jl b/src/Test/test_conic.jl index 6e898c5236..952b1ee15c 100644 --- a/src/Test/test_conic.jl +++ b/src/Test/test_conic.jl @@ -4753,7 +4753,7 @@ function _test_conic_PositiveSemidefiniteCone_helper( atol = config.atol rtol = config.rtol square = psdcone == MOI.PositiveSemidefiniteConeSquare - scaled = psdcone == MOI.ScaledPositiveSemidefiniteConeTriangle + scaled = psdcone == MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle} @requires MOI.supports_incremental_interface(model) @requires MOI.supports( model, @@ -4989,7 +4989,7 @@ function test_conic_ScaledPositiveSemidefiniteConeTriangle_VectorAffineFunction( _test_conic_PositiveSemidefiniteCone_helper( model, false, - MOI.ScaledPositiveSemidefiniteConeTriangle, + MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle}, config, ) return @@ -5009,7 +5009,7 @@ function setup_test( ones(T, 3), ( MOI.VectorAffineFunction{T}, - MOI.ScaledPositiveSemidefiniteConeTriangle, + MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle}, ) => [T[1, -sqrt(T(2)), 1]], (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}) => T[2], ), diff --git a/src/Utilities/Utilities.jl b/src/Utilities/Utilities.jl index 7a626ec7e6..c553545f4f 100644 --- a/src/Utilities/Utilities.jl +++ b/src/Utilities/Utilities.jl @@ -79,6 +79,7 @@ include("universalfallback.jl") include("print.jl") include("penalty_relaxation.jl") include("lazy_iterators.jl") +include("set_dot.jl") include("distance_to_set.jl") diff --git a/src/Utilities/matrix_of_constraints.jl b/src/Utilities/matrix_of_constraints.jl index 1f5e0418ca..10bdac499f 100644 --- a/src/Utilities/matrix_of_constraints.jl +++ b/src/Utilities/matrix_of_constraints.jl @@ -561,15 +561,14 @@ function set_with_dimension(::Type{S}, dim) where {S<:MOI.AbstractVectorSet} return S(dim) end +function set_with_dimension(::Type{MOI.Scaled{S}}, dim) where {S} + return MOI.Scaled(set_with_dimension(S, dim)) +end + function set_with_dimension( ::Type{S}, dim, -) where { - S<:Union{ - MOI.AbstractSymmetricMatrixSetTriangle, - MOI.ScaledPositiveSemidefiniteConeTriangle, - }, -} +) where {S<:MOI.AbstractSymmetricMatrixSetTriangle} side_dimension = side_dimension_for_vectorized_dimension(dim) return S(side_dimension) end diff --git a/src/Utilities/model.jl b/src/Utilities/model.jl index 8881039983..b385b9d376 100644 --- a/src/Utilities/model.jl +++ b/src/Utilities/model.jl @@ -792,7 +792,7 @@ const LessThanIndicatorZero{T} = MOI.PositiveSemidefiniteConeTriangle, MOI.PositiveSemidefiniteConeSquare, MOI.HermitianPositiveSemidefiniteConeTriangle, - MOI.ScaledPositiveSemidefiniteConeTriangle, + MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle}, MOI.RootDetConeTriangle, MOI.RootDetConeSquare, MOI.LogDetConeTriangle, diff --git a/src/Utilities/results.jl b/src/Utilities/results.jl index b56820d9b1..8c85b3ed33 100644 --- a/src/Utilities/results.jl +++ b/src/Utilities/results.jl @@ -488,143 +488,3 @@ function get_fallback( func = MOI.get(model, MOI.ConstraintFunction(), ci) return variable_dual(model, attr, ci, func) end - -# Scalar product. Any vector set defined that does not use the standard scalar -# product between vectors of ``R^n`` should redefine `set_dot` and -# `dot_coefficients`. - -""" - set_dot(x::AbstractVector, y::AbstractVector, set::AbstractVectorSet) - -Return the scalar product between a vector `x` of the set `set` and a vector -`y` of the dual of the set `s`. -""" -function set_dot(x::AbstractVector, y::AbstractVector, ::MOI.AbstractVectorSet) - return dot(x, y) -end - -""" - set_dot(x, y, set::AbstractScalarSet) - -Return the scalar product between a number `x` of the set `set` and a number -`y` of the dual of the set `s`. -""" -set_dot(x, y, ::MOI.AbstractScalarSet) = dot(x, y) - -function triangle_dot( - x::AbstractVector{S}, - y::AbstractVector{T}, - dim::Int, - offset::Int, -) where {S,T} - U = MA.promote_operation(MA.add_mul, S, T) - result = zero(U) - k = offset - for i in 1:dim - for j in 1:i - k += 1 - if i == j - result = MA.add_mul!!(result, x[k], y[k]) - else - result = MA.add_mul!!(result, 2, x[k], y[k]) - end - end - end - return result -end - -function set_dot( - x::AbstractVector, - y::AbstractVector, - set::MOI.AbstractSymmetricMatrixSetTriangle, -) - return triangle_dot(x, y, MOI.side_dimension(set), 0) -end - -function MOI.Utilities.set_dot( - x::AbstractVector, - y::AbstractVector, - set::MOI.HermitianPositiveSemidefiniteConeTriangle, -) - sym = MOI.PositiveSemidefiniteConeTriangle(set.side_dimension) - result = set_dot(x, y, sym) - for k in (MOI.dimension(sym)+1):MOI.dimension(set) - result = MA.add_mul!!(result, 2, x[k], y[k]) - end - return result -end - -function set_dot( - x::AbstractVector, - y::AbstractVector, - set::MOI.RootDetConeTriangle, -) - return x[1] * y[1] + triangle_dot(x, y, set.side_dimension, 1) -end - -function set_dot( - x::AbstractVector, - y::AbstractVector, - set::MOI.LogDetConeTriangle, -) - return x[1] * y[1] + x[2] * y[2] + triangle_dot(x, y, set.side_dimension, 2) -end - -""" - dot_coefficients(a::AbstractVector, set::AbstractVectorSet) - -Return the vector `b` such that for all vector `x` of the set `set`, -`set_dot(b, x, set)` is equal to `dot(a, x)`. -""" -function dot_coefficients(a::AbstractVector, ::MOI.AbstractVectorSet) - return a -end - -function triangle_coefficients!( - b::AbstractVector{T}, - dim::Int, - offset::Int, -) where {T} - k = offset - for i in 1:dim - for j in 1:i - k += 1 - if i != j - b[k] /= 2 - end - end - end -end - -function dot_coefficients( - a::AbstractVector, - set::MOI.AbstractSymmetricMatrixSetTriangle, -) - b = copy(a) - triangle_coefficients!(b, MOI.side_dimension(set), 0) - return b -end - -function dot_coefficients( - a::AbstractVector, - set::MOI.HermitianPositiveSemidefiniteConeTriangle, -) - sym = MOI.PositiveSemidefiniteConeTriangle(set.side_dimension) - b = dot_coefficients(a, sym) - for k in (MOI.dimension(sym)+1):MOI.dimension(set) - b[k] /= 2 - end - return b -end - -function dot_coefficients(a::AbstractVector, set::MOI.RootDetConeTriangle) - b = copy(a) - triangle_coefficients!(b, set.side_dimension, 1) - return b -end - -function dot_coefficients(a::AbstractVector, set::MOI.LogDetConeTriangle) - b = copy(a) - triangle_coefficients!(b, set.side_dimension, 2) - return b -end diff --git a/src/Utilities/set_dot.jl b/src/Utilities/set_dot.jl new file mode 100644 index 0000000000..828dd349d0 --- /dev/null +++ b/src/Utilities/set_dot.jl @@ -0,0 +1,352 @@ +# 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. + +# Scalar product. Any vector set defined that does not use the standard scalar +# product between vectors of ``R^n`` should redefine `set_dot` and +# `dot_coefficients`. + +""" + set_dot(x::AbstractVector, y::AbstractVector, set::AbstractVectorSet) + +Return the scalar product between a vector `x` of the set `set` and a vector +`y` of the dual of the set `s`. +""" +function set_dot(x::AbstractVector, y::AbstractVector, ::MOI.AbstractVectorSet) + return dot(x, y) +end + +""" + set_dot(x, y, set::AbstractScalarSet) + +Return the scalar product between a number `x` of the set `set` and a number +`y` of the dual of the set `s`. +""" +set_dot(x, y, ::MOI.AbstractScalarSet) = dot(x, y) + +function triangle_dot( + x::AbstractVector{S}, + y::AbstractVector{T}, + dim::Int, + offset::Int, +) where {S,T} + U = MA.promote_operation(MA.add_mul, S, T) + result = zero(U) + k = offset + for i in 1:dim + for j in 1:i + k += 1 + if i == j + result = MA.add_mul!!(result, x[k], y[k]) + else + result = MA.add_mul!!(result, 2, x[k], y[k]) + end + end + end + return result +end + +function set_dot( + x::AbstractVector, + y::AbstractVector, + set::MOI.AbstractSymmetricMatrixSetTriangle, +) + return triangle_dot(x, y, MOI.side_dimension(set), 0) +end + +function MOI.Utilities.set_dot( + x::AbstractVector, + y::AbstractVector, + set::MOI.HermitianPositiveSemidefiniteConeTriangle, +) + sym = MOI.PositiveSemidefiniteConeTriangle(set.side_dimension) + I = Base.OneTo(MOI.dimension(sym)) + result = set_dot(view(x, I), view(y, I), sym) + for k in (MOI.dimension(sym)+1):MOI.dimension(set) + result = MA.add_mul!!(result, 2, x[k], y[k]) + end + return result +end + +function set_dot( + x::AbstractVector, + y::AbstractVector, + set::MOI.RootDetConeTriangle, +) + return x[1] * y[1] + triangle_dot(x, y, set.side_dimension, 1) +end + +function set_dot( + x::AbstractVector, + y::AbstractVector, + set::MOI.LogDetConeTriangle, +) + return x[1] * y[1] + x[2] * y[2] + triangle_dot(x, y, set.side_dimension, 2) +end + +""" + dot_coefficients(a::AbstractVector, set::AbstractVectorSet) + +Return the vector `b` such that for all vector `x` of the set `set`, +`set_dot(b, x, set)` is equal to `dot(a, x)`. +""" +function dot_coefficients(a::AbstractVector, ::MOI.AbstractVectorSet) + return a +end + +function triangle_coefficients!( + b::AbstractVector{T}, + dim::Int, + offset::Int, +) where {T} + k = offset + for i in 1:dim + for j in 1:i + k += 1 + if i != j + b[k] /= 2 + end + end + end +end + +function dot_coefficients( + a::AbstractVector, + set::MOI.AbstractSymmetricMatrixSetTriangle, +) + b = copy(a) + triangle_coefficients!(b, MOI.side_dimension(set), 0) + return b +end + +function dot_coefficients( + a::AbstractVector, + set::MOI.HermitianPositiveSemidefiniteConeTriangle, +) + sym = MOI.PositiveSemidefiniteConeTriangle(set.side_dimension) + b = dot_coefficients(a, sym) + for k in (MOI.dimension(sym)+1):MOI.dimension(set) + b[k] /= 2 + end + return b +end + +function dot_coefficients(a::AbstractVector, set::MOI.RootDetConeTriangle) + b = copy(a) + triangle_coefficients!(b, set.side_dimension, 1) + return b +end + +function dot_coefficients(a::AbstractVector, set::MOI.LogDetConeTriangle) + b = copy(a) + triangle_coefficients!(b, set.side_dimension, 2) + return b +end + +# For `SetDotScalingVector`, we would like to compute the dot product +# of canonical vectors in O(1) instead of O(n) +# See https://github.com/jump-dev/Dualization.jl/pull/135 + +""" + struct ZeroVector{T} <: AbstractVector{T} + n::Int + end + +Vector of length `n` with only zero elements. +""" +struct ZeroVector{T} <: AbstractVector{T} + n::Int +end + +Base.eltype(::Type{ZeroVector{T}}) where {T} = T + +Base.length(v::ZeroVector) = v.n + +Base.size(v::ZeroVector) = (v.n,) + +function Base.getindex(::ZeroVector{T}, ::Integer) where {T} + return zero(T) +end + +""" + struct CanonicalVector{T} <: AbstractVector{T} + index::Int + n::Int + end + +Vector of length `n` with only zero elements except at index `index` where +the element is one. +""" +struct CanonicalVector{T} <: AbstractVector{T} + index::Int + n::Int +end + +Base.eltype(::Type{CanonicalVector{T}}) where {T} = T + +Base.length(v::CanonicalVector) = v.n + +Base.size(v::CanonicalVector) = (v.n,) + +function Base.getindex(v::CanonicalVector{T}, i::Integer) where {T} + return convert(T, i == v.index) +end + +function Base.view(v::CanonicalVector{T}, I::AbstractUnitRange) where {T} + if v.index in I + return CanonicalVector{T}(v.index - first(I) + 1, length(I)) + else + return ZeroVector{T}(length(I)) + end +end + +# This is much faster than the default implementation that goes +# through all entries even if only one is nonzero. +function LinearAlgebra.dot( + x::CanonicalVector{T}, + y::CanonicalVector{T}, +) where {T} + return convert(T, x.index == y.index) +end + +function triangle_dot( + x::CanonicalVector{T}, + y::CanonicalVector{T}, + ::Int, + offset::Int, +) where {T} + if x.index != y.index || x.index <= offset + return zero(T) + elseif is_diagonal_vectorized_index(x.index - offset) + return one(T) + else + return 2one(T) + end +end + +function _set_dot(i::Integer, s::MOI.AbstractVectorSet, T::Type) + vec = CanonicalVector{T}(i, MOI.dimension(s)) + return set_dot(vec, vec, s) +end + +function _set_dot(::Integer, ::MOI.AbstractScalarSet, T::Type) + return one(T) +end + +""" + struct SetDotScalingVector{T,S<:MOI.AbstractSet} <: AbstractVector{T} + set::S + len::Int + end + +Vector `s` of scaling for the entries of the vectorized form of +a vector `x` in `set` and `y` in `MOI.dual_set(set)` such that +`MOI.Utilities.set_dot(x, y) = LinearAlgebra.dot(s .* x, s .* y)`. + +## Examples + +Combined with `LinearAlgebra`, this vector can be used to scale +a [`MOI.AbstractVectorFunction`](@ref). + +``` +julia> import MathOptInterface as MOI + +julia> model = MOI.Utilities.Model{Float64}() +MOIU.Model{Float64} + +julia> x = MOI.add_variables(model, 3); + +julia> func = MOI.VectorOfVariables(x) +┌ ┐ +│MOI.VariableIndex(1)│ +│MOI.VariableIndex(2)│ +│MOI.VariableIndex(3)│ +└ ┘ + +julia> set = MOI.PositiveSemidefiniteConeTriangle(2) +MathOptInterface.PositiveSemidefiniteConeTriangle(2) + +julia> MOI.add_constraint(model, func, MOI.Scaled(set)) +MathOptInterface.ConstraintIndex{MathOptInterface.VectorOfVariables, MathOptInterface.Scaled{MathOptInterface.PositiveSemidefiniteConeTriangle}}(1) + +julia> a = MOI.Utilities.SetDotScalingVector{Float64}(set) +3-element MathOptInterface.Utilities.SetDotScalingVector{Float64, MathOptInterface.PositiveSemidefiniteConeTriangle}: + 1.0 + 1.4142135623730951 + 1.0 + +julia> using LinearAlgebra + +julia> MOI.Utilities.operate(*, Float64, Diagonal(a), func) +┌ ┐ +│0.0 + 1.0 MOI.VariableIndex(1) │ +│0.0 + 1.4142135623730951 MOI.VariableIndex(2)│ +│0.0 + 1.0 MOI.VariableIndex(3) │ +└ ┘ + +julia> MOI.Utilities.operate(*, Float64, Diagonal(a), ones(3)) +3-element Vector{Float64}: + 1.0 + 1.4142135623730951 + 1.0 +``` +""" +struct SetDotScalingVector{T,S<:MOI.AbstractSet} <: AbstractVector{T} + set::S +end + +function SetDotScalingVector{T}(s::MOI.AbstractSet) where {T} + return SetDotScalingVector{T,typeof(s)}(s) +end + +function Base.getindex(s::SetDotScalingVector{T}, i::Base.Integer) where {T} + return sqrt(_set_dot(i, s.set, T)) +end + +Base.size(x::SetDotScalingVector) = (MOI.dimension(x.set),) + +function symmetric_matrix_scaling_vector(::Type{T}, n) where {T} + d = side_dimension_for_vectorized_dimension(n) + set = MOI.PositiveSemidefiniteConeTriangle(d) + return SetDotScalingVector{T}(set) +end + +function symmetric_matrix_inverse_scaling_vector(::Type{T}, n) where {T} + return lazy_map(T, inv, symmetric_matrix_scaling_vector(T, n)) +end + +""" + struct SymmetricMatrixScalingVector{T} <: AbstractVector{T} + no_scaling::T + scaling::T + len::Int + end + +Vector of scaling for the entries of the vectorized form of +a symmetric matrix. The values `no_scaling` and `scaling` +are stored in the `struct` to avoid creating a new one for each entry. + +!!! warning + This type is deprecated, use `SetDotScalingVector` instead. +""" +struct SymmetricMatrixScalingVector{T} <: AbstractVector{T} + scaling::T + no_scaling::T + len::Int +end + +function SymmetricMatrixScalingVector{T}(scaling::T, len::Int) where {T} + return SymmetricMatrixScalingVector{T}(scaling, one(T), len) +end + +function Base.getindex(s::SymmetricMatrixScalingVector, i::Base.Integer) + if is_diagonal_vectorized_index(i) + return s.no_scaling + else + return s.scaling + end +end + +Base.size(x::SymmetricMatrixScalingVector) = (x.len,) diff --git a/src/Utilities/sets.jl b/src/Utilities/sets.jl index 84ec1d61b3..1875d87e7a 100644 --- a/src/Utilities/sets.jl +++ b/src/Utilities/sets.jl @@ -158,45 +158,6 @@ function trimap(row::Integer, column::Integer) return div((row - 1) * row, 2) + column end -""" - struct SymmetricMatrixScalingVector{T} <: AbstractVector{T} - no_scaling::T - scaling::T - len::Int - end - -Vector of scaling for the entries of the vectorized form of -a symmetric matrix. The values `no_scaling` and `scaling` -are stored in the `struct` to avoid creating a new one for each entry. -""" -struct SymmetricMatrixScalingVector{T} <: AbstractVector{T} - scaling::T - no_scaling::T - len::Int -end - -function SymmetricMatrixScalingVector{T}(scaling::T, len::Int) where {T} - return SymmetricMatrixScalingVector{T}(scaling, one(T), len) -end - -function symmetric_matrix_scaling_vector(::Type{T}, len::Int) where {T} - return SymmetricMatrixScalingVector{T}(sqrt(T(2)), len) -end - -function symmetric_matrix_inverse_scaling_vector(::Type{T}, len::Int) where {T} - return SymmetricMatrixScalingVector{T}(inv(sqrt(T(2))), len) -end - -function Base.getindex(s::SymmetricMatrixScalingVector, i::Base.Integer) - if is_diagonal_vectorized_index(i) - return s.no_scaling - else - return s.scaling - end -end - -Base.size(x::SymmetricMatrixScalingVector) = (x.len,) - similar_type(::Type{<:MOI.LessThan}, ::Type{T}) where {T} = MOI.LessThan{T} function similar_type(::Type{<:MOI.GreaterThan}, ::Type{T}) where {T} diff --git a/src/sets.jl b/src/sets.jl index de1d686f9d..b321c0e65b 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -1139,59 +1139,6 @@ dual_set_type(::Type{NormNuclearCone}) = NormSpectralCone dimension(s::NormNuclearCone) = 1 + s.row_dim * s.column_dim -""" - ScaledPositiveSemidefiniteConeTriangle(side_dimension::Int) <: AbstractVectorSet - -The *scaled* (vectorized) cone of symmetric positive semidefinite matrices, with -non-negative `side_dimension` rows and columns. - -Compared to the [`PositiveSemidefiniteConeTriangle`](@ref), the off-diagonal -entries are scaled by `√2`. Thanks to this scaling, [`Utilities.set_dot`](@ref) -is the simply the sum of the pairwise product, while for -[`PositiveSemidefiniteConeTriangle`](@ref), the pairwise product additionally -have to be multiplied by `2`. - -## Example - -```jldoctest -julia> import MathOptInterface as MOI - -julia> model = MOI.Utilities.Model{Float64}() -MOIU.Model{Float64} - -julia> x = MOI.add_variables(model, 3); - -julia> MOI.add_constraint( - model, - MOI.VectorOfVariables(x), - MOI.ScaledPositiveSemidefiniteConeTriangle(2), - ) -MathOptInterface.ConstraintIndex{MathOptInterface.VectorOfVariables, MathOptInterface.ScaledPositiveSemidefiniteConeTriangle}(1) -``` -""" -struct ScaledPositiveSemidefiniteConeTriangle <: AbstractVectorSet - side_dimension::Int - function ScaledPositiveSemidefiniteConeTriangle( - side_dimension::Base.Integer, - ) - if !(side_dimension >= 0) - throw( - DimensionMismatch( - "Side dimension of PositiveSemidefiniteConeTriangle must " * - "be >= 0, not $(side_dimension).", - ), - ) - end - return new(side_dimension) - end -end - -dual_set(s::ScaledPositiveSemidefiniteConeTriangle) = copy(s) - -function dual_set_type(::Type{ScaledPositiveSemidefiniteConeTriangle}) - return ScaledPositiveSemidefiniteConeTriangle -end - """ abstract type AbstractSymmetricMatrixSetTriangle <: AbstractVectorSet end @@ -1289,12 +1236,9 @@ products we have """ abstract type AbstractSymmetricMatrixSetTriangle <: AbstractVectorSet end -function dimension( - set::Union{ - AbstractSymmetricMatrixSetTriangle, - ScaledPositiveSemidefiniteConeTriangle, - }, -) +is_set_dot_scaled(::Type{<:AbstractSymmetricMatrixSetTriangle}) = true + +function dimension(set::AbstractSymmetricMatrixSetTriangle) d = side_dimension(set) return div(d * (d + 1), 2) end @@ -1332,33 +1276,6 @@ abstract type AbstractSymmetricMatrixSetSquare <: AbstractVectorSet end dimension(set::AbstractSymmetricMatrixSetSquare) = side_dimension(set)^2 -""" - side_dimension( - set::Union{ - ScaledPositiveSemidefiniteConeTriangle, - AbstractSymmetricMatrixSetTriangle, - AbstractSymmetricMatrixSetSquare, - }, - ) - -Side dimension of the matrices in `set`. - -## Convention - -By convention, the side dimension should be stored in the `side_dimension` -field. If this is not the case for a subtype of [`AbstractSymmetricMatrixSetTriangle`](@ref), -you must implement this method. -""" -function side_dimension( - set::Union{ - ScaledPositiveSemidefiniteConeTriangle, - AbstractSymmetricMatrixSetTriangle, - AbstractSymmetricMatrixSetSquare, - }, -) - return set.side_dimension -end - """ triangular_form(S::Type{<:AbstractSymmetricMatrixSetSquare}) triangular_form(set::AbstractSymmetricMatrixSetSquare) @@ -1509,6 +1426,33 @@ function dimension(set::HermitianPositiveSemidefiniteConeTriangle) return real_nnz + imag_nnz end +""" + side_dimension( + set::Union{ + AbstractSymmetricMatrixSetTriangle, + AbstractSymmetricMatrixSetSquare, + HermitianPositiveSemidefiniteConeTriangle, + }, + ) + +Side dimension of the matrices in `set`. + +## Convention + +By convention, the side dimension should be stored in the `side_dimension` +field. If this is not the case for a subtype of [`AbstractSymmetricMatrixSetTriangle`](@ref), +or [`AbstractSymmetricMatrixSetSquare`](@ref) you must implement this method. +""" +function side_dimension( + set::Union{ + AbstractSymmetricMatrixSetTriangle, + AbstractSymmetricMatrixSetSquare, + HermitianPositiveSemidefiniteConeTriangle, + }, +) + return set.side_dimension +end + """ LogDetConeTriangle(side_dimension::Int) @@ -1735,6 +1679,148 @@ function triangular_form(set::RootDetConeSquare) return RootDetConeTriangle(set.side_dimension) end +""" + is_set_dot_scaled(::Type{<:AbstractVectorFunction}) + +Return whether [`Utilities.set_dot(x, y)`](@ref Utilities.set_dot) is equivalent +to `x' * Diagonal(s) * y` for some scaling vector `s`. This scaling vector `s` +can be obtained with [`Utilities.SetDotScalingVector`](@ref). + +!!! note + For such sets `S`, we `Diagonal(s) * S` is [`MOI.Scaled(S)`](@ref Scaled). + This linear relationship between the two sets allows briding between them + with [`Bridges.Constraint.SetDotScalingBridge`](@ref) and + [`Bridges.Constraint.SetDotInverseScalingBridge`](@ref). This scaling vector + `s` is also used by Dualization.jl to compute the dual. +""" +is_set_dot_scaled(::Type{<:AbstractVectorSet}) = false + +function is_set_dot_scaled( + ::Type{ + <:Union{ + AbstractSymmetricMatrixSetTriangle, + HermitianPositiveSemidefiniteConeTriangle, + LogDetConeTriangle, + RootDetConeTriangle, + }, + }, +) + return true +end + +""" + struct Scaled{S<:AbstractVectorSet} <: AbstractVectorSet + set::S + end + +Given a vector ``a \\in \\mathbb{R}^d`` and a `set` representing the set +``\\mathcal{S} \\in \\mathbb{R}^d`` such that [`Utilities.set_dot`](@ref) for +``x \\in \\mathcal{S}`` and ``y \\in \\mathcal{S}^*`` is +```math +\\sum_{i=1}^d a_i x_i y_i +``` +the set `Scaled(set)` is defined as +```math +\\{ (\\sqrt{a_1} x_1, \\sqrt{a_2} x_2, \\ldots, \\sqrt{a_d} x_d) : x \\in S \\} +``` + +## Example + +This can be used to scale a vector of numbers +```jldoctest scaling +julia> import MathOptInterface as MOI + +julia> set = MOI.PositiveSemidefiniteConeTriangle(2) +MathOptInterface.PositiveSemidefiniteConeTriangle(2) + +julia> a = MOI.Utilities.SetDotScalingVector{Float64}(set) +3-element MathOptInterface.Utilities.SetDotScalingVector{Float64, MathOptInterface.PositiveSemidefiniteConeTriangle}: + 1.0 + 1.4142135623730951 + 1.0 + +julia> using LinearAlgebra + +julia> MOI.Utilities.operate(*, Float64, Diagonal(a), ones(3)) +3-element Vector{Float64}: + 1.0 + 1.4142135623730951 + 1.0 +``` + +It can be also used to scale a vector of function +```jldoctest scaling +julia> model = MOI.Utilities.Model{Float64}() +MOIU.Model{Float64} + +julia> x = MOI.add_variables(model, 3); + +julia> func = MOI.VectorOfVariables(x) +┌ ┐ +│MOI.VariableIndex(1)│ +│MOI.VariableIndex(2)│ +│MOI.VariableIndex(3)│ +└ ┘ + +julia> set = MOI.PositiveSemidefiniteConeTriangle(2) +MathOptInterface.PositiveSemidefiniteConeTriangle(2) + +julia> MOI.Utilities.operate(*, Float64, Diagonal(a), func) +┌ ┐ +│0.0 + 1.0 MOI.VariableIndex(1) │ +│0.0 + 1.4142135623730951 MOI.VariableIndex(2)│ +│0.0 + 1.0 MOI.VariableIndex(3) │ +└ ┘ +``` +""" +struct Scaled{S<:AbstractVectorSet} <: AbstractVectorSet + set::S +end + +dimension(s::Scaled) = dimension(s.set) + +side_dimension(s::Scaled) = side_dimension(s.set) + +dual_set(s::Scaled) = Scaled(dual_set(s.set)) + +dual_set_type(::Type{Scaled{S}}) where {S} = Scaled{dual_set_type(S)} + +""" + const ScaledPositiveSemidefiniteConeTriangle = Scaled{PositiveSemidefiniteConeTriangle} + +The [`Scaled`](@ref) (vectorized) cone of symmetric positive semidefinite matrices, with +non-negative `side_dimension` rows and columns. + +Compared to the [`PositiveSemidefiniteConeTriangle`](@ref), the off-diagonal +entries are scaled by `√2`. Thanks to this scaling, [`Utilities.set_dot`](@ref) +is the simply the sum of the pairwise product, while for +[`PositiveSemidefiniteConeTriangle`](@ref), the pairwise product additionally +have to be multiplied by `2`. + +!!! note + Using `MOI.ScaledPositiveSemidefiniteConeTriangle(d)` is deprecated, use + `MOI.Scaled(MOI.PositiveSemidefiniteConeTriangle(d))`. +""" +const ScaledPositiveSemidefiniteConeTriangle = + Scaled{PositiveSemidefiniteConeTriangle} + +# TODO remove in MOI v2 +function Scaled{PositiveSemidefiniteConeTriangle}(side_dimension::Int) + return Scaled(PositiveSemidefiniteConeTriangle(side_dimension)) +end + +# TODO remove in MOI v2 +function Base.getproperty( + set::Scaled{PositiveSemidefiniteConeTriangle}, + f::Symbol, +) + if f == :side_dimension + return getproperty(set.set, f) + else + return getfield(set, f) + end +end + """ SOS1{T<:Real}(weights::Vector{T}) diff --git a/test/Bridges/Constraint/symmetric_matrix_scaling.jl b/test/Bridges/Constraint/set_dot_scaling.jl similarity index 84% rename from test/Bridges/Constraint/symmetric_matrix_scaling.jl rename to test/Bridges/Constraint/set_dot_scaling.jl index 1c2a5693cf..134c8bafa9 100644 --- a/test/Bridges/Constraint/symmetric_matrix_scaling.jl +++ b/test/Bridges/Constraint/set_dot_scaling.jl @@ -4,7 +4,7 @@ # 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 TestConstraintSymmetricMatrixScaling +module TestConstraintSetDotScaling using Test @@ -23,7 +23,7 @@ end function test_scaling() MOI.Bridges.runtests( - MOI.Bridges.Constraint.SymmetricMatrixScalingBridge, + MOI.Bridges.Constraint.SetDotScalingBridge, """ variables: x, y, z [x, 1.0 * y, z] in PositiveSemidefiniteConeTriangle(2) @@ -38,7 +38,7 @@ end function test_scaling_vector_of_variables() MOI.Bridges.runtests( - MOI.Bridges.Constraint.SymmetricMatrixScalingBridge, + MOI.Bridges.Constraint.SetDotScalingBridge, """ variables: x, y, z [x, y, z] in PositiveSemidefiniteConeTriangle(2) @@ -53,7 +53,7 @@ end function test_scaling_quadratic() MOI.Bridges.runtests( - MOI.Bridges.Constraint.SymmetricMatrixScalingBridge, + MOI.Bridges.Constraint.SetDotScalingBridge, """ variables: x, y, z [x, 1.0 * y * y + 1.0 * y + 3.0, z] in PositiveSemidefiniteConeTriangle(2) @@ -68,7 +68,7 @@ end function test_inverse_scaling() MOI.Bridges.runtests( - MOI.Bridges.Constraint.SymmetricMatrixInverseScalingBridge, + MOI.Bridges.Constraint.SetDotInverseScalingBridge, """ variables: x, y, z [x, √2 * y, z] in ScaledPositiveSemidefiniteConeTriangle(2) @@ -83,7 +83,7 @@ end function test_inverse_scaling_vector_of_variables() MOI.Bridges.runtests( - MOI.Bridges.Constraint.SymmetricMatrixInverseScalingBridge, + MOI.Bridges.Constraint.SetDotInverseScalingBridge, """ variables: x, y, z [x, y, z] in ScaledPositiveSemidefiniteConeTriangle(2) @@ -98,7 +98,7 @@ end function test_inverse_scaling_quadratic() MOI.Bridges.runtests( - MOI.Bridges.Constraint.SymmetricMatrixInverseScalingBridge, + MOI.Bridges.Constraint.SetDotInverseScalingBridge, """ variables: x, y, z [x, √2 * y * y, z] in ScaledPositiveSemidefiniteConeTriangle(2) @@ -113,4 +113,4 @@ end end # module -TestConstraintSymmetricMatrixScaling.runtests() +TestConstraintSetDotScaling.runtests() diff --git a/test/Bridges/lazy_bridge_optimizer.jl b/test/Bridges/lazy_bridge_optimizer.jl index 319b3742f8..4695af4451 100644 --- a/test/Bridges/lazy_bridge_optimizer.jl +++ b/test/Bridges/lazy_bridge_optimizer.jl @@ -1324,7 +1324,7 @@ function _test_continuous_conic_with_NoVariableModel(T) "Bridge graph with ", "constrained variables in `MOI.NormCone` are not supported", "`MOI.VectorOfVariables`-in-`MOI.SecondOrderCone` constraints are bridged", - "`MOI.VectorOfVariables`-in-`MOI.ScaledPositiveSemidefiniteConeTriangle` constraints are not supported", + "`MOI.VectorOfVariables`-in-`MOI.Scaled{MOI.PositiveSemidefiniteConeTriangle}` constraints are not supported", ) @test occursin(needle, graph) end diff --git a/test/Utilities/matrix_of_constraints.jl b/test/Utilities/matrix_of_constraints.jl index b803a473d3..275e0aa72f 100644 --- a/test/Utilities/matrix_of_constraints.jl +++ b/test/Utilities/matrix_of_constraints.jl @@ -607,6 +607,12 @@ function test_set_with_dimension() @test_throws AssertionError MOI.Utilities.set_with_dimension(S, 2) @test_throws AssertionError MOI.Utilities.set_with_dimension(S, 4) end + S = MOI.PositiveSemidefiniteConeTriangle + set = MOI.Utilities.set_with_dimension(MOI.Scaled{S}, 3) + @test MOI.side_dimension(set) == 2 + S = MOI.HermitianPositiveSemidefiniteConeTriangle + set = MOI.Utilities.set_with_dimension(MOI.Scaled{S}, 4) + @test MOI.side_dimension(set) == 2 return end diff --git a/test/Utilities/set_dot.jl b/test/Utilities/set_dot.jl new file mode 100644 index 0000000000..7076f56fb5 --- /dev/null +++ b/test/Utilities/set_dot.jl @@ -0,0 +1,53 @@ +# 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 TestSetDot + +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 + +function test_set_dot(T = Int) + @test MOI.Utilities._set_dot(1, MOI.ZeroOne(), T) == 1 + set = MOI.PositiveSemidefiniteConeTriangle(2) + n = 3 + a = MOI.Utilities.ZeroVector{T}(n) + @test eltype(a) == T + @test length(a) == n + @test size(a) == (n,) + for i in 1:n + b = MOI.Utilities.CanonicalVector{T}(i, n) + @test eltype(b) == T + @test length(b) == n + @test size(b) == (n,) + c = MOI.Utilities.CanonicalVector{T}(mod1(i + 1, 3), n) + @test a[i] == 0 + @test b[i] == 1 + @test c[i] == 0 + @test iszero(MOI.Utilities.set_dot(a, b, set)) + @test iszero(MOI.Utilities.set_dot(b, a, set)) + @test iszero(MOI.Utilities.set_dot(b, c, set)) + @test iszero(MOI.Utilities.set_dot(c, b, set)) + expected = i == 2 ? 2 : 1 + @test MOI.Utilities.set_dot(b, b, set) == expected + @test MOI.Utilities._set_dot(i, set, T) == expected + end +end + +end + +TestSetDot.runtests() diff --git a/test/Utilities/sets.jl b/test/Utilities/sets.jl index e059475764..bb7102baaf 100644 --- a/test/Utilities/sets.jl +++ b/test/Utilities/sets.jl @@ -197,23 +197,29 @@ function test_trimap() return end -function test_symmetric_matrix_scaling() - n = 10 +function test_set_dot_scaling(n = 10) N = div(n * (n + 1), 2) + M = N + div(n * (n - 1), 2) v = MOI.Utilities.SymmetricMatrixScalingVector{Float64}(1.5, 0.5, N) s = MOI.Utilities.symmetric_matrix_scaling_vector(Float64, N) s32 = MOI.Utilities.symmetric_matrix_scaling_vector(Float32, N) is = MOI.Utilities.symmetric_matrix_inverse_scaling_vector(Float64, N) is32 = MOI.Utilities.symmetric_matrix_inverse_scaling_vector(Float32, N) + hpsd = MOI.HermitianPositiveSemidefiniteConeTriangle(n) + hermitian = MOI.Utilities.SetDotScalingVector{Float64}(hpsd) k = 0 + imag_k = 0 for j in 1:n for i in 1:(j-1) k += 1 + imag_k += 1 @test v[k] == 1.5 @test s[k] == √2 @test s32[k] == √Float32(2) @test is[k] == inv(√2) @test is32[k] == inv(√Float32(2)) + @test hermitian[k] == √2 + @test hermitian[N+imag_k] == √2 end k += 1 @test v[k] == 0.5 @@ -221,6 +227,7 @@ function test_symmetric_matrix_scaling() @test isone(s32[k]) @test isone(is[k]) @test isone(is32[k]) + @test isone(hermitian[k]) end return end From bf3177f3f164cc50959577fdb7b8a1b0eac92b18 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Wed, 6 Sep 2023 10:04:17 +1200 Subject: [PATCH 2/3] Apply suggestions from code review --- src/Utilities/set_dot.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Utilities/set_dot.jl b/src/Utilities/set_dot.jl index 828dd349d0..2f4b81a45d 100644 --- a/src/Utilities/set_dot.jl +++ b/src/Utilities/set_dot.jl @@ -15,7 +15,7 @@ Return the scalar product between a vector `x` of the set `set` and a vector `y` of the dual of the set `s`. """ function set_dot(x::AbstractVector, y::AbstractVector, ::MOI.AbstractVectorSet) - return dot(x, y) + return LinearAlgebra.dot(x, y) end """ @@ -24,7 +24,7 @@ end Return the scalar product between a number `x` of the set `set` and a number `y` of the dual of the set `s`. """ -set_dot(x, y, ::MOI.AbstractScalarSet) = dot(x, y) +set_dot(x, y, ::MOI.AbstractScalarSet) = LinearAlgebra.dot(x, y) function triangle_dot( x::AbstractVector{S}, @@ -250,7 +250,7 @@ a vector `x` in `set` and `y` in `MOI.dual_set(set)` such that Combined with `LinearAlgebra`, this vector can be used to scale a [`MOI.AbstractVectorFunction`](@ref). -``` +```jldoctest julia> import MathOptInterface as MOI julia> model = MOI.Utilities.Model{Float64}() From 28f9c19c6551a5eebf2c8eeb83f4076802617c54 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 6 Sep 2023 09:58:57 +0200 Subject: [PATCH 3/3] Improve code coverage --- test/Utilities/set_dot.jl | 19 ++++++++++++------- test/Utilities/sets.jl | 3 +++ 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/test/Utilities/set_dot.jl b/test/Utilities/set_dot.jl index 7076f56fb5..969e0442c4 100644 --- a/test/Utilities/set_dot.jl +++ b/test/Utilities/set_dot.jl @@ -23,8 +23,9 @@ end function test_set_dot(T = Int) @test MOI.Utilities._set_dot(1, MOI.ZeroOne(), T) == 1 - set = MOI.PositiveSemidefiniteConeTriangle(2) + psd = MOI.PositiveSemidefiniteConeTriangle(2) n = 3 + soc = MOI.SecondOrderCone(n) a = MOI.Utilities.ZeroVector{T}(n) @test eltype(a) == T @test length(a) == n @@ -38,13 +39,17 @@ function test_set_dot(T = Int) @test a[i] == 0 @test b[i] == 1 @test c[i] == 0 - @test iszero(MOI.Utilities.set_dot(a, b, set)) - @test iszero(MOI.Utilities.set_dot(b, a, set)) - @test iszero(MOI.Utilities.set_dot(b, c, set)) - @test iszero(MOI.Utilities.set_dot(c, b, set)) + for set in [psd, soc] + @test iszero(MOI.Utilities.set_dot(a, b, set)) + @test iszero(MOI.Utilities.set_dot(b, a, set)) + @test iszero(MOI.Utilities.set_dot(b, c, set)) + @test iszero(MOI.Utilities.set_dot(c, b, set)) + end + @test MOI.Utilities.set_dot(b, b, soc) == 1 + @test MOI.Utilities._set_dot(i, soc, T) == 1 expected = i == 2 ? 2 : 1 - @test MOI.Utilities.set_dot(b, b, set) == expected - @test MOI.Utilities._set_dot(i, set, T) == expected + @test MOI.Utilities.set_dot(b, b, psd) == expected + @test MOI.Utilities._set_dot(i, psd, T) == expected end end diff --git a/test/Utilities/sets.jl b/test/Utilities/sets.jl index bb7102baaf..789c29a86a 100644 --- a/test/Utilities/sets.jl +++ b/test/Utilities/sets.jl @@ -201,6 +201,7 @@ function test_set_dot_scaling(n = 10) N = div(n * (n + 1), 2) M = N + div(n * (n - 1), 2) v = MOI.Utilities.SymmetricMatrixScalingVector{Float64}(1.5, 0.5, N) + w = MOI.Utilities.SymmetricMatrixScalingVector{Float64}(1.5, N) s = MOI.Utilities.symmetric_matrix_scaling_vector(Float64, N) s32 = MOI.Utilities.symmetric_matrix_scaling_vector(Float32, N) is = MOI.Utilities.symmetric_matrix_inverse_scaling_vector(Float64, N) @@ -214,6 +215,7 @@ function test_set_dot_scaling(n = 10) k += 1 imag_k += 1 @test v[k] == 1.5 + @test w[k] == 1.5 @test s[k] == √2 @test s32[k] == √Float32(2) @test is[k] == inv(√2) @@ -223,6 +225,7 @@ function test_set_dot_scaling(n = 10) end k += 1 @test v[k] == 0.5 + @test w[k] == 1 @test isone(s[k]) @test isone(s32[k]) @test isone(is[k])