diff --git a/src/sd.jl b/src/sd.jl index 5f9e986d312..0dc684affa8 100644 --- a/src/sd.jl +++ b/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 " * diff --git a/src/variables.jl b/src/variables.jl index 64189327e34..329e3573bd8 100644 --- a/src/variables.jl +++ b/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) diff --git a/test/test_constraint.jl b/test/test_constraint.jl index c3880002d4f..d2b421e9a97 100644 --- a/test/test_constraint.jl +++ b/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)