Implement starting values for Variable.Hermitian bridge

src/Bridges/Variable/bridges/hermitian.jl

@@ -61,14 +61,15 @@ 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,
         MOI.PositiveSemidefiniteConeTriangle,
     }
     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

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}

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)