From 95da051b3be441999c62057fa4834a3bd8882329 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Tue, 6 Aug 2024 09:01:38 +1200 Subject: [PATCH] Fix querying dual of Symmetric and Hermitian equality constraints (#3797) --- src/sd.jl | 160 ++++++++++++++++++++++++++++++++++++++-- src/variables.jl | 4 +- test/test_constraint.jl | 11 +++ 3 files changed, 166 insertions(+), 9 deletions(-) diff --git a/src/sd.jl b/src/sd.jl index fdb370a5722..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( @@ -174,6 +209,39 @@ function reshape_set( return PSDCone() end +""" + SymmetricMatrixAdjointShape(side_dimension) + +The [`dual_shape`](@ref) of [`SymmetricMatrixShape`](@ref). + +This shape is not intended for regular use. +""" +struct SymmetricMatrixAdjointShape <: AbstractShape + side_dimension::Int +end + +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 + for i in 1:j-1 + k += 1 + matrix[j, i] = matrix[i, j] = 0.5 * v[k] + end + k += 1 + matrix[j, j] = v[k] + end + return LinearAlgebra.Symmetric(matrix) +end + """ triangle_vec(matrix::Matrix) @@ -432,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) @@ -479,6 +582,49 @@ function reshape_vector(v::Vector{T}, shape::HermitianMatrixShape) where {T} return LinearAlgebra.Hermitian(matrix) end +""" + HermitianMatrixAdjointShape(side_dimension) + +The [`dual_shape`](@ref) of [`HermitianMatrixShape`](@ref). + +This shape is not intended for regular use. +""" +struct HermitianMatrixAdjointShape <: AbstractShape + side_dimension::Int +end + +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) + real_k = 0 + imag_k = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n)) + for j in 1:n + for i in 1:(j-1) + real_k += 1 + imag_k += 1 + matrix[i, j] = (v[real_k] + im * v[imag_k]) / 2 + matrix[j, i] = (v[real_k] - im * v[imag_k]) / 2 + end + real_k += 1 + matrix[j, j] = v[real_k] + end + return LinearAlgebra.Hermitian(matrix) +end + function _vectorize_complex_variables(error_fn::Function, matrix::Matrix) if any(_is_binary, matrix) || any(_is_integer, matrix) # We would then need to fix the imaginary value to zero. Let's wait to @@ -552,7 +698,7 @@ 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 @@ -565,7 +711,7 @@ 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 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)