Quantum state discrimination

This tutorial solves the problem of quantum state discrimination.

The purpose of this tutorial to demonstrate how to solve problems involving complex-valued decision variables and the HermitianPSDCone. See Complex number support for more details.

Required packages

This tutorial makes use of the following packages:

using JuMP
import LinearAlgebra
import SCS

Formulation

A d-dimensional quantum state, $\rho$, can be defined by a complex-valued Hermitian matrix with a trace of 1. Assume we have N d-dimensional quantum states, $\{\rho_i\}_{i=1}^N$, each of which is equally likely.

The goal of the quantum state discrimination problem is to choose a positive operator-valued measure (POVM), $\{ E_i \}_{i=1}^N$, such that if we observe $E_i$ then the most probable state that we are in is $\rho_i$.

Each POVM element, $E_i$, is a complex-valued Hermitian matrix, and there is a requirement that $\sum\limits_{i=1}^N E_i = \mathbf{I}$.

To choose a POVM, we want to maximize the probability that we guess the quantum state correctly. This can be formulated as the following optimization problem:

\[\begin{aligned} \max\limits_{E} \;\; & \frac{1}{N} \sum\limits_{i=1}^N \operatorname{tr}(\rho_i E_i) \\ \text{s.t.} \;\; & \sum\limits_{i=1}^N E_i = \mathbf{I} \\ & E_i \succeq 0 \; \forall i = 1,\ldots,N. \end{aligned}\]

Data

To setup our problem, we need N d-dimensional quantum states. To keep the problem simple, we use N = 2 and d = 2.

N, d = 2, 2

(2, 2)

We then generated N random d-dimensional quantum states:

function random_state(d)
    x = randn(ComplexF64, (d, d))
    y = x * x'
    return LinearAlgebra.Hermitian(y / LinearAlgebra.tr(y))
end

ρ = [random_state(d) for i in 1:N]

2-element Vector{LinearAlgebra.Hermitian{ComplexF64, Matrix{ComplexF64}}}:
 [0.9317077745932991 + 0.0im 0.17785044289093424 - 0.11927897015052207im; 0.17785044289093424 + 0.11927897015052207im 0.06829222540670073 + 0.0im]
 [0.22453156773914584 + 0.0im 0.061840181904862485 - 0.052357134831832484im; 0.061840181904862485 + 0.052357134831832484im 0.7754684322608543 + 0.0im]

JuMP formulation

To model the problem in JuMP, we need a solver that supports positive semidefinite matrices:

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

Then, we construct our set of E variables:

E = [@variable(model, [1:d, 1:d] in HermitianPSDCone()) for i in 1:N]

2-element Vector{LinearAlgebra.Hermitian{GenericAffExpr{ComplexF64, VariableRef}, Matrix{GenericAffExpr{ComplexF64, VariableRef}}}}:
 [_[1] _[2] + (0.0 + 1.0im) _[4]; _[2] + (0.0 - 1.0im) _[4] _[3]]
 [_[5] _[6] + (0.0 + 1.0im) _[8]; _[6] + (0.0 - 1.0im) _[8] _[7]]

Here we have created a vector of matrices. This is different to other modeling languages such as YALMIP, which allow you to create a multi-dimensional array in which 2-dimensional slices of the array are Hermitian matrices.

We also need to enforce the constraint that $\sum\limits_{i=1}^N E_i = \mathbf{I}$:

@constraint(model, sum(E[i] for i in 1:N) .== LinearAlgebra.I)

2×2 Matrix{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{ComplexF64}, MathOptInterface.EqualTo{ComplexF64}}, ScalarShape}}:
 _[1] + _[5] = 1.0 - 0.0im                                            …  _[2] + (0.0 + 1.0im) _[4] + _[6] + (0.0 + 1.0im) _[8] = 0.0 - 0.0im
 _[2] + (0.0 - 1.0im) _[4] + _[6] + (0.0 - 1.0im) _[8] = 0.0 + 0.0im     _[3] + _[7] = 1.0 - 0.0im

This constraint is a complex-valued equality constraint. In the solver, it will be decomposed onto two types of equality constraints: one to enforce equality of the real components, and one to enforce equality of the imaginary components.

Our objective is to maximize the expected probability of guessing correctly:

@objective(
    model,
    Max,
    sum(real(LinearAlgebra.tr(ρ[i] * E[i])) for i in 1:N) / N,
)

\[ 0.46585388729664956 {\_}_{1} + 0.17785044289093424 {\_}_{2} - 0.11927897015052207 {\_}_{4} + 0.03414611270335036 {\_}_{3} + 0.11226578386957292 {\_}_{5} + 0.061840181904862485 {\_}_{6} - 0.052357134831832484 {\_}_{8} + 0.38773421613042713 {\_}_{7} \]

Now we optimize:

optimize!(model)
solution_summary(model)

* Solver : SCS

* Status
  Result count       : 1
  Termination status : OPTIMAL
  Message from the solver:
  "solved"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 8.59873e-01
  Dual objective value : 8.59872e-01

* Work counters
  Solve time (sec)   : 6.49411e-04

The probability of guessing correctly is:

objective_value(model)

0.8598731680110233

When N = 2, there is a known analytical solution of:

0.5 + 0.25 * sum(LinearAlgebra.svdvals(ρ[1] - ρ[2]))

0.859873276389449

proving that we found the optimal solution.

Finally, the optimal POVM is:

solution = [value.(e) for e in E]

2-element Vector{Matrix{ComplexF64}}:
 [0.9912673599117596 - 0.0im 0.08059102584033337 - 0.04648984764809457im; 0.08059102584033337 + 0.04648984764809457im 0.00873261686626009 - 0.0im]
 [0.008732639556408501 - 0.0im -0.0805910259158836 + 0.04648984770180667im; -0.0805910259158836 - 0.04648984770180667im 0.9912673827839558 - 0.0im]

Alternative formulation

The formulation above includes N Hermitian matrices and a set of linear equality constraints. We can simplify the problem by replacing $E_N$ with $E_N = I - \sum\limits_{i=1}^{N-1} E_i$. This results in:

model = Model(SCS.Optimizer)
set_silent(model)
E = [@variable(model, [1:d, 1:d] in HermitianPSDCone()) for i in 1:N-1]
E_N = LinearAlgebra.Hermitian(LinearAlgebra.I - sum(E))
@constraint(model, E_N in HermitianPSDCone())
push!(E, E_N)

2-element Vector{LinearAlgebra.Hermitian{GenericAffExpr{ComplexF64, VariableRef}, Matrix{GenericAffExpr{ComplexF64, VariableRef}}}}:
 [_[1] _[2] + (0.0 + 1.0im) _[4]; _[2] + (0.0 - 1.0im) _[4] _[3]]
 [(-1.0 - 0.0im) _[1] + (1.0 - 0.0im) (-1.0 - 0.0im) _[2] + (-0.0 - 1.0im) _[4]; (-1.0 + 0.0im) _[2] + (-0.0 + 1.0im) _[4] (-1.0 - 0.0im) _[3] + (1.0 - 0.0im)]

The objective can also be simplified, by observing that it is equivalent to:

@objective(model, Max, real(LinearAlgebra.dot(ρ, E)) / N)

\[ 0.35358810342707664 {\_}_{1} + 0.11601026098607176 {\_}_{2} - 0.06692183531868959 {\_}_{4} - 0.35358810342707675 {\_}_{3} + 0.5 \]

Then we can check that we get the same solution:

optimize!(model)
solution_summary(model)

* Solver : SCS

* Status
  Result count       : 1
  Termination status : OPTIMAL
  Message from the solver:
  "solved"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 8.59873e-01
  Dual objective value : 8.59872e-01

* Work counters
  Solve time (sec)   : 6.24710e-04

objective_value(model)

0.859873143250681