Hello!
Introduction
I am a PhD student in quantum information, and I am currently studying measurement-device independent quantum random number generation. Unfortunately I am not familiar with SDP, as I have only worked with linear programming in my previous studies, and I am somewhat stuck in a situation where I don’t know whether it is my problem formulation, or the way I am inputting it into CVX’s solver that is incorrect.
For random generation, a superposition state is prepared and measured, which is expected to yield unbiased, random measurement outcomes (either binary 1 or binary 0 with equal probability). In order to estimate the amount of information that could have leaked to an adversary, one prepares and measures states which are expected to yield deterministic outcome (depending on choice of test state, either 1 or 0 every time). These states are typically called test sates. If an attacker (Eve) is present, preparing and measuring a test state does not yield the expected outcomes, and by experimentally determining the probabilities Pr{a|w_x} where a is the measurement outcome (1/0) and w_x is the prepared state, one is able to bound the probability that an attacker would be able to guess the measurement outcome and obtain information about the generated sequence.
Problem formulation
I am trying to replicate the security analysis by Cariñe et al [Carñe et al. Optica 7, 542-550 (2020)], which looks like this:
I have been able to formulate this problem in 2D for one test state and one state for randomness generation, where the former corresponds to a horizontal polarization state, and the latter diagonal polarization. The optimization is over a POVM which is represented by four matrices, one for each possible outcome for Alice and Eve, where the POVM encodes all possible guessing strategies Eve could possibly choose.
For now, experimentally-derived values are simulated so that there is a high probability to measure correctly.
My issue
While I believe that I have deciphered the mathematical formulation of the general case and formulated it into a form applicable for my case, I fail to input the problem into the CVX solver. My problem seems to be unbounded, and I cannot see why that is the case, as I have numerical values for all constraints.
a, e are measurement outcomes.
The question
- Have I entered the problem properly, or is there anything immediate that I have missed?
- I do not know the value of q(e), would that be necessary for converegence?
Sincerely,
Joakim
Appendix: MATLAB-code
%% Compute parameters
clc; clear all;
epsilon = 1e-9; % by convention
det_H = 1000;
det_V = 5;
prob_correct_measurement = (det_H)/(det_H+det_V);
%% Input state
x_star = QNorm([1;0]);
omega_x_star = QKet(x_star)*QBra(x_star);
disp("Input state density matrix: ")
disp(omega_x_star)
x_0 = QNorm([1;0]);
x_1 = QNorm([1;1]);
w_0 = QKet(x_0)*QBra(x_0);
w_1 = QKet(x_1)*QBra(x_1);
%% Observed probabilities + Chernoff-Hoeffding bound
p_corr = prob_correct_measurement;
p_ones = 0.49;
% states = [w_0, w_1] w_0: teststate, w_1: random state
pr_obs = [1-p_corr, 1-p_ones
p_corr, p_ones];
disp("Observed probabilities")
disp(pr_obs)
disp(" |w0> |w1>")
disp("0: " + pr_obs(1,1) + " " + pr_obs(1,2))
disp("1: " + pr_obs(2,1) + " " + pr_obs(2,2))
t_x_helper = @(n_x, e) sqrt( (log(1/e))/(2*n_x));
n_tot = 1000;
p_teststate = 0.01;
n_w0 = p_teststate*n_tot;
n_w1 = (1-p_teststate)*n_tot;
t_x = [t_x_helper(n_w0, epsilon), t_x_helper(n_w1, epsilon)];
disp(" ")
disp("State 'slack'")
disp("w_0: " + t_x(1))
disp("w_1: " + t_x(2))
II = eye(2,2);
%% Print LHS & RHS of constrs
disp(" ")
disp("Constraint matrix")
for i=1:2
for j=1:2
disp("pr_obs(" + i +","+j +") - t_x(1): " + (pr_obs(i, j) - t_x(j)))
disp("pr_obs(" + i +","+j +") + t_x(1): " + (pr_obs(i, j) + t_x(j)))
disp(" ")
end
end
%% Solver
cvx_clear
% N_xy denotes the POVM elements
% q_x encodes the no-signalling constraint
cvx_begin SDP
variable N_00(2,2) hermitian
variable N_01(2,2) hermitian
variable N_10(2,2) hermitian
variable N_11(2,2) hermitian
variable q_0
variable q_1
maximize(trace(N_00*omega_x_star + N_11 * omega_x_star))
subject to
pr_obs(1, 1) - t_x(1) <= trace(N_00*w_0 + N_01*w_0) <= pr_obs(1,1) + t_x(1); % Prob. measure '0' when w_0 sent
pr_obs(1, 2) - t_x(2) <= trace(N_00*w_1 + N_01*w_1) <= pr_obs(1,2) + t_x(2); % Prob. measure '0' when w_1 sent
pr_obs(2, 1) - t_x(1) <= trace(N_10*w_0 + N_11*w_0) <= pr_obs(1,1) + t_x(1); % Prob. measure '1' when w_0 sent
pr_obs(2, 2) - t_x(2) <= trace(N_10*w_1 + N_11*w_1) <= pr_obs(1,2) + t_x(2); % Prob. measure '1' when w_1 sent
N_00 + N_10 == q_0 * II;
N_01 + N_11 == q_1 * II;
q_0 + q_1 == 1;
cvx_end
p_g_x_star = cvx_optval;
cert_entropy = -log2(p_g_x_star);
disp("Certified entropy: " + cert_entropy + " bits")
The helper functions QNorm, QKet and QBra are simply corresponding to normalization of a quantum state, and formatting of a quantum state vector to a ket and bra.
QNorm.m
function [outp] = QNorm(inp)
outp = (1/norm(inp)) * inp;
end
QKet.m
function [ket] = QKet(state)
ket = state;
end
QBra.m
function [bra] = QBra(ket)
bra = ket'; % Hermitian conjugate
end