Edit on GitHub

Examples

Let's load the package, an SDP solver, and define a simple plotting routine.

using GuessworkQuantumSideInfo
using SCS, Plots
get_sdp_solver() = SCSSolver(verbose=false)
plot_pmfN(data) = bar(pmfN(data); xlabel="Guess number", ylabel="Probability of guessing correctly", legend = false)
plot_pmfN (generic function with 1 method)

Next, we define some basic qubit states.

dB = 2
ketplus = (ket(1, dB) + ket(2,dB))/sqrt(2)
ketminus = (ket(1, dB) - ket(2,dB))/sqrt(2)
ketzero = ket(1, dB)
ketone = ket(2, dB)
2-element SparseArrays.SparseVector{Complex{Float64},Int64} with 1 stored entry:
  [2]  =  1.0+0.0im

Example 1: A warmup with trivial examples

Let's consider the case with $J=2$ and both states are the same. The side information is therefore completely uninformative, and the guesswork is

\[G(X|B) = \frac{1}{2}\cdot 1 + \frac{1}{2} \cdot 2 = 1.5.\]

We can check this:

julia> p = [0.5, 0.5]
2-element Array{Float64,1}:
 0.5
 0.5

julia> ρBs = dm.([ ketzero, ketzero  ])
2-element Array{Array{Complex{Float64},2},1}:
 [1.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im]
 [1.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im]

julia> output = guesswork(p, ρBs; solver = get_sdp_solver());

julia> output.optval
1.5000105244940878

We see the result agrees with 1.5, as we expected. Likewise, if we choose the two states as $ |0\rangle, |1\rangle$, we can get it in one guess every time, of course, since they are orthogonal:

julia> p = [0.5, 0.5]
2-element Array{Float64,1}:
 0.5
 0.5

julia> ρBs = dm.([ ketzero, ketone  ])
2-element Array{Array{Complex{Float64},2},1}:
 [1.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im]
 [0.0+0.0im 0.0+0.0im; 0.0+0.0im 1.0+0.0im]

julia> output = guesswork(p, ρBs; solver = get_sdp_solver());

julia> output.optval
0.9999976286547437

We can inspect the POVMs:

julia> output.Es
2-element Array{Array{Complex{Float64},2},1}:
 [1.00001+0.0im 0.0+0.0im; 0.0+0.0im -4.72791e-6+0.0im]
 [-4.72791e-6+0.0im 0.0+0.0im; 0.0+0.0im 1.00001+0.0im]

As we would expect, we (approximately) obtain the projection onto $|0 \rangle$ and the projection onto $|1 \rangle$.

Example 2: the BB84 states

Let's consider the four states $|+ \rangle, |-\rangle, |0\rangle, |1\rangle$.

julia> p = [0.25, 0.25, 0.25, 0.25]
4-element Array{Float64,1}:
 0.25
 0.25
 0.25
 0.25

julia> ρBs = dm.([ ketplus, ketminus, ketzero, ketone  ])
4-element Array{Array{Complex{Float64},2},1}:
 [0.5+0.0im 0.5+0.0im; 0.5+0.0im 0.5+0.0im]
 [0.5+0.0im -0.5-0.0im; -0.5+0.0im 0.5+0.0im]
 [1.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im]
 [0.0+0.0im 0.0+0.0im; 0.0+0.0im 1.0+0.0im]

julia> output = guesswork(p, ρBs; solver = get_sdp_solver());

julia> output.optval
1.7094306690270955
plot_pmfN(output)
1 2 3 4 0.0 0.1 0.2 0.3 0.4 0.5 Guess number Probability of guessing correctly

Let's try the same example, but imposing a steep cost for the fourth guess.

julia> c = [1.0, 2.0, 3.0, 5000.0]
4-element Array{Float64,1}:
    1.0
    2.0
    3.0
 5000.0

julia> output = guesswork(p, ρBs; c = c, solver = get_sdp_solver());

julia> output.optval
1.7499750650637333

We see that the average number of guesses to get a correct answer has gone up. However, inspecting the probability mass function for the number of guesses under the optimal strategy

plot_pmfN(output)
1 2 3 4 0.0 0.1 0.2 0.3 0.4 0.5 Guess number Probability of guessing correctly

we see that the probability that the probability of guessing correctly on the fourth guess goes to almost zero.

Example 3: two copies of the BB84 states

Let us consider two tensor copies of the BB84 states:

julia> p = ones(16)/16
16-element Array{Float64,1}:
 0.0625
 0.0625
 0.0625
 0.0625
 0.0625
 0.0625
 0.0625
 0.0625
 0.0625
 0.0625
 0.0625
 0.0625
 0.0625
 0.0625
 0.0625
 0.0625

julia> ρBs = iid_copies(BB84_states(), 2)
16-element Array{Array{Complex{Float64},2},1}:
 [1.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im]
 [0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 1.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im]
 [0.5+0.0im 0.0+0.0im -0.5-0.0im 0.0-0.0im; 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0-0.0im; -0.5+0.0im -0.0+0.0im 0.5+0.0im 0.0+0.0im; -0.0+0.0im -0.0+0.0im 0.0+0.0im 0.0+0.0im]
 [0.5+0.0im 0.0+0.0im 0.5+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.5+0.0im 0.0+0.0im 0.5+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im]
 [0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 1.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im]
 [0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 1.0+0.0im]
 [0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0-0.0im; 0.0+0.0im 0.5+0.0im 0.0-0.0im -0.5-0.0im; -0.0+0.0im -0.0+0.0im 0.0+0.0im 0.0+0.0im; -0.0+0.0im -0.5+0.0im 0.0+0.0im 0.5+0.0im]
 [0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.5+0.0im 0.0+0.0im 0.5+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.5+0.0im 0.0+0.0im 0.5+0.0im]
 [0.5+0.0im -0.5-0.0im 0.0+0.0im 0.0-0.0im; -0.5+0.0im 0.5+0.0im -0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0-0.0im 0.0+0.0im 0.0-0.0im; -0.0+0.0im 0.0+0.0im -0.0+0.0im 0.0+0.0im]
 [0.0+0.0im 0.0-0.0im 0.0+0.0im 0.0-0.0im; -0.0+0.0im 0.0+0.0im -0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0-0.0im 0.5+0.0im -0.5-0.0im; -0.0+0.0im 0.0+0.0im -0.5+0.0im 0.5+0.0im]
 [0.25+0.0im -0.25-0.0im -0.25-0.0im 0.25+0.0im; -0.25+0.0im 0.25+0.0im 0.25+0.0im -0.25-0.0im; -0.25+0.0im 0.25+0.0im 0.25+0.0im -0.25-0.0im; 0.25-0.0im -0.25+0.0im -0.25+0.0im 0.25+0.0im]
 [0.25+0.0im -0.25-0.0im 0.25+0.0im -0.25-0.0im; -0.25+0.0im 0.25+0.0im -0.25+0.0im 0.25+0.0im; 0.25+0.0im -0.25-0.0im 0.25+0.0im -0.25-0.0im; -0.25+0.0im 0.25+0.0im -0.25+0.0im 0.25+0.0im]
 [0.5+0.0im 0.5+0.0im 0.0+0.0im 0.0+0.0im; 0.5+0.0im 0.5+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im]
 [0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im; 0.0+0.0im 0.0+0.0im 0.5+0.0im 0.5+0.0im; 0.0+0.0im 0.0+0.0im 0.5+0.0im 0.5+0.0im]
 [0.25+0.0im 0.25+0.0im -0.25-0.0im -0.25-0.0im; 0.25+0.0im 0.25+0.0im -0.25-0.0im -0.25-0.0im; -0.25+0.0im -0.25+0.0im 0.25+0.0im 0.25+0.0im; -0.25+0.0im -0.25+0.0im 0.25+0.0im 0.25+0.0im]
 [0.25+0.0im 0.25+0.0im 0.25+0.0im 0.25+0.0im; 0.25+0.0im 0.25+0.0im 0.25+0.0im 0.25+0.0im; 0.25+0.0im 0.25+0.0im 0.25+0.0im 0.25+0.0im; 0.25+0.0im 0.25+0.0im 0.25+0.0im 0.25+0.0im]

In this case, there are $16! = 20922789888000$ possible guessing orders, and hence $16!$ variables in the primal formulation of the SDP, or $16!+1$ constraints in the dual form of the SDP. In either case, we can't even fit them all into our computer's memory. Instead, we resort to bounds:

julia> lb_output = guesswork_lower_bound(p, ρBs, solver = get_sdp_solver());

julia> lb_output.optval
2.4999988321870794

julia> ub_output = guesswork_upper_bound(p, ρBs; max_time = 30, make_solver = get_sdp_solver);

julia> ub_output.optval
3.7850304389673903

A closer look at guesswork_upper_bound

We can understand guesswork_upper_bound better by setting verbose=true. The algorithm computes a sequence of upper bounds by relaxing the dual problem. First, it removes all constraints, and chooses the dual variable $Y$ as the identity matrix. Then it uses a simulated annealing algorithm to heuristically minimize $\lambda_\text{min}(R_{\vec g} - Y)$ over the possible guessing orders $\vec g$ (which are permutations), to find a constraint that is "maximally" violated by this choice of $Y$. Then we add the corresponding constraint $Y \leq R_{\vec g}$ to the dual problem and solve it again. This is repeated until either a fixed number of constraints is added, some number of simulated annealing runs fails to find another violated constraint, or a time limit is reached. In the following, we set a time limit of 30 seconds.

julia> ub_output = guesswork_upper_bound(p, ρBs; verbose=true, max_time = 30, make_solver = get_sdp_solver);
┌ Info: Adding constraint 1 
│   violation = -15.04763530510284
│   upper_bound = 8.499998102437322
│   SA_time = 0.621630557
│   solve_time = 0.002665789
└   total_time_so_far = 0.621630557
┌ Info: Adding constraint 2 
│   violation = -2.1434304110950664
│   upper_bound = 6.109145197058703
│   SA_time = 0.712243777
│   solve_time = 0.003947884
└   total_time_so_far = 1.3338743339999999
┌ Info: Adding constraint 3 
│   violation = -1.2088573402823848
│   upper_bound = 4.900570958966567
│   SA_time = 0.666860672
│   solve_time = 0.00466978
└   total_time_so_far = 2.0007350059999998
┌ Info: Adding constraint 4 
│   violation = -1.0726936574886017
│   upper_bound = 3.827924218038989
│   SA_time = 0.669246768
│   solve_time = 0.006375173
└   total_time_so_far = 2.669981774
┌ Info: Adding constraint 5 
│   violation = -0.05409696222522099
│   upper_bound = 3.823675277468773
│   SA_time = 0.658922914
│   solve_time = 0.00696237
└   total_time_so_far = 3.328904688
┌ Info: Adding constraint 6 
│   violation = -0.007118057426946071
│   upper_bound = 3.8235711961433845
│   SA_time = 0.665609489
│   solve_time = 0.011104453
└   total_time_so_far = 3.994514177
┌ Info: Adding constraint 7 
│   violation = -0.027944430629653967
│   upper_bound = 3.8221360985039117
│   SA_time = 0.805764203
│   solve_time = 0.012375348
└   total_time_so_far = 4.80027838
┌ Info: Adding constraint 8 
│   violation = -0.019340699705592386
│   upper_bound = 3.8164187336703246
│   SA_time = 0.641728599
│   solve_time = 0.017942524
└   total_time_so_far = 5.442006979
┌ Info: Adding constraint 9 
│   violation = -0.039438844730537194
│   upper_bound = 3.8138184204432792
│   SA_time = 0.660692822
│   solve_time = 0.022773004
└   total_time_so_far = 6.102699801
┌ Info: Adding constraint 10 
│   violation = -0.019450788017574952
│   upper_bound = 3.812580967407884
│   SA_time = 0.655473448
│   solve_time = 0.022099407
└   total_time_so_far = 6.758173249
┌ Info: Adding constraint 11 
│   violation = -0.013354042454588617
│   upper_bound = 3.810294648288276
│   SA_time = 0.678446155
│   solve_time = 0.023253002
└   total_time_so_far = 7.436619404
┌ Info: Adding constraint 12 
│   violation = -0.006991120419104037
│   upper_bound = 3.8088446535796536
│   SA_time = 0.684400733
│   solve_time = 0.034521456
└   total_time_so_far = 8.121020137
┌ Info: Adding constraint 13 
│   violation = -0.010225420414434987
│   upper_bound = 3.8032318783602004
│   SA_time = 0.665278817
│   solve_time = 0.065247527
└   total_time_so_far = 8.786298954000001
┌ Info: Adding constraint 14 
│   violation = -0.00854797035176246
│   upper_bound = 3.8007075536226553
│   SA_time = 0.681930152
│   solve_time = 0.052101883
└   total_time_so_far = 9.468229106
┌ Info: Adding constraint 15 
│   violation = -0.020972059066668457
│   upper_bound = 3.795250646358601
│   SA_time = 1.36211752
│   solve_time = 0.072227499
└   total_time_so_far = 10.830346626
┌ Info: Adding constraint 16 
│   violation = -0.005417103821982476
│   upper_bound = 3.7934412873668957
│   SA_time = 1.316569325
│   solve_time = 0.058500157
└   total_time_so_far = 12.146915951
┌ Info: Adding constraint 17 
│   violation = -0.005734892762903272
│   upper_bound = 3.791735647425536
│   SA_time = 0.651098197
│   solve_time = 0.061050447
└   total_time_so_far = 12.798014148
┌ Info: Adding constraint 18 
│   violation = -0.0063616445714986085
│   upper_bound = 3.7878883216192785
│   SA_time = 0.679417382
│   solve_time = 0.092434717
└   total_time_so_far = 13.47743153
┌ Info: Adding constraint 19 
│   violation = -0.005431778550949503
│   upper_bound = 3.7877448070168134
│   SA_time = 0.681983077
│   solve_time = 0.084447351
└   total_time_so_far = 14.159414607
┌ Info: Adding constraint 20 
│   violation = -0.0017854551482699562
│   upper_bound = 3.7872340576765624
│   SA_time = 0.683947872
│   solve_time = 0.141246216
└   total_time_so_far = 14.843362479
┌ Info: Adding constraint 21 
│   violation = -0.0019491406748739794
│   upper_bound = 3.786289756311934
│   SA_time = 0.717835535
│   solve_time = 0.141211417
└   total_time_so_far = 15.561198014
┌ Info: Adding constraint 22 
│   violation = -0.004416909999571353
│   upper_bound = 3.7845998292930707
│   SA_time = 1.334771398
│   solve_time = 0.574838534
└   total_time_so_far = 16.895969412
┌ Info: Adding constraint 23 
│   violation = -0.0015878002498591042
│   upper_bound = 3.7835397473893293
│   SA_time = 3.334307019
│   solve_time = 0.371798179
└   total_time_so_far = 20.230276431
┌ Info: Adding constraint 24 
│   violation = -0.0009149900494965233
│   upper_bound = 3.783400416165109
│   SA_time = 1.349366087
│   solve_time = 0.336789226
└   total_time_so_far = 21.579642518
┌ Info: Adding constraint 25 
│   violation = -0.0006151044229729564
│   upper_bound = 3.782909203680918
│   SA_time = 2.080223126
│   solve_time = 0.336263933
└   total_time_so_far = 23.659865644
┌ Info: Adding constraint 26 
│   violation = -0.001641335103922728
│   upper_bound = 3.7825483496186707
│   SA_time = 1.356265292
│   solve_time = 0.491323406
└   total_time_so_far = 25.016130936
┌ Info: Adding constraint 27 
│   violation = -0.008160637987595719
│   upper_bound = 3.7810997747615143
│   SA_time = 2.726521266
│   solve_time = 0.542392309
└   total_time_so_far = 27.742652202
┌ Info: Next iteration is projected to exceed maximum time 30; terminating here. Change `max_time` to adjust this behavior.
│   total_time_so_far = 27.742652202
└   projected_time = 31.338457134499997