Add HermitianPositiveSemidefiniteConeTriangle with bridges (#1962)

* Add bridge from hermitian PSD to PSD * Update src/Bridges/Variable/bridges/hermitian.jl Co-authored-by: Oscar Dowson <[email protected]> * Update src/Bridges/Variable/bridges/hermitian.jl Co-authored-by: Oscar Dowson <[email protected]> * Update src/Test/test_conic.jl Co-authored-by: Oscar Dowson <[email protected]> * Update src/Bridges/Variable/bridges/hermitian.jl Co-authored-by: Oscar Dowson <[email protected]> * Update src/Bridges/Variable/bridges/hermitian.jl Co-authored-by: Oscar Dowson <[email protected]> * Update src/Test/test_conic.jl Co-authored-by: Oscar Dowson <[email protected]> * Update src/Test/test_conic.jl Co-authored-by: Oscar Dowson <[email protected]> * Fix format * Address review comments * Add docstring * Add bridge tests * Tidy new tests * Use import LinearAlgebra in tests * Simplifications * Update coverage * Add explanation back * Clean up test * Fix * Update src/Test/test_conic.jl Co-authored-by: Oscar Dowson <[email protected]> Co-authored-by: odow <[email protected]>
jump-dev · Jul 26, 2022 · 8c8636f · 8c8636f
1 parent db2c00f
commit 8c8636f
Show file tree

Hide file tree

Showing 10 changed files with 731 additions and 1 deletion.
diff --git a/docs/src/manual/standard_form.md b/docs/src/manual/standard_form.md
@@ -93,7 +93,9 @@ The matrix-valued set types implemented in MathOptInterface.jl are:
 | [`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} \}``
+| [`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. |
 
 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
@@ -141,6 +141,7 @@ List of recognized matrix sets.
 ```@docs
 PositiveSemidefiniteConeTriangle
 PositiveSemidefiniteConeSquare
+HermitianPositiveSemidefiniteConeTriangle
 LogDetConeTriangle
 LogDetConeSquare
 RootDetConeTriangle

diff --git a/docs/src/submodules/Bridges/list_of_bridges.md b/docs/src/submodules/Bridges/list_of_bridges.md
@@ -84,4 +84,5 @@ Bridges.Variable.RSOCtoSOCBridge
 Bridges.Variable.SOCtoRSOCBridge
 Bridges.Variable.VectorizeBridge
 Bridges.Variable.ZerosBridge
+Bridges.Variable.HermitianToSymmetricPSDBridge
 ```
diff --git a/src/Bridges/Variable/Variable.jl b/src/Bridges/Variable/Variable.jl
@@ -22,6 +22,7 @@ include("bridges/rsoc_soc.jl")
 include("bridges/soc_rsoc.jl")
 include("bridges/vectorize.jl")
 include("bridges/zeros.jl")
+include("bridges/hermitian.jl")
 
 """
     add_all_bridges(model, ::Type{T}) where {T}
@@ -38,6 +39,7 @@ function add_all_bridges(model, ::Type{T}) where {T}
     MOI.Bridges.add_bridge(model, SOCtoRSOCBridge{T})
     MOI.Bridges.add_bridge(model, RSOCtoSOCBridge{T})
     MOI.Bridges.add_bridge(model, RSOCtoPSDBridge{T})
+    MOI.Bridges.add_bridge(model, HermitianToSymmetricPSDBridge{T})
     return
 end
 

diff --git a/src/Bridges/Variable/bridges/hermitian.jl b/src/Bridges/Variable/bridges/hermitian.jl
@@ -0,0 +1,334 @@
+# 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} <: Bridges.Variable.AbstractBridge
+
+`HermitianToSymmetricPSDBridge` implements the following reformulation:
+
+  * Hermitian positive semidefinite `n x n` complex matrix to a symmetric
+    positive semidefinite `2n x 2n` real matrix satisfying equality constraints
+    described below.
+
+## Source node
+
+`HermitianToSymmetricPSDBridge` supports:
+
+  * [`MOI.VectorOfVariables`](@ref) in
+    [`MOI.HermitianPositiveSemidefiniteConeTriangle`](@ref)
+
+## Target node
+
+`HermitianToSymmetricPSDBridge` creates:
+
+  * [`MOI.VectorOfVariables`](@ref) in [`MOI.PositiveSemidefiniteConeTriangle`](@ref)
+  * [`MOI.ScalarAffineFunction{T}`](@ref) in [`MOI.EqualTo{T}`](@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 adding a new set of variables in
+`MOI.PositiveSemidefiniteConeTriangle(6)`, and then adding three groups of
+equality constraints to:
+
+ * constrain the two `x` blocks to be equal
+ * force the diagonal of the `y` blocks to be `0`
+ * force the lower triangular of the `y` block to be the negative of the upper
+   triangle.
+"""
+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}}}
+end
+
+const HermitianToSymmetricPSD{T,OT<:MOI.ModelLike} =
+    SingleBridgeOptimizer{HermitianToSymmetricPSDBridge{T},OT}
+
+function bridge_constrained_variable(
+    ::Type{HermitianToSymmetricPSDBridge{T}},
+    model::MOI.ModelLike,
+    set::MOI.HermitianPositiveSemidefiniteConeTriangle,
+) where {T}
+    n = set.side_dimension
+    variables, psd_ci = MOI.add_constrained_variables(
+        model,
+        MOI.PositiveSemidefiniteConeTriangle(2n),
+    )
+    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
+            k11 += 1
+            k12 += 1
+            k22 += 1
+            f_x = MOI.Utilities.operate(-, T, variables[k11], variables[k22])
+            push!(ceq, MOI.add_constraint(model, f_x, MOI.EqualTo(zero(T))))
+            if i == j  # y_{ii} = 0
+                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])
+                push!(ceq, MOI.add_constraint(model, f_y, MOI.EqualTo(zero(T))))
+            end
+        end
+        k12 += n
+    end
+    return HermitianToSymmetricPSDBridge(variables, psd_ci, n, ceq)
+end
+
+function supports_constrained_variable(
+    ::Type{<:HermitianToSymmetricPSDBridge},
+    ::Type{MOI.HermitianPositiveSemidefiniteConeTriangle},
+)
+    return true
+end
+
+function MOI.Bridges.added_constrained_variable_types(
+    ::Type{<:HermitianToSymmetricPSDBridge},
+)
+    return Tuple{Type}[(MOI.PositiveSemidefiniteConeTriangle,)]
+end
+
+function MOI.Bridges.added_constraint_types(
+    ::Type{HermitianToSymmetricPSDBridge{T}},
+) where {T}
+    return Tuple{Type,Type}[(MOI.ScalarAffineFunction{T}, MOI.EqualTo{T})]
+end
+
+function MOI.get(bridge::HermitianToSymmetricPSDBridge, ::MOI.NumberOfVariables)
+    return length(bridge.variables)
+end
+
+function MOI.get(
+    bridge::HermitianToSymmetricPSDBridge,
+    ::MOI.ListOfVariableIndices,
+)
+    return copy(bridge.variables)
+end
+
+function MOI.get(
+    ::HermitianToSymmetricPSDBridge,
+    ::MOI.NumberOfConstraints{
+        MOI.VectorOfVariables,
+        MOI.PositiveSemidefiniteConeTriangle,
+    },
+)::Int64
+    return 1
+end
+
+function MOI.get(
+    bridge::HermitianToSymmetricPSDBridge,
+    ::MOI.ListOfConstraintIndices{
+        MOI.VectorOfVariables,
+        MOI.PositiveSemidefiniteConeTriangle,
+    },
+)
+    return [bridge.psd]
+end
+
+function MOI.get(
+    bridge::HermitianToSymmetricPSDBridge{T},
+    ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}},
+) where {T}
+    return length(bridge.ceq)
+end
+
+function MOI.get(
+    bridge::HermitianToSymmetricPSDBridge{T},
+    ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}},
+) where {T}
+    return copy(bridge.ceq)
+end
+
+function MOI.delete(model::MOI.ModelLike, bridge::HermitianToSymmetricPSDBridge)
+    MOI.delete(model, bridge.ceq)
+    MOI.delete(model, bridge.variables)
+    return
+end
+
+function MOI.get(
+    ::MOI.ModelLike,
+    ::MOI.ConstraintSet,
+    bridge::HermitianToSymmetricPSDBridge,
+)
+    return MOI.HermitianPositiveSemidefiniteConeTriangle(bridge.n)
+end
+
+function _matrix_indices(k)
+    # If `k` is a diagonal index, `s(k)` is odd and 1 + 8k is a perfect square.
+    n = 1 + 8k
+    s = isqrt(n)
+    j = if s^2 == n
+        div(s, 2)
+    else
+        # Otherwise, if it is after the diagonal index `k` but before the
+        # diagonal index `k'` with `s(k') = s(k) + 2`, we have
+        # `s(k) <= s < s(k) + 2`.
+        # By shifting by `+1` before `div`, we make sure to have the right
+        # column.
+        div(s + 1, 2)
+    end
+    i = k - MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(j - 1))
+    return i, j
+end
+
+function _variable_map(idx::MOI.Bridges.IndexInVector, n)
+    N = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n))
+    if idx.value <= N
+        return idx.value
+    end
+    i, j = _matrix_indices(idx.value - N)
+    d = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(j))
+    return N + j * n + d + i
+end
+
+function _variable(
+    bridge::HermitianToSymmetricPSDBridge,
+    i::MOI.Bridges.IndexInVector,
+)
+    return bridge.variables[_variable_map(i, bridge.n)]
+end
+
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::MOI.ConstraintPrimal,
+    bridge::HermitianToSymmetricPSDBridge{T},
+) where {T}
+    values = MOI.get(model, attr, bridge.psd)
+    M = MOI.dimension(MOI.get(model, MOI.ConstraintSet(), bridge))
+    n = bridge.n
+    return [values[_variable_map(MOI.Bridges.IndexInVector(i), n)] for i in 1:M]
+end
+
+# We don't need to take into account the equality constraints. We just need to
+# sum up (with appropriate +/-) each dual variable associated with the original
+# x or y element.
+# The reason for this is as follows:
+# Suppose for simplicity that the elements of a `2n x 2n` matrix are ordered as:
+# ```
+# \\ 1 |\\ 2
+#  \\  | 3
+#   \\ |4_\\
+#    \\  5
+#     \\
+#      \\
+# ```
+# Let `H = HermitianToSymmetricPSDBridge(n)`,
+# `S = PositiveSemidefiniteConeTriangle(2n)` and `ceq` be the linear space of
+# `2n x 2n` symmetric matrices such that the block `1` and `5` are equal, `2` and `4` are opposite and `3` is zero.
+# We consider the cone `P = S ∩ ceq`.
+# We have `P = A * H` where
+# ```
+#     [I  0]
+#     [0  I]
+# A = [0  0]
+#     [0 -I]
+#     [I  0]
+# ```
+# Therefore, `H* = A* * P*` where
+# ```
+#      [I 0 0  0 I]
+# A* = [0 I 0 -I 0]
+# ```
+# Moreover, as `(S ∩ T)* = S* + T*` for cones `S` and `T`, we have
+# ```
+# P* = S* + ceq*
+# ```
+# the dual vector of `P*` is the dual vector of `S*` for which we add in the corresponding
+# entries the dual of the three constraints, multiplied by the coefficients for the `EqualTo` constraints.
+# Note that these contributions cancel out when we multiply them by `A*`:
+# A* * (S* + ceq*) = A* * S*
+# so we can just ignore them.
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::MOI.ConstraintDual,
+    bridge::HermitianToSymmetricPSDBridge{T},
+) where {T}
+    dual = MOI.get(model, attr, bridge.psd)
+    M = MOI.dimension(MOI.get(model, MOI.ConstraintSet(), bridge))
+    result = zeros(T, M)
+    n = bridge.n
+    N = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n))
+    k11, k12, k22 = 0, N, N
+    k21 = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(2n)) + 1
+    k = 0
+    for j in 1:n
+        k21 -= n + 1 - j
+        k22 += n
+        for i in 1:j
+            k11 += 1
+            k12 += 1
+            k21 -= 1
+            k22 += 1
+            result[k11] += dual[k11] + dual[k22]
+            if i != j
+                k += 1
+                result[N+k] += dual[k12] - dual[k21]
+            end
+        end
+        k12 += n
+        k21 -= n - j
+    end
+    return result
+end
+
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::MOI.VariablePrimal,
+    bridge::HermitianToSymmetricPSDBridge{T},
+    i::MOI.Bridges.IndexInVector,
+) where {T}
+    return MOI.get(model, attr, _variable(bridge, i))
+end
+
+function MOI.Bridges.bridged_function(
+    bridge::HermitianToSymmetricPSDBridge{T},
+    i::MOI.Bridges.IndexInVector,
+) where {T}
+    return convert(MOI.ScalarAffineFunction{T}, _variable(bridge, i))
+end
+
+function unbridged_map(
+    bridge::HermitianToSymmetricPSDBridge{T},
+    vi::MOI.VariableIndex,
+    i::MOI.Bridges.IndexInVector,
+) where {T}
+    return (_variable(bridge, i) => convert(MOI.ScalarAffineFunction{T}, vi),)
+end
diff --git a/src/Test/Test.jl b/src/Test/Test.jl
@@ -6,6 +6,8 @@
 
 module Test
 
+import LinearAlgebra
+
 using MathOptInterface
 const MOI = MathOptInterface
 const MOIU = MOI.Utilities