QDimSum

Logo

Symmetric SDP relaxations for qudits systems

View the Project on GitHub denisrosset/qdimsum

QDimSum

This package was written by Denis Rosset, Armin Tavakoli and Marc-Olivier Renou.

It implements the algorithms described in

and is based on the Navascués-Vértesi hierarchy described in

We also mention the first use of symmetrisation applied to the Navascués-Vértesi hierarchy in the independent work:

The project is licensed under the 3-Clause BSD License.

If you use our software to produce research outcomes, please cite the papers mentioned above.

Installation

  1. Verify that you have a recent Matlab version (any version >= 2015b should work).

  2. Download and install YALMIP.

  3. Download and install the SDP solver of your choice. We recommend MOSEK (free for academia), or SeDuMi (open source).

  4. Run yalmiptest and verify that SDP problems can be solved in your configuration.

  5. Download the software package QDimSum from this URL.

  6. Unpack the qdimsum-master.zip file in a folder of your choice.

  7. Add your chosen folder location (/qdimsum-master) to your Matlab path (only that folder and not its subfolders).

  8. Run the script TestRAC22.m present in the qdimsum-master folder. The script should finish without any error.

You are now ready to start learning about QDimSum.

Simplest example (no symmetrisation)

We now consider the RAC example of our paper for $n=2$ and $d=2$.

For that, we create a file named RAC22a.m, which is a class that extends the NVProblem class provided by QDimSum.

classdef RAC22a < NVProblem
    properties
       forceReal = true;
    end
    methods
        function X = sampleOperators(self, rank)
        % the first 4 operators are the states, next 2 a projective measurement
        % the last 2 another projective measurement
            dim = 2;
            X = cell(1, 8);
            for i = 1:4
                X{i} = ...
  qdimsum.Random.pureNormalizedDensityMatrix(dim);
            end
            U = qdimsum.Random.unitary(2);
            X{5} = U*[1 0; 0 0]*U';
            X{6} = U*[0 0; 0 1]*U';
            U = qdimsum.Random.unitary(2);
            X{7} = U*[1 0; 0 0]*U';
            X{8} = U*[0 0; 0 1]*U';
        end
        function K = sampleStateKraus(self)
            K = eye(2); % dimension is 2
        end
        function obj = computeObjective(self, X, K)
            obj = 0;
            for x1 = 1:2
                for x2 = 1:2
                    for y = 1:2
                        if y == 1
                            b = x1;
                        else
                            b = x2;
                        end
                        rho = X{x1+(x2-1)*2};
                        M = X{4+b+(y-1)*2};
                        obj = obj + real(trace(M * rho)/8);
                    end
                end
            end
        end
    end
end

Now, a few explanations.

Now, to optimize over that problem, we write the following in a TestRAC22a.m script:

settings = NVSettings;
problem = RAC22a;
monomials = {'npa' 2};
disp('The objective value is:')
upperBoundSDP = nvOptimize(problem, monomials, ...
  'none', settings)
disp('While the analytical value is:')
sol = 1/2*(1+1/sqrt(2))

Again, a few explanations:

Using symmetrisation

To use symmetrisation, you need to specify the symmetry group of the problem. For that, add the method below to the class RAC22 (just before or after computeObjective). For our purposes, we renamed the class RAC22b for that second step.

        function generators = symmetryGroupGenerators(self)
            swapX1X2 = [1 3 2 4 7 8 5 6];
            flipX1 = [2 1 4 3 6 5 7 8];
            generators = [swapX1X2
                          flipX1];
        end

The symmetryGroupGenerators method specifies the symmetries one wishes to exploit. A symmetry corresponds to a permutation of the sequence of physical operators which leaves the objective function invariant. If there are $N$ physical operators in the problem, a symmetry group element is specified by providing the image of the list $(1, \ldots N)$ under that symmetry group element. If the symmetry group has $M$ generators, the method symmetryGroupGenerators should return a $M \times N$ matrix containing one such permutation per row.

Here, The generator swapX1X2 permutes the indices $x_1$ and $x_2$ in $\rho_{x1,x2}$, while flipX1 permutes $x_1$.

Now, you have access to all the symmetrisation methods of QDimSum, we use the following script:

settings = NVSettings;
problem = RAC22a;
monomials = {'npa' 2};
nvOptimize(problem, monomials, 'none', settings)
nvOptimize(problem, monomials, 'reynolds', settings)
nvOptimize(problem, monomials, 'isotypic', settings)
nvOptimize(problem, monomials, 'irreps', settings)
nvOptimize(problem, monomials, 'blocks', settings)

The application of summarization to a problem amounts to two parts. First we reduce the number of linearly independent basis elements found through the sampling procedure. This number will decrease as the number of applied independent symmetries increases. Second, the final SDP matrix can be block-diagonalized. The size and number of these blocks depends on the symmetry group identified.

The software package offers several methods for symmetrized implementations of the NV hierarchy.

Optionally, if the isotypic or irreducible decomposition of the group action is already known, it can be passed to the nvOptimize method (see its documentation).

Each of these methods offers a different degree of symmetrization which may be more or less relevant depending on the specific problem of interest. In particular, due to a temporary limitation of our code, irreps and blocks require all the irreducible representations appearing to be of real type.

Note: we allow permutations to flip the sign of the operator. For that, just flip the sign of the integer in the image (see I3322.m for an example).

Controlling monomials: families of operators

Instead of the NPA type hierarchy, one may want to use local levels to decide on the monomials involved.

For that purpose, one has to provide an additional method operatorTypes:

        function types = operatorTypes(self)
            types = {1:4 5:8};
        end

where types is a cell array marking each physical operator (states, measurements etc.) that appear in the physical problem of interest. For example, in a prepare and measure scenario with N preparations and M measurement operators, the cell array reads {1:N N+(1:M)}, i.e., the first N positions will correspond to preparations and the final M to the measurement operators. Similarly, in a bipartite Bell experiment with N measurement operators for each party, the cell array reads {1:N N+(1:N)}.

Now, one can specify that the monomials involved are $1$, $\rho$, $M$, $\rho \cdot M$ by running the code as follows:

settings = NVSettings;
problem = RAC22c;
monomials = {'families' [] [1] [2] [1 2]};
nvOptimize(problem, monomials, 'blocks', settings)

See the files RAC22c.m and TestRAC22c.m.

Controlling settings

We describe below the most important options in NVSettings. To change them from the default options, one calls NVSettings as such:

settings = NVSettings( ...
  'yalmipSettings', ...
  sdpsettings('solver', 'sedumi', 'verbose', 0), ...
  'verbosityLevel', 0);

This selects SeDuMi as a solver (by passing the relevant options to YALMIP), and prevents the code to output its progress, by silencing both QDimSum (with verbosityLevel) and SeDuMi (verbose = 0).

The documentation of all options can be found in the NVSettings.m source file.

Sanity checks

In the options above, the setting checkLevel is set to 1 by default. Then, QDimSum will perform consistency checks at every step of the computation. Those checks are not expensive and should not be disabled without good reason. Moreover, one can help QDimSum in checking the user provided data by explaining to the toolbox the constraints that the operators have to satisfy.

        function C = operatorSDPConstraints(self, X)
            C = X; % all operators are SDP
        end
        function C = operatorEqualityConstraints(self, X)
            dim = 2;
            % X{5}, X{6} form a projective measurement
            % X{7}, X{8} form a projective measurement
            C = {eye(dim) - X{5} - X{6}
                 eye(dim) - X{7} - X{8}};
        end
        function C = scalarEqualityConstraints(self, X)
            dim = 2;
            mm = eye(dim)/dim; % states have the same
			% trace as the maximally mixed state
            C = {mm - X{1}
                 mm - X{2}
                 mm - X{3}
                 mm - X{4}};
        end

The documentation of those constraints is given in NVProblem.m, here they correspond respectively to: the positive semidefiniteness of the operators, to the projector elements summing to identity, and the states being normalized.

In the file RAC22d.m, we introduced an error on purpose in the symmetry group generators. When we run TestRAC22d.m, we get the error message:

Error using qdimsum.Check.sampleObeysConstraints (line 68)
For generator 1  3  2  4  7  6  5  8: operator
equality constraint #1 violated, maximal
singular value 0.960222

that identifies the problem: run TestRAC22d.m to reproduce that error and its detection.

Additional tools