From f0f1cfce53b4675f8c75aef56ad3b04a5f8cfbbb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sun, 29 Oct 2023 22:55:01 +0100 Subject: [PATCH] Implement starting values for Variable.Hermitian bridge --- src/Bridges/Variable/bridges/hermitian.jl | 75 ++++++++++++++++++++--- src/Utilities/sets.jl | 33 +++++++++- test/Utilities/sets.jl | 11 ++++ 3 files changed, 108 insertions(+), 11 deletions(-) diff --git a/src/Bridges/Variable/bridges/hermitian.jl b/src/Bridges/Variable/bridges/hermitian.jl index f5566b4380..8af98dfa47 100644 --- a/src/Bridges/Variable/bridges/hermitian.jl +++ b/src/Bridges/Variable/bridges/hermitian.jl @@ -61,7 +61,7 @@ equality constraints to: * force the lower triangular of the `y` block to be the negative of the upper triangle. """ -struct HermitianToSymmetricPSDBridge{T} <: AbstractBridge +mutable struct HermitianToSymmetricPSDBridge{T} <: AbstractBridge variables::Vector{MOI.VariableIndex} psd::MOI.ConstraintIndex{ MOI.VectorOfVariables, @@ -69,6 +69,7 @@ struct HermitianToSymmetricPSDBridge{T} <: AbstractBridge } n::Int ceq::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}}} + imag_diag_start_set::Bool end const HermitianToSymmetricPSD{T,OT<:MOI.ModelLike} = @@ -87,11 +88,6 @@ function bridge_constrained_variable( ceq = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}}[] k11 = 0 k12 = k22 = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n)) - function X21(i, j) - I, J = j, n + i - k21 = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(J - 1)) + I - return variables[k21] - end for j in 1:n k22 += n for i in 1:j @@ -104,13 +100,14 @@ function bridge_constrained_variable( f_0 = convert(MOI.ScalarAffineFunction{T}, variables[k12]) push!(ceq, MOI.add_constraint(model, f_0, MOI.EqualTo(zero(T)))) else # y_{ij} = -y_{ji} - f_y = MOI.Utilities.operate(+, T, X21(i, j), variables[k12]) + k21 = MOI.Utilities.trimap(j, n + i) + f_y = MOI.Utilities.operate(+, T, variables[k21], variables[k12]) push!(ceq, MOI.add_constraint(model, f_y, MOI.EqualTo(zero(T)))) end end k12 += n end - return HermitianToSymmetricPSDBridge(variables, psd_ci, n, ceq) + return HermitianToSymmetricPSDBridge(variables, psd_ci, n, ceq, false) end function supports_constrained_variable( @@ -311,7 +308,7 @@ end function MOI.get( model::MOI.ModelLike, - attr::MOI.VariablePrimal, + attr::Union{MOI.VariablePrimal,MOI.VariablePrimalStart}, bridge::HermitianToSymmetricPSDBridge{T}, i::MOI.Bridges.IndexInVector, ) where {T} @@ -332,3 +329,63 @@ function unbridged_map( ) where {T} return (_variable(bridge, i) => convert(MOI.ScalarAffineFunction{T}, vi),) end + +function MOI.supports( + ::MOI.ModelLike, + ::MOI.VariablePrimalStart, + ::Type{<:HermitianToSymmetricPSDBridge}, +) + return true +end + +# Set the starting value of the diagonal elements of the imaginary part +function _set_imag_diag_start( + model::MOI.ModelLike, + attr::MOI.VariablePrimalStart, + bridge::HermitianToSymmetricPSDBridge, + value, +) + for i in 1:bridge.n + k = MOI.Utilities.trimap(i, bridge.n + i) + MOI.set(model, attr, bridge.variables[k], value) + end +end + +_minus(a) = -a +_minus(::Nothing) = nothing + +function MOI.set( + model::MOI.ModelLike, + attr::MOI.VariablePrimalStart, + bridge::HermitianToSymmetricPSDBridge, + value, + index::MOI.Bridges.IndexInVector, +) + n = bridge.n + d = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n)) + if index.value > d + # Imaginary part + i, j = MOI.Utilities.inverse_trimap(index.value - d) + # Increment `j` by `1` to account for the zero diagonal + j += 1 + k12 = MOI.Utilities.trimap(i, n + j) + k21 = MOI.Utilities.trimap(j, n + i) + MOI.set(model, attr, bridge.variables[k12], value) + MOI.set(model, attr, bridge.variables[k21], _minus(value)) + else + # Real part + i, j = MOI.Utilities.inverse_trimap(index.value) + k11 = MOI.Utilities.trimap(i, j) + k22 = MOI.Utilities.trimap(n + i, n + j) + MOI.set(model, attr, bridge.variables[k11], value) + MOI.set(model, attr, bridge.variables[k22], value) + end + if isnothing(value) && bridge.imag_diag_start_set + bridge.imag_diag_start_set = false + _set_imag_diag_start(model, attr, bridge, nothing) + elseif !isnothing(value) && !bridge.imag_diag_start_set + bridge.imag_diag_start_set = true + _set_imag_diag_start(model, attr, bridge, zero(value)) + end + return +end diff --git a/src/Utilities/sets.jl b/src/Utilities/sets.jl index 1875d87e7a..b58a02a0e6 100644 --- a/src/Utilities/sets.jl +++ b/src/Utilities/sets.jl @@ -136,8 +136,8 @@ Return the dimension `d` such that """ function side_dimension_for_vectorized_dimension(n::Base.Integer) # We have `d*(d+1)/2 = n` so - # `d² + d - 2n = 0` hence `d = (-1 ± √(1 + 8d)) / 2` - # The integer `√(1 + 8d)` is odd and `√(1 + 8d) - 1` is even. + # `d² + d - 2n = 0` hence `d = (-1 ± √(1 + 8n)) / 2` + # The integer `√(1 + 8n)` is odd and `√(1 + 8n) - 1` is even. # We can drop the `- 1` as `div` already discards it. return div(isqrt(1 + 8n), 2) end @@ -150,6 +150,9 @@ the corresponding element in the triangular representation. This is most useful when mapping between `ConeSquare` and `ConeTriangle` sets, e.g., as part of an [`MOI.AbstractSymmetricMatrixSetTriangle`](@ref) set. + +!!! note + Use [`inverse_trimap`](@ref) for the reverse mapping. """ function trimap(row::Integer, column::Integer) if row < column @@ -158,6 +161,32 @@ function trimap(row::Integer, column::Integer) return div((row - 1) * row, 2) + column end +""" + inverse_trimap(row::Integer, column::Integer) + +Convert between the the linear index of a +[`MathOptInterface.AbstractSymmetricMatrixSetTriangle`] to the row and column +indices of upper triangular part of the corresponding matrix. + +!!! note + Use [`trimap`](@ref) for the reverse mapping. +""" +function inverse_trimap(index::Integer) + # Because `isqrt` and `div` are monotone functions + # because `index` is between `trimap(j - 1, j - 1)` and + # `trimap(j, j)`,` + # `side_dimension_for_vectorized_dimension(index)` is between `j - 1` + # and `j`. + j = side_dimension_for_vectorized_dimension(index) + n = trimap(j, j) + if index <= n + j -= 1 + n = trimap(j, j) + end + i = index - n + return i, j + 1 +end + 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/test/Utilities/sets.jl b/test/Utilities/sets.jl index 789c29a86a..bbbc1b8ca8 100644 --- a/test/Utilities/sets.jl +++ b/test/Utilities/sets.jl @@ -197,6 +197,17 @@ function test_trimap() return end +function test_inverse_trimap(n = 10) + for j in 1:n + for i in 1:j + I, J = MOI.Utilities.inverse_trimap(MOIU.trimap(i, j)) + @test I == i + @test J == j + end + end + return +end + function test_set_dot_scaling(n = 10) N = div(n * (n + 1), 2) M = N + div(n * (n - 1), 2)