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.0imExample 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.5000105244940878We 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.9999976286547437We 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.7094306690270955plot_pmfN(output)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.7499750650637333We 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)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.776910138848433A 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.04680680979341
│ upper_bound = 8.500000024707353
│ SA_time = 0.849038939
│ solve_time = 0.003512435
└ total_time_so_far = 0.849038939
┌ Info: Adding constraint 2
│ violation = -2.0523990104184393
│ upper_bound = 6.220440713086008
│ SA_time = 0.886167186
│ solve_time = 0.004042241
└ total_time_so_far = 1.735206125
┌ Info: Adding constraint 3
│ violation = -1.1779629749393707
│ upper_bound = 5.045704670043538
│ SA_time = 0.853867632
│ solve_time = 0.005380355
└ total_time_so_far = 2.589073757
┌ Info: Adding constraint 4
│ violation = -1.2060853680043822
│ upper_bound = 3.839619058361676
│ SA_time = 0.911973793
│ solve_time = 0.006682667
└ total_time_so_far = 3.50104755
┌ Info: Adding constraint 5
│ violation = -0.06437541958623592
│ upper_bound = 3.8324170077570567
│ SA_time = 0.826765412
│ solve_time = 0.009661397
└ total_time_so_far = 4.327812962
┌ Info: Adding constraint 6
│ violation = -0.04000065884375521
│ upper_bound = 3.8287756984457646
│ SA_time = 0.853450958
│ solve_time = 0.011282013
└ total_time_so_far = 5.18126392
┌ Info: Adding constraint 7
│ violation = -0.0591901262392904
│ upper_bound = 3.8190928671224924
│ SA_time = 0.859450695
│ solve_time = 0.013190332
└ total_time_so_far = 6.040714615000001
┌ Info: Adding constraint 8
│ violation = -0.05996059564871384
│ upper_bound = 3.814536396317005
│ SA_time = 0.883892516
│ solve_time = 0.014979349
└ total_time_so_far = 6.924607131
┌ Info: Adding constraint 9
│ violation = -0.026520451319771314
│ upper_bound = 3.8117409221714427
│ SA_time = 0.866112314
│ solve_time = 0.017572674
└ total_time_so_far = 7.790719445000001
┌ Info: Adding constraint 10
│ violation = -0.02486525932573949
│ upper_bound = 3.8101190516016934
│ SA_time = 0.913927865
│ solve_time = 0.021904017
└ total_time_so_far = 8.70464731
┌ Info: Adding constraint 11
│ violation = -0.016884146017112127
│ upper_bound = 3.8033329724683047
│ SA_time = 0.901609317
│ solve_time = 0.033197027
└ total_time_so_far = 9.606256627
┌ Info: Adding constraint 12
│ violation = -0.043674004455128396
│ upper_bound = 3.795218912712105
│ SA_time = 1.772686255
│ solve_time = 0.034401838
└ total_time_so_far = 11.378942882
┌ Info: Adding constraint 13
│ violation = -0.01064128877882862
│ upper_bound = 3.7949078158003617
│ SA_time = 1.739156132
│ solve_time = 0.038358375
└ total_time_so_far = 13.118099014
┌ Info: Adding constraint 14
│ violation = -0.01188960435198699
│ upper_bound = 3.7849921814142764
│ SA_time = 0.910068576
│ solve_time = 0.042626415
└ total_time_so_far = 14.02816759
┌ Info: Adding constraint 15
│ violation = -0.002195968730972395
│ upper_bound = 3.7849734630404948
│ SA_time = 0.912103971
│ solve_time = 0.062756709
└ total_time_so_far = 14.940271561000001
┌ Info: Adding constraint 16
│ violation = -0.006707546957262572
│ upper_bound = 3.784522372223587
│ SA_time = 1.780438145
│ solve_time = 0.053103013
└ total_time_so_far = 16.720709706
┌ Info: Adding constraint 17
│ violation = -0.0056704150754510255
│ upper_bound = 3.7831459543223698
│ SA_time = 0.883361521
│ solve_time = 0.065860034
└ total_time_so_far = 17.604071227000002
┌ Info: Adding constraint 18
│ violation = -0.004352621904843135
│ upper_bound = 3.7826680958032544
│ SA_time = 1.691506454
│ solve_time = 0.528323456
└ total_time_so_far = 19.295577681
┌ Info: Adding constraint 19
│ violation = -0.004593414079917464
│ upper_bound = 3.7804753991598643
│ SA_time = 2.57377804
│ solve_time = 0.090630259
└ total_time_so_far = 21.869355721
┌ Info: Adding constraint 20
│ violation = -0.006562062208066941
│ upper_bound = 3.7803383491328706
│ SA_time = 1.851010132
│ solve_time = 0.097861823
└ total_time_so_far = 23.720365853
┌ Info: Adding constraint 21
│ violation = -0.0003872766418948155
│ upper_bound = 3.780346212652025
│ SA_time = 1.841163743
│ solve_time = 0.099216331
└ total_time_so_far = 25.561529596
┌ Info: Adding constraint 22
│ violation = -0.00038883909794259586
│ upper_bound = 3.780330642925389
│ SA_time = 2.645417752
│ solve_time = 0.115292876
└ total_time_so_far = 28.206947348
┌ Info: Next iteration is projected to exceed maximum time 30; terminating here. Change `max_time` to adjust this behavior.
│ total_time_so_far = 28.206947348
└ projected_time = 31.2437290388