Example: quantum state discrimination

This tutorial was generated using Literate.jl. Download the source as a .jl file.

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 uses the following packages:

using JuMP
import Clarabel
import LinearAlgebra

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.9049983143615443 + 0.0im -0.11984167294191471 + 0.2224268161763913im; -0.11984167294191471 - 0.2224268161763913im 0.09500168563845571 + 0.0im]
 [0.2503060704439256 + 0.0im -0.16997258971668044 - 0.09226624832975384im; -0.16997258971668044 + 0.09226624832975384im 0.7496939295560743 + 0.0im]

JuMP formulation

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

model = Model(Clarabel.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] + _[4] im; _[2] - _[4] im _[3]]
 [_[5] _[6] + _[8] im; _[6] - _[8] im _[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) == LinearAlgebra.I)

\[ \begin{bmatrix} {\_}_{1} + {\_}_{5} - 1 & {\_}_{2} + {\_}_{6} + {\_}_{4} im + {\_}_{8} im\\ {\_}_{2} + {\_}_{6} - {\_}_{4} im - {\_}_{8} im & {\_}_{3} + {\_}_{7} - 1\\ \end{bmatrix} \in \text{Zeros()} \]

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.45249915718077216 {\_}_{1} - 0.11984167294191471 {\_}_{2} + 0.2224268161763913 {\_}_{4} + 0.047500842819227854 {\_}_{3} + 0.1251530352219628 {\_}_{5} - 0.16997258971668044 {\_}_{6} - 0.09226624832975384 {\_}_{8} + 0.37484696477803714 {\_}_{7} \]

Now we optimize:

optimize!(model)
assert_is_solved_and_feasible(model)
solution_summary(model)

solution_summary(; result = 1, verbose = false)
├ solver_name          : Clarabel
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ └ raw_status         : SOLVED
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : FEASIBLE_POINT
│ ├ objective_value      : 8.64063e-01
│ └ dual_objective_value : 8.64063e-01
└ Work counters
  ├ solve_time (sec)   : 1.30670e-03
  └ barrier_iterations : 8

The probability of guessing correctly is:

objective_value(model)

0.8640627425181733

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

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

0.8640627582954737

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.9496066849311735 + 0.0im 0.034413986810280475 + 0.21603121721998986im; 0.034413986810280475 - 0.21603121721998986im 0.050393307139617316 + 0.0im]
 [0.05039331506882663 + 0.0im -0.034413986810280475 - 0.21603121721998986im; -0.034413986810280475 + 0.21603121721998986im 0.9496066928603828 + 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(Clarabel.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] + _[4] im; _[2] - _[4] im _[3]]
 [-_[1] + 1 -_[2] - _[4] im; -_[2] + _[4] im -_[3] + 1]

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

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

\[ 0.32734612195880936 {\_}_{1} + 0.050130916774765735 {\_}_{2} + 0.31469306450614515 {\_}_{4} - 0.3273461219588093 {\_}_{3} + 0.49999999999999994 \]

Then we can check that we get the same solution:

optimize!(model)
assert_is_solved_and_feasible(model)
solution_summary(model)

solution_summary(; result = 1, verbose = false)
├ solver_name          : Clarabel
├ Termination
│ ├ termination_status : OPTIMAL
│ ├ result_count       : 1
│ └ raw_status         : SOLVED
├ Solution (result = 1)
│ ├ primal_status        : FEASIBLE_POINT
│ ├ dual_status          : FEASIBLE_POINT
│ ├ objective_value      : 8.64063e-01
│ └ dual_objective_value : 8.64063e-01
└ Work counters
  ├ solve_time (sec)   : 1.27033e-03
  └ barrier_iterations : 9

objective_value(model)

0.8640627545817933