Cumulative update.

- Various documentation changes for function descriptions.

- decoyAnalysisIndependentLP:
  - forceSep now defaults to True. This will help ensure bounds don't cross over each other.
  - Minor update to sdp formatting.
  - yields are now clipped to the range [0,1] to reduce effects from cvx precision.

- Frank wolfe 2 step solver:
  - An improved method for perturbing the  objective function was found that does not require a continuity bound.
  - mustFollowPerturbationTheorem was removed and various functions were updated for this analysis.
  - Appendix B of the software manual was updated to describe this new method.

- minor changes to the step 1 solver's cvx problem formatting.

- updated ind2subPlus to fix compatibility with Matlab version 2024b after ind2sub was updated.

Author: JBurniston

BasicProtocols/BasicBB84WCPDecoy/decoyAnalysisIndependentLP.m

@@ -1,29 +1,39 @@
 function [conExpL,conExpU] = decoyAnalysisIndependentLP(conExp,decoys,options)
 % decoyAnalysisIndependentLP: A decoy analysis module for assymptotic data.
 % Returns the upper and lower bounds on the expectations condition on
-% single photons being sent. This function assumes that the source uses a
-% weak coherent pulsed laser to generate signals.
+% single photons being sent and signal choice. This function assumes that
+% the source uses a weak coherent pulsed laser to generate signals.
 %
 % Input:
 % * conExp: A table of expectation values of Bob's measurement result
 %   (columns) conditioned on Alice's measurement result (rows) and decoy
 %   intensity used (pages deep / third array dimension). Each row must be
-%   non-negative and sum to 1. Values of conditioned on intensity.
+%   non-negative and sum to 1 (ie. values are conditioned on intensity and
+%   signal sent by Alice). The first page is the key generating intensity.
 % * decoys: vector of intensities used for each page in conExp. Ordered by
 %   the page number. The intensities must be non-negative. Warning: an
 %   intensity of 0 is rather unstabe in numerical decoy analysis. We
 %   recommend you try something slightly above.
+%
 % Name-value pair options:
 % * decoyTolerance (1e-12): A small extra tolerance term to loosen the
 %   linear constraints by when solving.
 % * decoySolver ("SDPT3"): The solver you want CVX to use.
 % * decoyPrecision ("high"): The precision you want CVX to use.
 % * photonCutOff (10): The photon number cut off you want the decoy
 %   analysis to solve up to. 
-% * forceSep (false): For each expectation value, the single photon lower
+% * forceSep (true): For each expectation value, the single photon lower
 %   bound shouldn't be higher than the upper bound. Sometimes CVX does not
 %   satisfy this trivial bound. If set to true, an extra constraint will be
 %   added to force the single photon components to respect the bound.
+%
+% Output:
+% * conExpL: Lower bound on the probability of Bob's measurment outcomes
+%   (rows) conditioned on Alice's signal choice (columns) and single photon
+%   sent.
+% * conExpU: Upper bound on the probability of Bob's measurment outcomes
+%   (rows) conditioned on Alice's signal choice (columns) and single photon
+%   sent.
 %   
 %   Reviewed by Devashish Tupkary 2023/09/23
 arguments
@@ -33,7 +43,7 @@
     options.decoySolver (1,1) string = "SDPT3";
     options.decoyPrecision (1,1) string = "high";
     options.photonCutOff (1,1) double {mustBeInteger,mustBePositive} = 10;
-    options.forceSep (1,1) logical = false;
+    options.forceSep (1,1) logical = true;
 end
 
 
@@ -114,17 +124,18 @@
     cvx_solver(convertStringsToChars(decoySolver));
     cvx_precision(convertStringsToChars(decoyPrecision));
 
-    variables yieldsL(numPhotons,numOutcomes) yieldsU(numPhotons,numOutcomes)
+    yieldsL = variable('yieldsL(numPhotons,numOutcomes)','nonnegative');
+    yieldsU = variable('yieldsU(numPhotons,numOutcomes)','nonnegative');
 
     minimize(sum(yieldsL(2,:) - yieldsU(2,:))) % 2 corresponds to single-photon yields
     subject to
 
-    0 <= yieldsL <= 1;
+    yieldsL <= 1;
 
     poissonCache*yieldsL <= conExp                +decoyTolerance;
     poissonCache*yieldsL >= conExp -probRemaining -decoyTolerance;
 
-    0 <= yieldsU <= 1;
+    yieldsU <= 1;
 
     poissonCache*yieldsU <= conExp                +decoyTolerance;
     poissonCache*yieldsU >= conExp -probRemaining -decoyTolerance;
@@ -141,8 +152,10 @@ variables yieldsL(numPhotons,numOutcomes) yieldsU(numPhotons,numOutcomes)
     throwAsCaller(exception);
 end
 
-conExpLowerBound = yieldsL(2,:);
-conExpUpperBound = yieldsU(2,:);
+%read out single photons from the yields and clip between 0 and 1 to remove
+%numerical errors.
+conExpLowerBound = min(max(yieldsL(2,:),0),1);
+conExpUpperBound = min(max(yieldsU(2,:),0),1);
 
 end
 

MathSolverModules/OQCTFrankWolfe2StepSolver/FW2StepSolver.m

@@ -68,7 +68,9 @@
 %   more details (the second argument, x1, in the function).
 % * linearConstraintTolerance (1e-10): constraint tolerance on the
 %   equalityConstraints, inequalityConstraints, vectorOneNormConstraints
-%   and matrixOneNormConstraints. A bit broader than the name suggests.
+%   and matrixOneNormConstraints for step 1 of the solver. This is to help
+%   the solver find feasible points during step 1 and play no role in step
+%   2.
 % * initMethod (1): Integer selected from {1,2,3}. For the Frank Wolfe
 %   algorithm, the initial point must satisfy the constraints. This selects
 %   which technique is used, from 3 choices:

MathSolverModules/OQCTFrankWolfe2StepSolver/HelperFunctions/mustFollowPerturbationTheorem.m

@@ -42,13 +42,10 @@
     end
 
     epsilon = dim*(minEigenvalue-eigMin)/(real(trace(rho))-eigMin*dim);
-    try
-        %check 
-        if options.perturbationCheck
-            mustFollowPerturbationTheorem(epsilon,rho);
-        end
-    catch ME
-        rethrow(ME)
+    %check
+    if options.perturbationCheck && (epsilon < 0 || epsilon > 1)
+        throw(MException("perturbationCHannelEpsilon:UnphysicalPerturbation", ...
+            "The perturbation parameter calculated was not in the range [0,1]."))
     end
 end
 

MathSolverModules/OQCTFrankWolfe2StepSolver/HelperFunctions/perturbationChannelEpsilon.m

@@ -15,12 +15,12 @@
 % * Dfval: The value of the gradient at the given value of rho. This uses
 %   numerator convention.
 %
-% See also primalDf, primalfep, FW2StepSolver, perturbationChannelEpsilon
+% See also primalDf, primalfep, FW2StepSolver
 arguments
     %minimial checks just to make sure cells are formatted in the correct
     %orientation.
-    perturbation (1,1) double
-    rho (:,:) double {mustBeHermitian,mustFollowPerturbationTheorem(perturbation,rho)}
+    perturbation (1,1) double {mustBeInRange(perturbation,0,1)}
+    rho (:,:) double {mustBeHermitian}
     keyProj (:,1) cell %checks are too complex for this.
     krausOperators (:,1) cell %checks are too complex for this.
     safeCutOff (1,1) double {mustBePositive} = 1e-14;

MathSolverModules/OQCTFrankWolfe2StepSolver/HelperFunctions/primalDfep.m

@@ -17,12 +17,12 @@
 % * realEpsilon: The epsilon value that was used to compute fval. Must be
 %   between 0 and 1.
 %
-% See also primalf, primalDfep, FW2StepSolver, perturbationChannelEpsilon
+% See also primalf, primalDfep, FW2StepSolver
 arguments
     %minimial checks just to make sure cells are formatted in the correct
     %orientation.
-    perturbation (1,1) double
-    rho (:,:) double {mustBeHermitian, mustFollowPerturbationTheorem(perturbation,rho)}
+    perturbation (1,1) double {mustBeInRange(perturbation,0,1)}
+    rho (:,:) double {mustBeHermitian}
     keyProj (:,1) cell
     krausOperators (:,1) cell
     safeCutOff (1,1) double {mustBePositive} = 1e-14;

MathSolverModules/OQCTFrankWolfe2StepSolver/HelperFunctions/primalfep.m

@@ -157,7 +157,8 @@
         variable(convertStringsToChars(rhoStrings(index)),'hermitian','semidefinite')
     end
     % Get a single cell array to capture those variables.
-    rhoBlocks = eval("{"+strjoin(compose("rho%d",1:numel(newDims)),",")+"}");
+    rhoBlocksString = "{"+strjoin(compose("rho%d",1:numel(newDims)),",")+"}";
+    rhoBlocks = eval(rhoBlocksString);
 
 
     % there is a bug in cvx's blkdiag override that messes up some times.
@@ -192,6 +193,7 @@ minimize norm(rho0-rho)
 % reconstruct rho from blocks incase it didn't get filled in by CVX
 if options.blockDiagonal
     rho = zeros(dim);
+    rhoBlocks = eval(rhoBlocksString);
     rho = directSumWorkArround(rho, rhoBlocks);
     rho = P'*rho*P;
 end
@@ -227,7 +229,8 @@ minimize norm(rho0-rho)
         variable(convertStringsToChars(deltaRhoStrings(index)),'hermitian') %NOT semidefinite
     end
     % Get a single cell array to capture those variables.
-    deltaRhoBlocks = eval("{"+strjoin(compose("deltaRho%d",1:numel(newDims)),",")+"}");
+    deltaRhoBlocksString = "{"+strjoin(compose("deltaRho%d",1:numel(newDims)),",")+"}";
+    deltaRhoBlocks = eval(deltaRhoBlocksString);
 
 
     % there is a bug in cvx's blkdiag override that messes up some times.
@@ -255,6 +258,7 @@ minimize real(trace(gradf*deltaRho)) % numerator form
 % reconstruct rho from blocks incase it didn't get filled in by CVX
 if options.blockDiagonal
     deltaRho = zeros(dim);
+    deltaRhoBlocks = eval(deltaRhoBlocksString);
     deltaRho = directSumWorkArround(deltaRho, deltaRhoBlocks);
     deltaRho = P'*deltaRho*P;
 end

MathSolverModules/OQCTFrankWolfe2StepSolver/HelperFunctions/step1Solver.m

@@ -20,7 +20,7 @@
 [perturbation,safeCutOff] = computePerturbationEpsilon(rho, krausOps, keyProj);
 debugInfo.storeInfo("perturbationValue",perturbation);
 
-% 2. Calculate f_epsilon_p(rho) and gradftranspose of that
+% 2. Calculate f_epsilon_p(rho) and its gradient at the same point.
 fval = primalfep(perturbation, rho, keyProj, krausOps,safeCutOff);
 debugInfo.storeInfo("relEntStep2Linearization",fval/log(2));
 
@@ -37,16 +37,12 @@
 % compute key rate lower bound from the result of step 2
 beta = dualSolution + fval - real(trace(gradf*rho)); %real() deals with small numerical instability
 
-% We compute the penalty zetaEp due to pertubation
-if perturbation == 0
-    zetaEp = 0; %limit as perturbation -> 0
-else
-    dimOut = size(krausOps{1},1);
-    zetaEp = 2 * 2*perturbation*(dimOut-1)/dimOut...
-        *log(dimOut^2/(2*perturbation*(dimOut-1)));
-end
+% Because f(maxMixed) = 0, we can use convex arguments to get
+% f_\epsilon(\rho)/(1-\epsilon) \leq f(\rho)$. This is way cleaner than
+% other perturbation bounds
+
 
-lowerBound = beta - zetaEp;
+lowerBound = beta/(1-perturbation);
 end
 
 %% Other Functions

MathSolverModules/OQCTFrankWolfe2StepSolver/HelperFunctions/step2Solver.m

@@ -39,7 +39,7 @@ Our software requires *at least version 2020b* for full functionality, but insta
 Our software has the following dependencies for its default settings:
 
 - [CVX](https://cvxr.com/cvx/download/) v2.2, a library for solving convex optimization problems in MATLAB.
-- [QETLAB](https://github.com/nathanieljohnston/QETLAB) *above* v0.9, a MATLAB toolbox for operations involving quantum channels and entanglement. Note that you cannot use the version from their website as it has a bugs associated with implementing Choi matrices. *You must use their latest copy on GitHub*. At the time of writing, this means downloading their code with the green "code" button and *not* the v0.9 release.
+- [QETLAB](https://github.com/nathanieljohnston/QETLAB) *above* v0.9, a MATLAB toolbox for operations involving quantum channels and entanglement. Note that you cannot use the version from their website as it has a bugs associated with implementing Choi matrices. *You must use their latest copy on Github*. At the time of writing, this means downloading their code with the green "code" button and *not* the v0.9 release.
 - [ZGNQKD](https://www.math.uwaterloo.ca/~hwolkowi/henry/reports/ZGNQKDmainsolverUSEDforPUBLCNJuly31/) solver (optional), an alternative to our Frank-Wolfe solver for quantum relative entropy. Currently, it only supports equality constraints. Download the zip and the "Solver" folder and sub-folders to your path.
 - [MOSEK](https://www.mosek.com/) (optional), a more advanced semidefinite programming (SDP) solver than the default (SDPT3) used in CVX. Note that the MOSEK solver can be downloaded together with CVX, but requires a separate license to use. See [this page](https://cvxr.com/cvx/doc/mosek.html) for more information.