Update

src/sd.jl

@@ -141,15 +141,50 @@ MathOptInterface.PositiveSemidefiniteConeTriangle(2)
 struct PSDCone end
 
 """
-    SymmetricMatrixShape
+    SymmetricMatrixShape(
+        side_dimension::Int;
+        needs_adjoint_dual::Bool = false,
+    )
+
+The shape object for a symmetric square matrix of `side_dimension` rows and
+columns.
 
-Shape object for a symmetric square matrix of `side_dimension` rows and columns.
 The vectorized form contains the entries of the upper-right triangular part of
 the matrix given column by column (or equivalently, the entries of the
 lower-left triangular part given row by row).
+
+## `needs_adjoint_dual`
+
+By default, the [`dual_shape`](@ref) of [`SymmetricMatrixShape`](@ref) is also
+[`SymmetricMatrixShape`](@ref). This is true for cases such as a
+`LinearAlgebra.Symmetric` matrix in [`PSDCone`](@ref).
+
+However, JuMP also supports `LinearAlgebra.Symmetric` matrix in `Zeros`, which
+is interpreted as an element-wise equality constraint. By exploiting symmetry,
+we pass only the upper triangle of the equality constraints. This works for the
+primal, but it leads to a factor of 2 difference in the off-diagonal dual
+elements. (The dual value of the `(i, j)` element in the triangle formulation
+should be divided by 2 when spread across the `(i, j)` and `(j, i)` elements in
+the square matrix formulation.) If the constraint has this dual inconsistency,
+set `needs_adjoint_dual = true`.
 """
 struct SymmetricMatrixShape <: AbstractShape
     side_dimension::Int
+    needs_adjoint_dual::Bool
+
+    function SymmetricMatrixShape(
+        side_dimension::Int;
+        needs_adjoint_dual::Bool = false,
+    )
+        return new(side_dimension, needs_adjoint_dual)
+    end
+end
+
+function dual_shape(s::SymmetricMatrixShape)
+    if s.needs_adjoint_dual
+        return SymmetricMatrixAdjointShape(s.side_dimension)
+    end
+    return s
 end
 
 function reshape_vector(
@@ -175,15 +210,25 @@ function reshape_set(
 end
 
 """
-    DualSymmetricMatrixShape(side_dimension)
+    SymmetricMatrixAdjointShape(side_dimension)
 
-A shape used when quering the dual of a symmetic matrix.
+The [`dual_shape`](@ref) of [`SymmetricMatrixShape`](@ref).
+
+This shape is not intended for regular use.
 """
-struct DualSymmetricMatrixShape <: AbstractShape
+struct SymmetricMatrixAdjointShape <: AbstractShape
     side_dimension::Int
 end
 
-function reshape_vector(v::Vector{T}, shape::DualSymmetricMatrixShape) where {T}
+function vectorize(matrix::AbstractMatrix, s::SymmetricMatrixAdjointShape)
+    n = LinearAlgebra.checksquare(matrix)
+    return [((i == j) ? 1 : 2) * matrix[i, j] for j in 1:n for i in 1:j]
+end
+
+function reshape_vector(
+    v::Vector{T},
+    shape::SymmetricMatrixAdjointShape,
+) where {T}
     matrix = Matrix{T}(undef, shape.side_dimension, shape.side_dimension)
     k = 0
     for j in 1:shape.side_dimension
@@ -455,14 +500,49 @@ parametrize a Hermitian matrix.
 struct HermitianPSDCone end
 
 """
-    HermitianMatrixShape
+    HermitianMatrixShape(
+        side_dimension::Int;
+        needs_adjoint_dual::Bool = false,
+    )
+
+The shape object for a Hermitian square matrix of `side_dimension` rows and
+columns.
 
-Shape object for a Hermitian square matrix of `side_dimension` rows and
-columns. The vectorized form corresponds to
+The vectorized form corresponds to
 [`MOI.HermitianPositiveSemidefiniteConeTriangle`](@ref).
+
+## `needs_adjoint_dual`
+
+By default, the [`dual_shape`](@ref) of [`HermitianMatrixShape`](@ref) is also
+[`HermitianMatrixShape`](@ref). This is true for cases such as a
+`LinearAlgebra.Hermitian` matrix in [`HermitianPSDCone`](@ref).
+
+However, JuMP also supports `LinearAlgebra.Hermitian` matrix in `Zeros`, which
+is interpreted as an element-wise equality constraint. By exploiting symmetry,
+we pass only the upper triangle of the equality constraints. This works for the
+primal, but it leads to a factor of 2 difference in the off-diagonal dual
+elements. (The dual value of the `(i, j)` element in the triangle formulation
+should be divided by 2 when spread across the `(i, j)` and `(j, i)` elements in
+the square matrix formulation.) If the constraint has this dual inconsistency,
+set `needs_adjoint_dual = true`.
 """
 struct HermitianMatrixShape <: AbstractShape
     side_dimension::Int
+    needs_adjoint_dual::Bool
+
+    function HermitianMatrixShape(
+        side_dimension::Int;
+        needs_adjoint_dual::Bool = false,
+    )
+        return new(side_dimension, needs_adjoint_dual)
+    end
+end
+
+function dual_shape(s::HermitianMatrixShape)
+    if s.needs_adjoint_dual
+        return HermitianMatrixAdjointShape(s.side_dimension)
+    end
+    return s
 end
 
 function vectorize(matrix, ::HermitianMatrixShape)
@@ -503,15 +583,30 @@ function reshape_vector(v::Vector{T}, shape::HermitianMatrixShape) where {T}
 end
 
 """
-    DualHermitianMatrixShape(side_dimension)
+    HermitianMatrixAdjointShape(side_dimension)
 
-A shape used when quering the dual of a symmetic matrix.
+The [`dual_shape`](@ref) of [`HermitianMatrixShape`](@ref).
+
+This shape is not intended for regular use.
 """
-struct DualHermitianMatrixShape <: AbstractShape
+struct HermitianMatrixAdjointShape <: AbstractShape
     side_dimension::Int
 end
 
-function reshape_vector(v::Vector{T}, shape::DualHermitianMatrixShape) where {T}
+function vectorize(matrix, ::HermitianMatrixAdjointShape)
+    n = LinearAlgebra.checksquare(matrix)
+    real_shape = SymmetricMatrixAdjointShape(n)
+    imag_shape = SymmetricMatrixShape(n - 1)
+    return vcat(
+        vectorize(_real.(matrix), real_shape),
+        vectorize(2 * _imag.(matrix[1:(end-1), 2:end]), imag_shape),
+    )
+end
+
+function reshape_vector(
+    v::Vector{T},
+    shape::HermitianMatrixAdjointShape,
+) where {T}
     NewType = _MA.promote_operation(_MA.add_mul, T, Complex{Bool}, T)
     n = shape.side_dimension
     matrix = Matrix{NewType}(undef, n, n)
@@ -603,25 +698,11 @@ function build_constraint(
     ::Zeros,
 )
     n = LinearAlgebra.checksquare(H)
-    shape = HermitianMatrixShape(n)
+    shape = HermitianMatrixShape(n; needs_adjoint_dual = true)
     x = vectorize(H, shape)
     return VectorConstraint(x, MOI.Zeros(length(x)), shape)
 end
 
-function dual(
-    con_ref::ConstraintRef{
-        <:AbstractModel,
-        MOI.ConstraintIndex{F,MOI.Zeros},
-        HermitianMatrixShape,
-    };
-    result::Int = 1,
-) where {F}
-    return reshape_vector(
-        _constraint_dual(con_ref, result),
-        DualHermitianMatrixShape(con_ref.shape.side_dimension),
-    )
-end
-
 reshape_set(s::MOI.Zeros, ::HermitianMatrixShape) = Zeros()
 
 function build_constraint(
@@ -630,25 +711,11 @@ function build_constraint(
     ::Zeros,
 )
     n = LinearAlgebra.checksquare(f)
-    shape = SymmetricMatrixShape(n)
+    shape = SymmetricMatrixShape(n; needs_adjoint_dual = true)
     x = vectorize(f, shape)
     return VectorConstraint(x, MOI.Zeros(length(x)), shape)
 end
 
-function dual(
-    con_ref::ConstraintRef{
-        <:AbstractModel,
-        MOI.ConstraintIndex{F,MOI.Zeros},
-        SymmetricMatrixShape,
-    };
-    result::Int = 1,
-) where {F}
-    return reshape_vector(
-        _constraint_dual(con_ref, result),
-        DualSymmetricMatrixShape(con_ref.shape.side_dimension),
-    )
-end
-
 function build_constraint(error_fn::Function, ::AbstractMatrix, ::Nonnegatives)
     return error_fn(
         "Unsupported matrix in vector-valued set. Did you mean to use the " *

src/variables.jl

@@ -2304,10 +2304,10 @@ function _imag(s::String)
 end
 
 _real(v::ScalarVariable) = _mapinfo(real, v)
-_real(scalar::AbstractJuMPScalar) = real(scalar)
+_real(scalar) = real(scalar)
 
 _imag(v::ScalarVariable) = _mapinfo(imag, v)
-_imag(scalar::AbstractJuMPScalar) = imag(scalar)
+_imag(scalar) = imag(scalar)
 
 _conj(v::ScalarVariable) = _mapinfo(conj, v)
 

test/test_constraint.jl

@@ -1635,6 +1635,12 @@ function test_hermitian_in_zeros()
     model = Model()
     @variable(model, H[1:2, 1:2] in HermitianPSDCone())
     c = @constraint(model, H == LinearAlgebra.I)
+    @test c.shape == HermitianMatrixShape(2; needs_adjoint_dual = true)
+    set_dual_start_value(c, [1 2+1im; 2-1im 3])
+    @test dual_start_value(c) ==
+          LinearAlgebra.Hermitian([1.0 2.0+1.0im; 2.0-1.0im 3.0])
+    @test MOI.get(backend(model), MOI.ConstraintDualStart(), index(c)) ==
+          [1.0, 4.0, 3.0, 2.0]
     con = constraint_object(c)
     @test isequal_canonical(con.func[1], real(H[1, 1]) - 1)
     @test isequal_canonical(con.func[2], real(H[1, 2]))
@@ -1650,6 +1656,11 @@ function test_symmetric_in_zeros()
     model = Model()
     @variable(model, H[1:2, 1:2], Symmetric)
     c = @constraint(model, H == LinearAlgebra.I)
+    @test c.shape == SymmetricMatrixShape(2; needs_adjoint_dual = true)
+    set_dual_start_value(c, [1 2; 2 3])
+    @test dual_start_value(c) == LinearAlgebra.Symmetric([1.0 2.0; 2.0 3.0])
+    @test MOI.get(backend(model), MOI.ConstraintDualStart(), index(c)) ==
+          [1.0, 4.0, 3.0]
     con = constraint_object(c)
     @test isequal_canonical(con.func[1], H[1, 1] - 1)
     @test isequal_canonical(con.func[2], H[1, 2] - 0)