[Bridges] Add hermitian bridge for constraints (#1977)

blegat · web-flow · commit ae0ad01a27c3 · 2022-08-15T09:00:07.000+12:00
diff --git a/docs/src/submodules/Bridges/list_of_bridges.md b/docs/src/submodules/Bridges/list_of_bridges.md
@@ -47,6 +47,7 @@ Bridges.Constraint.RelativeEntropyBridge
 Bridges.Constraint.NormSpectralBridge
 Bridges.Constraint.NormNuclearBridge
 Bridges.Constraint.SquareBridge
+Bridges.Constraint.HermitianToSymmetricPSDBridge
 Bridges.Constraint.RootDetBridge
 Bridges.Constraint.LogDetBridge
 Bridges.Constraint.IndicatorActiveOnFalseBridge
diff --git a/src/Bridges/Constraint/Constraint.jl b/src/Bridges/Constraint/Constraint.jl
@@ -50,6 +50,7 @@ include("bridges/soc_to_nonconvex_quad.jl") # do not add these bridges by defaul
 include("bridges/soc_to_psd.jl")
 include("bridges/split_complex_equalto.jl")
 include("bridges/split_complex_zeros.jl")
+include("bridges/hermitian.jl")
 include("bridges/square.jl")
 include("bridges/table.jl")
 include("bridges/vectorize.jl")
diff --git a/src/Bridges/Constraint/bridges/hermitian.jl b/src/Bridges/Constraint/bridges/hermitian.jl
@@ -0,0 +1,289 @@
+# 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.
+
+"""
+    HermitianToSymmetricPSDBridge{T,F,G} <: Bridges.Constraint.AbstractBridge
+
+`HermitianToSymmetricPSDBridge` implements the following reformulation:
+
+  * Hermitian positive semidefinite `n x n` complex matrix to a symmetric
+    positive semidefinite `2n x 2n` real matrix.
+
+See also [`MOI.Bridges.Variable.HermitianToSymmetricPSDBridge`](@ref).
+
+## Source node
+
+`HermitianToSymmetricPSDBridge` supports:
+
+  * `G` in [`MOI.HermitianPositiveSemidefiniteConeTriangle`](@ref)
+
+## Target node
+
+`HermitianToSymmetricPSDBridge` creates:
+
+  * `F` in [`MOI.PositiveSemidefiniteConeTriangle`](@ref)
+
+## Reformulation
+
+The reformulation is best described by example.
+
+The Hermitian matrix:
+```math
+\\begin{bmatrix}
+  x_{11}            & x_{12} + y_{12}im & x_{13} + y_{13}im\\\\
+  x_{12} - y_{12}im & x_{22}            & x_{23} + y_{23}im\\\\
+  x_{13} - y_{13}im & x_{23} - y_{23}im & x_{33}
+\\end{bmatrix}
+```
+is positive semidefinite if and only if the symmetric matrix:
+```math
+\\begin{bmatrix}
+    x_{11} & x_{12} & x_{13} & 0       & y_{12}  & y_{13} \\\\
+           & x_{22} & x_{23} & -y_{12} & 0       & y_{23} \\\\
+           &        & x_{33} & -y_{13} & -y_{23} & 0      \\\\
+           &        &        & x_{11}  & x_{12}  & x_{13} \\\\
+           &        &        &         & x_{22}  & x_{23} \\\\
+           &        &        &         &         & x_{33}
+\\end{bmatrix}
+```
+is positive semidefinite.
+
+The bridge achieves this reformulation by constraining the above matrix to
+belong to the `MOI.PositiveSemidefiniteConeTriangle(6)`.
+"""
+struct HermitianToSymmetricPSDBridge{T,F,G} <: SetMapBridge{
+    T,
+    MOI.PositiveSemidefiniteConeTriangle,
+    MOI.HermitianPositiveSemidefiniteConeTriangle,
+    F,
+    G,
+}
+    constraint::MOI.ConstraintIndex{F,MOI.PositiveSemidefiniteConeTriangle}
+end
+
+const HermitianToSymmetricPSD{T,OT<:MOI.ModelLike} =
+    SingleBridgeOptimizer{HermitianToSymmetricPSDBridge{T},OT}
+
+function _promote_minus_vcat(::Type{T}, ::Type{G}) where {T,G}
+    S = MOI.Utilities.scalar_type(G)
+    M = MOI.Utilities.promote_operation(-, T, S)
+    F = MOI.Utilities.promote_operation(vcat, T, S, M, T)
+    return F
+end
+
+function concrete_bridge_type(
+    ::Type{<:HermitianToSymmetricPSDBridge{T}},
+    G::Type{<:MOI.AbstractVectorFunction},
+    ::Type{MOI.HermitianPositiveSemidefiniteConeTriangle},
+) where {T}
+    F = _promote_minus_vcat(T, G)
+    return HermitianToSymmetricPSDBridge{T,F,G}
+end
+
+function MOI.Bridges.map_set(
+    ::Type{<:HermitianToSymmetricPSDBridge},
+    set::MOI.HermitianPositiveSemidefiniteConeTriangle,
+)
+    return MOI.PositiveSemidefiniteConeTriangle(2set.side_dimension)
+end
+
+function MOI.Bridges.inverse_map_set(
+    ::Type{<:HermitianToSymmetricPSDBridge},
+    set::MOI.PositiveSemidefiniteConeTriangle,
+)
+    dim = set.side_dimension
+    @assert iseven(dim)
+    return MOI.HermitianPositiveSemidefiniteConeTriangle(div(dim, 2))
+end
+
+function MOI.Bridges.map_function(
+    ::Type{<:HermitianToSymmetricPSDBridge{T}},
+    func,
+) where {T}
+    complex_scalars = MOI.Utilities.eachscalar(func)
+    S = MOI.Utilities.scalar_type(_promote_minus_vcat(T, typeof(func)))
+    complex_dim = length(complex_scalars)
+    complex_set = MOI.Utilities.set_with_dimension(
+        MOI.HermitianPositiveSemidefiniteConeTriangle,
+        complex_dim,
+    )
+    n = complex_set.side_dimension
+    real_set = MOI.PositiveSemidefiniteConeTriangle(2n)
+    real_dim = MOI.dimension(real_set)
+    real_scalars = Vector{S}(undef, real_dim)
+    complex_index = 0
+    half_real_set = MOI.PositiveSemidefiniteConeTriangle(n)
+    half_real_dim = MOI.dimension(half_real_set)
+    real_index_1 = 0
+    real_index_2 = half_real_dim + n
+    for j in 1:n
+        for i in 1:j
+            complex_index += 1
+            real_index_1 += 1
+            real_scalars[real_index_1] = complex_scalars[complex_index]
+            real_index_2 += 1
+            real_scalars[real_index_2] = complex_scalars[complex_index]
+            if i == j
+                real_index_2 += n
+            end
+        end
+    end
+    real_index_1 = half_real_dim
+    real_index_2 = real_dim - n + 1
+    for j in 1:n
+        for i in 1:(j-1)
+            complex_index += 1
+            real_index_1 += 1
+            real_scalars[real_index_1] = complex_scalars[complex_index]
+            real_index_2 -= 1
+            real_scalars[real_index_2] =
+                MOI.Utilities.operate(-, T, complex_scalars[complex_index])
+        end
+        real_scalars[real_index_1+1] = zero(S)
+        real_index_1 += n + 1
+        real_index_2 -= 2 * (n - j) + 1
+    end
+    @assert length(complex_scalars) == complex_index
+    return MOI.Utilities.vectorize(real_scalars)
+end
+
+function MOI.Bridges.inverse_map_function(
+    BT::Type{<:HermitianToSymmetricPSDBridge},
+    func,
+)
+    real_scalars = MOI.Utilities.eachscalar(func)
+    real_set = MOI.Utilities.set_with_dimension(
+        MOI.PositiveSemidefiniteConeTriangle,
+        length(real_scalars),
+    )
+    @assert iseven(real_set.side_dimension)
+    n = div(real_set.side_dimension, 2)
+    complex_set = MOI.HermitianPositiveSemidefiniteConeTriangle(n)
+    complex_scalars =
+        Vector{eltype(real_scalars)}(undef, MOI.dimension(complex_set))
+    real_index = 0
+    complex_index = 0
+    for j in 1:n
+        for i in 1:j
+            complex_index += 1
+            real_index += 1
+            complex_scalars[complex_index] = real_scalars[real_index]
+        end
+    end
+    for j in 1:n
+        for i in 1:(j-1)
+            complex_index += 1
+            real_index += 1
+            complex_scalars[complex_index] = real_scalars[real_index]
+        end
+        real_index += n + 1
+    end
+    @assert length(complex_scalars) == complex_index
+    return MOI.Utilities.vectorize(complex_scalars)
+end
+
+function MOI.Bridges.adjoint_map_function(
+    BT::Type{<:HermitianToSymmetricPSDBridge},
+    func,
+)
+    real_scalars = MOI.Utilities.eachscalar(func)
+    real_dim = length(real_scalars)
+    real_set = MOI.Utilities.set_with_dimension(
+        MOI.PositiveSemidefiniteConeTriangle,
+        real_dim,
+    )
+    @assert iseven(real_set.side_dimension)
+    n = div(real_set.side_dimension, 2)
+    complex_set = MOI.HermitianPositiveSemidefiniteConeTriangle(n)
+    complex_scalars =
+        Vector{eltype(real_scalars)}(undef, MOI.dimension(complex_set))
+    complex_index = 0
+    half_real_set = MOI.PositiveSemidefiniteConeTriangle(n)
+    half_real_dim = MOI.dimension(half_real_set)
+    real_index_1 = 0
+    real_index_2 = half_real_dim + n
+    for j in 1:n
+        for i in 1:j
+            complex_index += 1
+            real_index_1 += 1
+            real_index_2 += 1
+            complex_scalars[complex_index] =
+                real_scalars[real_index_1] + real_scalars[real_index_2]
+            if i == j
+                real_index_2 += n
+            end
+        end
+    end
+    real_index_1 = half_real_dim
+    real_index_2 = real_dim - n + 1
+    for j in 1:n
+        for i in 1:(j-1)
+            complex_index += 1
+            real_index_1 += 1
+            real_index_2 -= 1
+            complex_scalars[complex_index] =
+                real_scalars[real_index_1] - real_scalars[real_index_2]
+        end
+        real_index_1 += n + 1
+        real_index_2 -= 2 * (n - j) + 1
+    end
+    @assert length(complex_scalars) == complex_index
+    return MOI.Utilities.vectorize(complex_scalars)
+end
+
+# FIXME
+# It's not so clear how to do this one since the adjoint is not invertible
+# and it's not obvious how to generate a PSD matrix in the preimage.
+# The following heuristic may not be the best:
+function MOI.Bridges.inverse_adjoint_map_function(
+    BT::Type{<:HermitianToSymmetricPSDBridge{T}},
+    func,
+) where {T}
+    complex_scalars = MOI.Utilities.eachscalar(func)
+    S = MOI.Utilities.scalar_type(_promote_minus_vcat(T, typeof(func)))
+    complex_dim = length(complex_scalars)
+    complex_set = MOI.Utilities.set_with_dimension(
+        MOI.HermitianPositiveSemidefiniteConeTriangle,
+        complex_dim,
+    )
+    n = complex_set.side_dimension
+    real_set = MOI.PositiveSemidefiniteConeTriangle(2n)
+    real_dim = MOI.dimension(real_set)
+    real_scalars = Vector{S}(undef, real_dim)
+    complex_index = 0
+    half_real_set = MOI.PositiveSemidefiniteConeTriangle(n)
+    half_real_dim = MOI.dimension(half_real_set)
+    real_index_1 = 0
+    real_index_2 = half_real_dim + n
+    for j in 1:n
+        for i in 1:j
+            complex_index += 1
+            real_index_1 += 1
+            real_scalars[real_index_1] = complex_scalars[complex_index] / 2
+            real_index_2 += 1
+            real_scalars[real_index_2] = complex_scalars[complex_index] / 2
+            if i == j
+                real_index_2 += n
+            end
+        end
+    end
+    real_index_1 = half_real_dim
+    real_index_2 = real_dim - n + 1
+    for j in 1:n
+        for i in 1:(j-1)
+            complex_index += 1
+            real_index_1 += 1
+            real_scalars[real_index_1] = complex_scalars[complex_index]
+            real_index_2 -= 1
+            real_scalars[real_index_2] = zero(S)
+        end
+        real_scalars[real_index_1+1] = zero(S)
+        real_index_1 += n + 1
+        real_index_2 -= 2 * (n - j) + 1
+    end
+    @assert length(complex_scalars) == complex_index
+    return MOI.Utilities.vectorize(real_scalars)
+end
diff --git a/src/Utilities/functions.jl b/src/Utilities/functions.jl
@@ -403,6 +403,10 @@ end
 Type of functions obtained by indexing objects obtained by calling `eachscalar`
 on functions of type `F`.
 """
+function scalar_type end
+
+scalar_type(::Type{<:AbstractVector{T}}) where {T} = T
+
 scalar_type(::Type{MOI.VectorOfVariables}) = MOI.VariableIndex
 
 function scalar_type(::Type{MOI.VectorAffineFunction{T}}) where {T}
@@ -1527,10 +1531,14 @@ const TypedLike{T} = Union{TypedScalarLike{T},TypedVectorLike{T}}
 ###################################### +/- #####################################
 ## promote_operation
 
+function promote_operation(::typeof(-), ::Type{T}, ::Type{T}) where {T}
+    return T
+end
+
 function promote_operation(
     ::typeof(-),
     ::Type{T},
-    ::Type{<:ScalarAffineLike{T}},
+    ::Type{<:Union{MOI.VariableIndex,MOI.ScalarAffineFunction{T}}},
 ) where {T}
     return MOI.ScalarAffineFunction{T}
 end
@@ -2895,6 +2903,13 @@ function fill_constant(
     return
 end
 
+"""
+    vectorize(x::AbstractVector{<:Number})
+
+Returns `x`.
+"""
+vectorize(x::AbstractVector{<:Number}) = x
+
 """
     vectorize(x::AbstractVector{MOI.VariableIndex})
 
diff --git a/src/Utilities/matrix_of_constraints.jl b/src/Utilities/matrix_of_constraints.jl
@@ -576,6 +576,17 @@ function set_with_dimension(
     return S(isqrt(dim))
 end
 
+function set_with_dimension(
+    ::Type{MOI.HermitianPositiveSemidefiniteConeTriangle},
+    dim,
+)
+    # We have `n*(n+1)/2 + n*(n-1)/2 = dim` so
+    # `n² + n + n² - n = 2dim` hence `n² = dim`
+    # This can be seen geometrically as the vectorization
+    # contains the upper triangular followed by the strictly upper triangular.
+    return MOI.HermitianPositiveSemidefiniteConeTriangle(isqrt(dim))
+end
+
 function set_with_dimension(::Type{MOI.LogDetConeTriangle}, dim)
     side_dimension = side_dimension_for_vectorized_dimension(dim - 2)
     return MOI.LogDetConeTriangle(side_dimension)
diff --git a/src/Utilities/sets.jl b/src/Utilities/sets.jl
@@ -130,6 +130,10 @@ Return the dimension `d` such that
 `MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(d))` is `n`.
 """
 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.
+    # We can drop the `- 1` as `div` already discards it.
     return div(isqrt(1 + 8n), 2)
 end
 
diff --git a/test/Bridges/Constraint/hermitian.jl b/test/Bridges/Constraint/hermitian.jl
