Wram start of HermitianPSDCone in JuMP

Specific Domains Optimization (Mathematical)
1

I’m having trouble figuring out how to use set_start_value with HermitianPSDCone.

Here’s a minimal script based on the the quantum discrimination example in the docs:

using JuMP
using LinearAlgebra
import COSMO

ρs = [ [[1 0]; [0 0]],
       [[1 1]; [1 1]]/2 ]
ρs = map(Hermitian, ρs)

model = Model(COSMO.Optimizer)
set_silent(model)

@variable(model, E1[1:2, 1:2] in HermitianPSDCone())
@variable(model, E2[1:2, 1:2] in HermitianPSDCone())
@constraint(model, E1+E2 == LinearAlgebra.I)
@objective(model, Max, real(dot(ρs, [E1,E2])) / 2)

set_start_value.(E1, Hermitian(I(2)/2))
set_start_value.(E2, Hermitian(I(2)/2))

optimize!(model)
solution_summary(model)
println("value: ", objective_value(model))
println("true: ", 0.5 + 0.25 * sum(svdvals(ρs[1] - ρs[2])))

The set_start_value line gives the error:

MethodError: no method matching set_start_value(::GenericAffExpr{ComplexF64, VariableRef}, ::Float64)

I wasn’t sure if this was a bug or if there’s something special I need to do to set the initial value of a hermitian variable.

2

Hi @akirakyle, welcome to the forum!

This is a bug (really, a missing feature). I’ve opened an issue: set_start_value for complex variables · Issue #3550 · jump-dev/JuMP.jl · GitHub

For now, you can set a start value using the start keyword in @variable:

julia> using JuMP, LinearAlgebra

julia> model = Model();

julia> E1_start = Hermitian(I(2)/2)
2×2 Hermitian{Float64, Diagonal{Float64, Vector{Float64}}}:
 0.5  0.0
 0.0  0.5

julia> @variable(model, E1[i=1:2, j=1:2] in HermitianPSDCone(), start=E1_start[i,j])
2×2 Hermitian{GenericAffExpr{ComplexF64, VariableRef}, Matrix{GenericAffExpr{ComplexF64, VariableRef}}}:
 real(E1[1,1])                     real(E1[1,2]) + imag(E1[1,2]) im
 real(E1[1,2]) - imag(E1[1,2]) im  real(E1[2,2])

julia> all_variables(model)
4-element Vector{VariableRef}:
 real(E1[1,1])
 real(E1[1,2])
 real(E1[2,2])
 imag(E1[1,2])

julia> start_value.(all_variables(model))
4-element Vector{Float64}:
 0.5
 0.0
 0.5
 0.0
3

Thanks @odow for replying to me! Good to know that this was a bug and that I can use the start keyword.
However it seems to now error when calling optimize! with COSMO here’s what I ran:

using JuMP, LinearAlgebra, COSMO

model = Model(COSMO.Optimizer)
set_silent(model)

ρs = [ [[1 0]; [0 0]],
       [[1 1]; [1 1]]/2 ]
ρs = map(Hermitian, ρs)

E_ = Hermitian(I(2)/2)

@variable(model, E1[i=1:2, j=1:2] in HermitianPSDCone(), start = E_[i,j])
@variable(model, E2[i=1:2, j=1:2] in HermitianPSDCone(), start = E_[i,j])
@constraint(model, E1+E2 == LinearAlgebra.I)
@objective(model, Max, real(dot(ρs, [E1,E2])) / 2)

optimize!(model)
solution_summary(model)

Then I get the error:

MathOptInterface.UnsupportedAttribute{MathOptInterface.VariablePrimalStart}: Attribute MathOptInterface.VariablePrimalStart() is not supported by the model.

However this seems to be the result of using COSMO as it seems to work with SCS. I’m not sure if this is also a bug?

1 Like
4

Thanks @odow for replying to me!

No problem :smile:

However this seems to be the result of using COSMO as it seems to work with SCS. I’m not sure if this is also a bug?

Hmm. That shouldn’t happen. Let me take a look.

5

Okay, this is indeed a bug: Starting value in unsupported bridge errors instead of skipping · Issue #2327 · jump-dev/MathOptInterface.jl · GitHub.

Here’s what is happening:

  • If a solver does not support HermitianPSDCone explicitly, then JuMP will reformulate it into an equivalent PositiveSemidefiniteConeTriangle constraint through a bridge.
  • When you set the starting value of E1 and E2, the bridge needs to convert those starting values into equivalent starting values for the reformulated variables inside the solver
  • Some bridges (such as the hermitian to PSD bridge) do not yet have support for mapping starting values (Add support for starting values in bridges · Issue #684 · jump-dev/MathOptInterface.jl · GitHub)
  • Since starting values are an optional hint that can improve performance but not correctness, JuMP should silently ignore them if they are not supported by the solver.
  • COSMO and SCS supports starting values, but the Hermitan to PSD bridge does not, so JuMP should skip copying the starting values.
  • JuMP correctly skips the values for SCS (even though it doesn’t warn the user, this is a choice we made), but it errors for COSMO.
  • I don’t know why it errors, which is the bug.

So the conclusion is: don’t use starting values for now if you are using HermitianPSDCone :cry:

6

The starting values are implemented in Implement starting values for Variable.Hermitian bridge by blegat · Pull Request #2330 · jump-dev/MathOptInterface.jl · GitHub
You can try it out with:

using Pkg
pkg"add MathOptInterface#bl/start_herm"

restart Julia for the changes to take effect.

1 Like
7

Thanks @odow and @blegat for all the fast fixes related to this issue, I’m really quite impressed how quickly this has been addressed!

I’ve been trying to follow along the linked github issues and PRs as best I can but I may not be understanding exactly what’s been fixed yet.

The example I’ve posted now works for me with the latest master versions of JuMP and MathOptInterface when using the start = keyword of @variable. I can even get set_start_value working when using real starting values like so: set_start_value.(real.(E2), real.(E_)). However I can’t seem to figure out how to use set_starting_value with a hermitian starting value that has imaginary components. For example

using JuMP, LinearAlgebra, COSMO

model = Model(COSMO.Optimizer)
set_silent(model)

ρs = [ [[1 0]; [0 0]],
       [[1 1]; [1 1]]/2 ]

E_ = [[1 1im]; [-1im 1]]

@variable(model, E1[1:2, 1:2] in HermitianPSDCone())
@variable(model, E2[1:2, 1:2] in HermitianPSDCone())
@constraint(model, E1+E2 == LinearAlgebra.I)
@objective(model, Max, real(dot(ρs, [E1,E2])) / 2)

set_start_value.(real.(E2), real.(E_))
set_start_value.(imag.(E2), imag.(E_))

optimize!(model)
solution_summary(model)

gives me the error

Cannot call set_start_value with 0 because it is not an affine expression of one variable.

Using the start= keyword instead works just fine for this example so I can just use that for now, but I assume for repeatedly solving same problem with different coefficients, it would be more efficient to update the start value of with set_start_value rather than having to remove then re-add the variable?

1 Like
8

The issue is that the off-diagonal entries in E1 and E2 are actually the same variable. So there is a problem setting the imag(E2[2, 1]) element, which is actually -imag(E2[1, 2])

I didn’t test, but you could do something like this:

function set_hermitian_start(x, start)
    for j in 1:size(x, 2), i in 1:size(x, 1)
        set_start_value(real(x[i, j]), real(start[i, j]))
        if i < j
            set_start_value(imag(x[i, j]), imag(start[i, j]))
        end
    end
    return
end
set_hermitian_start(E2, E_)
9

The other issue is that there is no imaginary variable associated to the diagonal entries, so imag returns 0 hence the error message saying that you cannot call set_start_value with `0 .

10

@blegat, I guess we could handle that case if the start is also 0?

11

Yes, setting a starting value to a constant could be ignored, especially if they are the same.

12

@odow the set_hermitian_start function works for me. Thanks for explaining why what I was doing wasn’t working!

1 Like