Error when constraining variable to be Hermitian PSD

Specific Domains Optimization (Mathematical)
,
1

I’m writing an SDP for a traditional problem in quantum information, state discrimination. When I constrain the variables to real symmetric, it works. But the problem really needs them to be complex Hermitian. The code is

using JuMP
import LinearAlgebra
import Random
import SCS

function random_state(d)
	x = randn(ComplexF64, (d,d))
	y = x*adjoint(x)
	rho = Hermitian(y/tr(y))
	return rho
end
	
function discriminate(d,n)
    model = Model(SCS.Optimizer)
    rho = Array{ComplexF64}(undef,d, d, n)
    for i=1:n
    	rho[:,:,i] = random_state(d)
    end
	E = model[:E] = reshape(
           hcat([
               @variable(model, [1:d, 1:d] in PSDCone(), base_name = "E$(i)")
               for i in 1:n-1
           ]...),
           d, d, n-1,
       )
	lastE = I - sum(E[:,:,i] for i in 1:n-1)
	@constraint(model, lastE in PSDCone())
	obj = real(sum(dot(rho[:,:,i],E[:,:,i]) for i in 1:n-1)+dot(rho[:,:,n],lastE))/n
    @objective(model, Max, obj)
    optimize!(model)
    print(value(obj))
end

When I change “PSDCone()” to “HermitianPSDCone()” I get the error “Unrecognized constraint building format.” Is this a bug? Is there a way to make it work?

A bit unrelated, but is there a more convenient way to declare a set of Hermitian variables? With YALMIP I can just do

E = sdpvar(d,d,n-1,'hermitian','complex');

which is much more concise and readable.

2

I should improve the documentation, Complex number support · JuMP, to explain that H needs to be explicitly marked as Hermitian. (See [docs] clarify that matrix must be Hermitian by odow · Pull Request #3241 · jump-dev/JuMP.jl · GitHub)

Here’s how I would write your model:

using JuMP
import LinearAlgebra
import SCS

function random_state(d)
    x = randn(ComplexF64, (d, d))
    y = x * x'
    return LinearAlgebra.Hermitian(y / LinearAlgebra.tr(y))
end
    
function discriminate(d, n)
    model = Model(SCS.Optimizer)
    rho = [random_state(d) for i in 1:n]
    E = [
        @variable(model, [1:d, 1:d] in HermitianPSDCone(), base_name = "E$(i)")
        for i in 1:n-1
    ]
    lastE = LinearAlgebra.I - sum(E)
    @constraint(model, LinearAlgebra.Hermitian(lastE) in HermitianPSDCone())
    obj = real(
        sum(LinearAlgebra.dot(rho[i], E[i]) for i in 1:n-1) +
        LinearAlgebra.dot(rho[n], lastE)
    ) / n
    @objective(model, Max, obj)
    optimize!(model)
    return objective_value(model)
end
3

Thanks a lot! Something very weird happened, though. If instead of your version

	lastE = I - sum(E[i] for i in 1:n-1)
	@constraint(model, Hermitian(lastE) in HermitianPSDCone())

I write

	lastE = Hermitian(I - sum(E[i] for i in 1:n-1))
	@constraint(model, lastE in HermitianPSDCone())

I don’t get an error message either, but the SDP gives me the wrong result! Any idea what is going on?

And, since you’re improving the documentation anyway, perhaps it would be worth mentioning this much nicer way of declaring a set of Hermitian PSD variables that you wrote down?

4

Hmm. Let me take a look.

You could also write it like this, which is perhaps even simpler:

function discriminate(d, n)
    model = Model(SCS.Optimizer)
    set_silent(model)
    rho = [random_state(d) for i in 1:n]
    E = [@variable(model, [1:d, 1:d] in HermitianPSDCone()) for i in 1:n]
    @constraint(model, sum(E) .== LinearAlgebra.I)
    @objective(model, Max, real(LinearAlgebra.dot(rho, E)) / n)
    optimize!(model)
    return objective_value(model)
end
5

I’m aware of this formulation, it’s the most straightforward way to code state discrimination. I avoided it because it has unnecessary variables and constraints, and as such is slower and less precise

1 Like
6

It looks like your incorrect case is a bug in JuMP: Incorrect summation with Hermitian matrices · Issue #3242 · jump-dev/JuMP.jl · GitHub

We’ve only recently added support for HermitianPSDCone and the associated Complex support, but this is a pretty major nasty on our part. It took me a while to find!

cc @blegat

7

Wow, you tracked down and fixed the bug on a Sunday. That’s some amazing dedication, thanks!

May I suggest to add this code as an example in the documentation for complex numbers? There are no examples there, and this is a simple and well-known problem.

1 Like
8

Is there an accessible reference one could link to for this?

9

For @odow, it was already Monday in New-Zealand ^^
@araujoms Can you check that your model now works with JuMP#master ?
I agree it would make a nice tutorial, quantum information was one of the motivation for adding Complex number support. A reference would be helpful indeed.

1 Like
10

There’s always Wikipedia, but that article is not so nice. A better reference would be this review or Section 11.3 of Bergou’s book, that I tried to link but got blocked by the spam filter.

1 Like
11

Just tested with the git version, now the following code works:

function discriminate(d,n)
    model = Model(SCS.Optimizer)
    rho = [random_state(d) for i in 1:n]
	E = [@variable(model, [1:d, 1:d] in HermitianPSDCone()) for i in 1:n-1]
	lastE = Hermitian(I - sum(E[i] for i in 1:n-1))
	@constraint(model, lastE in HermitianPSDCone())
	obj = real(sum(dot(rho[i],E[i]) for i in 1:n-1)+dot(rho[n],lastE))/n
    @objective(model, Max, obj)
    optimize!(model)
    print(value(obj))
end
12

Thanks for the quick reply, this should be fixed in JuMP v1.8.2

13

@araujoms I’ve started a tutorial, but I’m not a domain expert so any comments/suggestions would be helpful: