Constraining Symmetric matrix to be equal to a Hermitian matrix is not Hermitian

[araujoms]

Bug introduced in #3797: before if I tried to constrain A == B for Symmetric A and Hermitian B it would give me an error. Now it silently fails.

using JuMP
using LinearAlgebra
import Hypatia
import Ket
import Random

function dual_test(d)
    Random.seed!(1)
    ρ = Ket.random_state(Float64, d^2, 1)

    model = Model()
    @variable(model, σ[1:d^2, 1:d^2] in PSDCone())
    @variable(model, λ)
    noisy_state = Hermitian(ρ + λ * I(d^2))
    @constraint(model, witness_constraint, σ == noisy_state)
    @objective(model, Min, λ)
    @constraint(model, Ket.partial_transpose(σ, 2, [d, d]) in PSDCone())

    set_optimizer(model, Hypatia.Optimizer)
    optimize!(model)

    W = dual(witness_constraint)
    return W
end

[odow]

This is expected behavior from #3778

Previously:

julia> using JuMP, LinearAlgebra

julia> model = Model()
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.

julia> @variable(model, A[1:2, 1:2] in PSDCone())
2×2 Symmetric{VariableRef, Matrix{VariableRef}}:
 A[1,1]  A[1,2]
 A[1,2]  A[2,2]

julia> B = LinearAlgebra.Hermitian(rand(ComplexF64, 2, 2))
2×2 Hermitian{ComplexF64, Matrix{ComplexF64}}:
 0.0952809+0.0im       0.0199259+0.447844im
 0.0199259-0.447844im  0.0929153+0.0im

julia> @constraint(model, A == B)
ERROR: At REPL[7]:1: `@constraint(model, A == B)`: Unsupported matrix in vector-valued set. Did you mean to use the broadcasting syntax `.==` for element-wise equality? Alternatively, this syntax is supported in the special case that the matrices are `LinearAlgebra.Symmetric` or `LinearAlgebra.Hermitian`.
Stacktrace:
 [1] error(::String, ::String)

Now

julia> using JuMP, LinearAlgebra

julia> model = Model()
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.

julia> @variable(model, A[1:2, 1:2] in PSDCone())
2×2 Symmetric{VariableRef, Matrix{VariableRef}}:
 A[1,1]  A[1,2]
 A[1,2]  A[2,2]

julia> B = LinearAlgebra.Hermitian(rand(ComplexF64, 2, 2))
2×2 Hermitian{ComplexF64, Matrix{ComplexF64}}:
 0.317083+0.0im       0.432169+0.203233im
 0.432169-0.203233im  0.198036+0.0im

julia> @constraint(model, A == B)
[A[1,1] - 0.31708259602015787                            …  A[1,2] + (-0.43216931547981996 - 0.2032326360515233im)
 A[1,2] + (-0.43216931547981996 + 0.2032326360515233im)     A[2,2] - 0.19803569737893978] ∈ Zeros()

julia> print(backend(model))
Feasibility

Subject to:

VectorOfVariables-in-PositiveSemidefiniteConeTriangle
 ┌      ┐
 │A[1,1]│
 │A[1,2]│
 │A[2,2]│
 └      ┘ ∈ PositiveSemidefiniteConeTriangle(2)

VectorAffineFunction{ComplexF64}-in-Zeros
 ┌                                                                    ┐
 │(-0.31708259602015787 - 0.0im) + (1.0 + 0.0im) A[1,1]               │
 │(-0.43216931547981996 + 0.2032326360515233im) + (1.0 + 0.0im) A[1,2]│
 │(-0.43216931547981996 - 0.2032326360515233im) + (1.0 + 0.0im) A[1,2]│
 │(-0.19803569737893978 - 0.0im) + (1.0 + 0.0im) A[2,2]               │
 └                                                                    ┘ ∈ Zeros(4)

The constraint A == B seems like a perfectly reasonable constraint to write.

[araujoms]

In your example you are constraining A[1,2] to be simultaneously equal to 0.432169+0.203233im and 0.432169-0.203233im. If you try to optimize that you'll get an error. I think it's much more helpful for the user to get the error when they write a nonsensical constraint, than letting them figure out where the error was after they tried to solve and got that a PrimalInconsistent.

In my example the problem was solved correctly, but a lot of redundant constraints were added, and the dual variable recovered is incorrect.

[odow]

If you try to optimize that you'll get an error

It will tell you the problem is infeasible. That's not quite the same as an error.

when they write a nonsensical constraint

It's not JuMP's job to decide what a nonsensical constraint is. In some cases, you might want to add A == B to begin with, and then modify it later.

the dual variable recovered is incorrect

Can you clarify? Is the dual still wrong after we fixed #3797?

[araujoms]

Yes, the dual is still wrong. Just run the example I pasted, it returns the dual variable. To get the correct one you need to write noisy_state = Symmetric(ρ + λ * I(d^2)) instead of noisy_state = Hermitian(ρ + λ * I(d^2)).

[odow]

What makes you think the dual is wrong?

I get:

julia> dot(ρ, W), objective_value(model)
(0.03661524167501384, 0.03661524512473435)

Now, the Hermitian dual is not symmetric, but that's because we're not adding symmetric constraints.

You could fix it with (W + W') / 2.

We don't add symmetric constraints because we loose the type of the matrix:

julia> σ
4×4 Symmetric{VariableRef, Matrix{VariableRef}}:
 σ[1,1]  σ[1,2]  σ[1,3]  σ[1,4]
 σ[1,2]  σ[2,2]  σ[2,3]  σ[2,4]
 σ[1,3]  σ[2,3]  σ[3,3]  σ[3,4]
 σ[1,4]  σ[2,4]  σ[3,4]  σ[4,4]

julia> noisy_state
4×4 Hermitian{AffExpr, Matrix{AffExpr}}:
 λ + 0.0007142349938333907  -0.005378044977694615    …  -0.02486244188463285
 -0.005378044977694615      λ + 0.04049559043147808     0.18720915646155625
 0.008164586821029125       -0.06147768665340706        -0.2842083727379334
 -0.02486244188463285       0.18720915646155625         λ + 0.8654588781055151

julia> σ - noisy_state
4×4 Matrix{AffExpr}:
 σ[1,1] - λ - 0.0007142349938333907  …  σ[1,4] + 0.02486244188463285
 σ[1,2] + 0.005378044977694615          σ[2,4] - 0.18720915646155625
 σ[1,3] - 0.008164586821029125          σ[3,4] + 0.2842083727379334
 σ[1,4] + 0.02486244188463285           σ[4,4] - λ - 0.8654588781055151

I guess we could perhaps expect that -(::Symmetric, ::Hermitian) is Hermitian?

[odow]

So perhaps your underlying issue is:

julia> using JuMP, LinearAlgebra

julia> model = Model();

julia> @variable(model, A[1:2, 1:2] in PSDCone());

julia> @variable(model, B[1:2, 1:2] in HermitianPSDCone());

julia> A + B
2×2 Matrix{GenericAffExpr{ComplexF64, VariableRef}}:
 A[1,1] + real(B[1,1])                    A[1,2] + real(B[1,2]) + imag(B[1,2]) im
 A[1,2] + real(B[1,2]) - imag(B[1,2]) im  A[2,2] + real(B[2,2])

julia> C = LinearAlgebra.Symmetric(rand(Float64, 2, 2));

julia> D = LinearAlgebra.Hermitian(rand(ComplexF64, 2, 2));

julia> C + D
2×2 Hermitian{ComplexF64, Matrix{ComplexF64}}:
  1.66933+0.0im       0.626978+0.438058im
 0.626978-0.438058im  0.427326+0.0im

[odow]

We're missing overloads for these base methods: https://github.com/JuliaLang/julia/blob/48d4fd48430af58502699fdf3504b90589df3852/stdlib/LinearAlgebra/src/symmetric.jl#L516-L519

[araujoms]

Fair enough. It's not the design I would have chosen, but there's no silent failure anymore, which is the part I found really objectionable. Now it gives either the correct answer or PrimalInconsistent. I'm just a bit puzzled about why the eltype of the return dual variable is suddenly ComplexF64.

[odow]

I'm just a bit puzzled about why the eltype of the return dual variable is suddenly ComplexF64.

Because we assumed that Hermitian matrices are Complex-valued.

I've added a new method to #3805 where we now add Real-valued Hermitian matrices as SymmetricMatrixShape instead.

[araujoms]

Cool, now it's real. Thanks a lot!

[blegat]

Thanks for catching the bugs, issues in the dual are not easy to spot!