Can't constrain affine function to be Hermitian PSD

[araujoms]

If I declare some variable a, and then b = k - a, where k is a constant, I can't constrain b to be Hermitian PSD. I don't get an error message, just wrong results. This only happens with Hypatia, not other solvers such as SCS or SeDuMi, so I assume it must be related to the recent addition of native support for Hermitian PSD constraints in #822. Now I don't know if the bug is in MOI or Hypatia, I'm guessing it is here.

The code below reproduces the problem. For d=2 you see that the variable E2 is not PSD, and for higher d you just get PrimalInfeasible.

function random_state(d)
    x = randn(ComplexF64, (d,d))
    y = x*x'
    rho = Hermitian(y/tr(y))
    return rho
end
	
function discriminate(d)
    model = Model(Hypatia.Optimizer)
    rho = [random_state(d) for i in 1:2]
    E1 = @variable(model, [1:d, 1:d] in HermitianPSDCone())
    E2 = I-E1
    @constraint(model, E2 in HermitianPSDCone())
    obj = real(dot(rho[1],E1) + dot(rho[2],E2))/2
    @objective(model, Max, obj)
    optimize!(model)
    println(value(obj))
    println(0.5 + 0.25*sum(svdvals(rho[1]-rho[2])))
    display(eigvals(value.(E1)))
    display(eigvals(value.(E2)))
end

[araujoms]

@odow can you tell whether the bug is in MOI or Hypatia?

[odow]

I'll take a look. But just as a note for future reference, you can improve the reproducibility of your examples by including the using JuMP, Hypatia, LinearAlgebra line, and by pasting the output of actually running this in the REPL.

[odow]

julia> using JuMP, Hypatia, LinearAlgebra

julia> function random_state(d)
           x = randn(ComplexF64, (d,d))
           y = x*x'
           rho = Hermitian(y/tr(y))
           return rho
       end
random_state (generic function with 1 method)

julia> function discriminate(d)
           model = Model(Hypatia.Optimizer)
           rho = [random_state(d) for i in 1:2]
           E1 = @variable(model, [1:d, 1:d] in HermitianPSDCone())
           E2 = I-E1
           @constraint(model, E2 in HermitianPSDCone())
           obj = real(dot(rho[1],E1) + dot(rho[2],E2))/2
           @objective(model, Max, obj)
           optimize!(model)
           @show value(obj)
           @show 0.5 + 0.25*sum(svdvals(rho[1]-rho[2]))
           @show eigvals(value.(E1))
           @show eigvals(value.(E2))
           return
       end
discriminate (generic function with 1 method)

julia> discriminate(2)

 iter        p_obj        d_obj |  abs_gap    x_feas    z_feas |      tau       kap        mu | dir_res     prox  step     alpha
    0  -5.5910e-01  -2.5301e+00 | 4.00e+00  2.22e-01  2.65e-01 | 1.00e+00  1.00e+00  1.00e+00 |
    1  -7.1917e-01  -1.0292e+00 | 7.64e-01  3.28e-02  3.91e-02 | 1.35e+00  1.74e-01  2.00e-01 | 2.2e-16  8.0e-01  co-a  8.00e-01
    2  -9.2930e-01  -9.8060e-01 | 9.19e-02  5.98e-03  7.13e-03 | 1.11e+00  3.20e-02  2.55e-02 | 1.5e-16  7.6e-01  co-a  8.50e-01
    3  -9.3984e-01  -9.4575e-01 | 1.12e-02  5.62e-04  6.71e-04 | 1.18e+00  1.91e-03  2.69e-03 | 2.9e-15  6.9e-01  co-a  9.00e-01
    4  -9.4322e-01  -9.4387e-01 | 1.03e-03  6.16e-05  7.35e-05 | 1.08e+00  1.89e-04  2.46e-04 | 1.2e-15  5.6e-01  co-a  9.00e-01
    5  -9.4362e-01  -9.4365e-01 | 4.75e-05  3.38e-06  4.03e-06 | 9.85e-01  8.67e-06  1.12e-05 | 1.6e-14  9.0e-01  co-a  9.50e-01
    6  -9.4364e-01  -9.4365e-01 | 4.04e-06  3.85e-07  4.59e-07 | 8.65e-01  1.02e-06  9.84e-07 | 5.8e-13  5.4e-01  co-a  9.00e-01
    7  -9.4364e-01  -9.4364e-01 | 3.86e-07  4.10e-08  4.89e-08 | 8.12e-01  9.34e-08  9.24e-08 | 1.6e-12  5.4e-01  co-a  9.00e-01
    8  -9.4364e-01  -9.4364e-01 | 1.75e-08  2.29e-09  2.69e-09 | 7.38e-01  4.80e-09  4.20e-09 | 4.2e-11  8.0e-01  co-a  9.50e-01
    9  -9.4364e-01  -9.4364e-01 | 1.60e-09  2.49e-10  2.97e-10 | 6.69e-01  4.29e-10  3.77e-10 | 7.4e-11  9.6e-01  co-a  9.00e-01
optimal solution found; terminating

status is Optimal after 9 iterations and 0.006 seconds

value(obj) = 0.9436436459949311
0.5 + 0.25 * sum(svdvals(rho[1] - rho[2])) = 0.8195179862102644
eigvals(value.(E1)) = [1.2997080495180668e-9, 1.39263144270254]
eigvals(value.(E2)) = [-0.3926314427025398, 0.9999999987002921]

julia> discriminate(3)

 iter        p_obj        d_obj |  abs_gap    x_feas    z_feas |      tau       kap        mu | dir_res     prox  step     alpha
    0  -3.9204e-01  -1.3409e+00 | 6.00e+00  1.39e-01  4.57e-01 | 1.00e+00  1.00e+00  1.00e+00 |
    1  -3.8423e-01   1.7077e-01 | 1.13e+00  6.49e-02  2.13e-01 | 4.28e-01  6.27e-01  2.00e-01 | 2.2e-16  7.9e-01  co-a  8.00e-01
    2  -4.0590e-01   1.2378e+01 | 1.75e-01  1.12e-01  3.69e-01 | 3.71e-02  5.33e-01  2.78e-02 | 3.8e-15  5.3e-01  co-a  8.50e-01
    3  -3.7882e-01   7.3046e+01 | 1.52e-02  5.81e-02  1.91e-01 | 7.18e-03  5.33e-01  2.72e-03 | 4.3e-16  8.1e-01  co-a  9.00e-01
    4  -4.1657e-01   7.9823e+02 | 2.24e-03  1.15e-01  3.79e-01 | 5.43e-04  4.34e-01  3.54e-04 | 1.9e-14  5.9e-01  co-a  8.50e-01
    5  -3.7099e-01   2.7371e+03 | 2.88e-04  5.81e-02  1.91e-01 | 1.62e-04  4.42e-01  5.13e-05 | 2.8e-16  8.3e-01  co-a  8.50e-01
    6  -4.3619e-01   3.0030e+04 | 4.22e-05  1.16e-01  3.81e-01 | 1.21e-05  3.64e-01  6.66e-06 | 9.4e-15  6.4e-01  co-a  8.50e-01
    7  -3.5277e-01   1.0043e+05 | 5.36e-06  5.74e-02  1.89e-01 | 3.68e-06  3.70e-01  9.60e-07 | 2.2e-16  9.3e-01  co-a  8.50e-01
    8  -4.9956e-01   1.0214e+06 | 7.66e-07  1.09e-01  3.58e-01 | 2.90e-07  2.97e-01  1.22e-07 | 1.8e-14  7.6e-01  co-a  8.50e-01
    9  -3.0777e-01   4.1645e+06 | 9.55e-08  6.37e-02  2.09e-01 | 7.46e-08  3.11e-01  1.69e-08 | 5.0e-16  8.4e-01  co-a  8.50e-01
   10  -4.9738e-01   1.8165e+07 | 2.76e-08  9.69e-02  3.18e-01 | 1.47e-08  2.67e-01  4.50e-09 | 3.3e-15  7.1e-01  co-a  7.00e-01
   11  -2.8480e-01   8.5288e+07 | 3.64e-09  7.03e-02  2.31e-01 | 3.04e-09  2.59e-01  6.33e-10 | 2.0e-16  8.3e-01  co-a  8.50e-01
   12  -5.4760e-01   5.9478e+08 | 7.09e-10  1.05e-01  3.46e-01 | 4.07e-10  2.42e-01  1.15e-10 | 1.0e-15  9.6e-01  co-a  8.00e-01
   13  -3.5192e-01   2.4474e+09 | 8.54e-11  6.91e-02  2.27e-01 | 9.28e-11  2.27e-01  1.52e-11 | 4.8e-16  8.3e-01  co-a  8.50e-01
   14  -4.5819e-01   1.4968e+10 | 1.81e-11  1.02e-01  3.34e-01 | 1.26e-11  1.89e-01  2.92e-12 | 9.0e-16  8.1e-01  co-a  8.00e-01
   15  -3.4630e-01   6.9563e+10 | 2.29e-12  7.05e-02  2.32e-01 | 2.73e-12  1.90e-01  4.02e-13 | 2.9e-16  6.5e-01  co-a  8.50e-01
primal infeasibility detected; terminating

status is PrimalInfeasible after 15 iterations and 0.011 seconds

value(obj) = 0.4999999999995805
0.5 + 0.25 * sum(svdvals(rho[1] - rho[2])) = 0.6793815656525367
eigvals(value.(E1)) = [-1.929835268936554e-12, 6.517232220696772e-13, 1.7843677078198655e-12]
eigvals(value.(E2)) = [0.9999999999982159, 0.9999999999993484, 1.00000000000193]

This just seems like a straight up bug.

[odow]

Reading: #822

well this was a bit of a pain because MOI chose a different definition of the Hermitian PSD cone to what we had implemented
anyway, there's a chance a bug slipped in here, @araujoms, so please be vigilant as you try out using the JuMP HermitianPSDCone with Hypatia

So yeah. Probably a bug somewhere.

[araujoms]

Thanks for taking a look. I have no doubt that it is a bug, I'm just unsure whether it is in Hypatia or MOI. I thought it could be in MOI because this code uses a constraint MOI.VectorAffineFunction{Float64}-in-MOI.HermitianPositiveSemidefiniteConeTriangle, that as far as I can see is never used otherwise, as all other solvers bridge it.

Sorry for leaving out the preamble and output. I was trying to be brief, but on a second thought brevity is of no benefit here.

[lkapelevich]

@chriscoey will get this a lot faster but I can take a look if no one can get around to it.
It looks like Hypatia's internal model.h ended up with scaling in the wrong indices.

[chriscoey]

fixed by #826