Use adjoint map for Hermitian PSD variable bridge

src/Bridges/Constraint/bridges/HermitianToSymmetricPSDBridge.jl

@@ -188,32 +188,36 @@ function MOI.Bridges.inverse_map_function(
 end
 
 function MOI.Bridges.adjoint_map_function(
-    BT::Type{<:HermitianToSymmetricPSDBridge},
+    ::Type{<:HermitianToSymmetricPSDBridge{T}},
     func,
-)
-    real_scalars = MOI.Utilities.eachscalar(func)
-    real_dim = length(real_scalars)
+) where {T}
+    real_dim = MOI.Utilities.output_dim(T, func)
     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))
+    S = MOI.Utilities.scalar_type(_promote_minus_vcat(T, typeof(func)))
+    scalars = Vector{S}(undef, n^2)
+    real_scalars = MOI.Utilities.eachscalar(func)
+    # Traverse the full triangle sequentially; inverse index lookups are only
+    # needed when requesting an individual component.
     complex_index = 0
-    half_real_set = MOI.PositiveSemidefiniteConeTriangle(n)
-    half_real_dim = MOI.dimension(half_real_set)
+    half_real_dim = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n))
     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]
+            scalars[complex_index] = MOI.Utilities.operate(
+                +,
+                T,
+                real_scalars[real_index_1],
+                real_scalars[real_index_2],
+            )
             if i == j
                 real_index_2 += n
             end
@@ -224,63 +228,46 @@ function MOI.Bridges.adjoint_map_function(
         for i in 1:(j-1)
             complex_index += 1
             real_index_1 += 1
-            complex_scalars[complex_index] =
-                real_scalars[real_index_1] -
-                real_scalars[MOI.Utilities.trimap(j, i + n)]
+            scalars[complex_index] = MOI.Utilities.operate(
+                -,
+                T,
+                real_scalars[real_index_1],
+                real_scalars[MOI.Utilities.trimap(j, i + n)],
+            )
         end
         real_index_1 += n + 1
     end
-    @assert length(complex_scalars) == complex_index
-    return MOI.Utilities.vectorize(complex_scalars)
+    return MOI.Utilities.vectorize(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}},
+function MOI.Bridges.adjoint_map_function(
+    ::Type{<:HermitianToSymmetricPSDBridge{T}},
     func,
+    index::MOI.Bridges.IndexInVector,
 ) 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,
+    scalars = MOI.Utilities.eachscalar(func)
+    n = div(
+        MOI.Utilities.side_dimension_for_vectorized_dimension(length(scalars)),
+        2,
     )
-    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
+    N = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n))
+    if index.value <= N
+        i, j = MOI.Utilities.inverse_trimap(index.value)
+        k = MOI.Utilities.trimap(n + i, n + j)
+        return MOI.Utilities.operate(+, T, scalars[index.value], scalars[k])
     end
-    real_index_1 = half_real_dim
-    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_scalars[MOI.Utilities.trimap(j, i + n)] = zero(S)
-        end
-        real_scalars[real_index_1+1] = zero(S)
-        real_index_1 += n + 1
-    end
-    @assert length(complex_scalars) == complex_index
-    return MOI.Utilities.vectorize(real_scalars)
+    i, j = MOI.Utilities.inverse_trimap(index.value - N)
+    j += 1  # The imaginary part excludes the diagonal.
+    k12 = MOI.Utilities.trimap(i, n + j)
+    k21 = MOI.Utilities.trimap(j, n + i)
+    return MOI.Utilities.operate(-, T, scalars[k12], scalars[k21])
+end
+
+function MOI.Bridges.inverse_adjoint_map_function(
+    BT::Type{<:HermitianToSymmetricPSDBridge{T}},
+    func,
+) where {T}
+    # F* F = 2I, so F / 2 is a right inverse of F* that preserves PSD.
+    mapped = MOI.Bridges.map_function(BT, func)
+    return MOI.Utilities.operate(/, T, mapped, T(2))
 end