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
This tutorial was generated using Literate.jl. View the source .jl
file on GitHub.