{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Advanced\n", "\n", "In this section we will consider some advanced aspects related to the package." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solving the dual versus the primal formulation of the SDP\n", "\n", "For semidefinite programs that appear often in causal compatibility problems, using the dual formulation speeds up the solve time and significantly lowers RAM usage.\n", "\n", "Consider the following example, where we use the MOSEK Fusion API to solve the primal version of a program, and then the dual:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SDP solver took 3.70 seconds\n", "The primal formulation was solved in 3.823288917541504 seconds.\n", "SDP solver took 0.46 seconds\n", "The dual formulation was solved in 0.5328986644744873 seconds.\n" ] } ], "source": [ "from inflation import InflationProblem, InflationSDP\n", "from time import time\n", "import numpy as np\n", "\n", "qtriangle = InflationProblem(dag={\"rho_AB\": [\"A\", \"B\"],\n", " \"rho_BC\": [\"B\", \"C\"],\n", " \"rho_AC\": [\"A\", \"C\"]}, \n", " outcomes_per_party=[2, 2, 2],\n", " settings_per_party=[1, 1, 1],\n", " inflation_level_per_source=[2, 2, 2])\n", "sdprelax = InflationSDP(qtriangle, verbose=0)\n", "sdprelax.generate_relaxation('npa2')\n", "\n", "P_W = np.zeros((2, 2, 2, 1, 1, 1))\n", "for a, b, c in np.ndindex((2, 2, 2)):\n", " if a + b + c == 1:\n", " P_W[a, b, c, 0, 0, 0] = 1 / 3\n", "\n", "sdprelax.set_distribution(P_W)\n", "\n", "time0 = time()\n", "sdprelax.solve(dualise=False)\n", "print(\"The primal formulation was solved in\", time() - time0, \"seconds.\")\n", "time0 = time()\n", "sdprelax.solve(dualise=True)\n", "print(\"The dual formulation was solved in\", time() - time0, \"seconds.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that there is an order of magnitude difference between the primal and dual formulations of the same problem. This is not true for all problems, but for the semidefinite programming relaxations generated for causal compatibility, almost always the dual formulation is more efficient. This should be taken into account when attempting to solve a relaxation. In what follows, we recompile some useful information for different interfaces.\n", "\n", "- [CVXPY](https://www.cvxpy.org/). If you export the problem to CVXPY, the behaviour depends on the solver you choose to use. When choosing MOSEK, note that CVXPY [dualises by default](https://www.cvxpy.org/tutorial/advanced/index.html?highlight=dualization) all continuous problems. There is [no automatic dualisation option](https://github.com/cvxpy/cvxpy/issues/1403). There is no option to specify whether to solve the primal or dual problem. Thus if you wanted to solve the primal with MOSEK, you would need to write the dual formulation manually, which when dualised would solve the primal (it is not expected that the user will need to do this!).\n", "- [PICOS 2.4](https://picos-api.gitlab.io/picos/). Picos [supports dualisation](https://picos-api.gitlab.io/picos/api/picos.modeling.options.html#option-dualize) with the `dualise=True` options flag. See [this issue](https://gitlab.com/picos-api/picos/-/issues/280) for more details. \n", "- [YALMIP](https://yalmip.github.io/). Like CVXPY, YALMIP [automatically dualises](https://yalmip.github.io/tutorial/automaticdualization) problems, however there is a flag, `dualize`, in `sdpsettings` to disable this feature if so desired.\n", "- MOSEK Fusion API. Our implementation of the semidefinite programming relaxation supports both the primal and dual formulations, as seen in the example above. This is done manually, as MOSEK Fusion API does not have functionality to change from the primal to the dual formulations.\n", "\n", "\n", "## Large scale problems\n", "\n", "For solving large scale semidefinite programs, it is recommended to use the MOSEK Fusion API, as going through interfaces for conic problems, such as PICOS or CVXPY, usually has an overhead in the pre-processing state (for example, there can be a higher RAM usage in the preprocessing stage than when solving the problem, which can lead to out-of-memory errors). There does not seem to be such an overhead when using YALMIP. For using YALMIP, the user can export the problem to `.dat-s` format using `InflationSDP.write_to_file()`, and load it in MATLAB using YALMIP's `loadsdpafile`.\n", "\n", "For large problems, it is recommended to try using a first-order SDP solver, such as [SCS](https://www.cvxgrp.org/scs/), if using second-order SDP solvers, such as MOSEK, is too slow or too memory-consuming. To use SCS the problem needs to be exported to the user's interface of choice and have SCS installed.\n", "\n", "It is also worth considering using symmetries to block-diagonalise the semidefinite program. This can be done with [RepLAB](https://replab.github.io/web/) in MATLAB. Symmetries arising from inflation can be calculated with `InflationSDP._discover_inflation_symmetries()`, and they are encoded as permutations of the list of generating monomials which leave the SDP invariant. This then can be used in RepLAB to block-diagonalise the problem, such as in the following example.\n", "\n", "## Using RepLAB in MATLAB to block-diagonalise an inflation SDP\n", "\n", "In the following example, we generate the SDP relaxation for a problem, in this case, optimisation of the Mermin inequality over the triangle scenario with quantum sources, and show how to use RepLAB to block diagonalise the SDP and solve it in MATLAB.\n", "\n", "First, we generate the SDP relaxation and write it to file. We generate a particularly large SDP (moment matrix with around 1700 columns) to showcase the advantages of block diagonalisation." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of columns in the moment matrix: 1729\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Calculating moment matrix: 100%|██████████| 1495585/1495585 [04:03<00:00, 6131.07it/s]\n", "Applying symmetries : 100%|██████████| 215/215 [00:00<00:00, 770.54it/s]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Number of variables after symmetrization: 6476\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Initializing monomials : 100%|██████████| 6476/6476 [00:08<00:00, 762.65it/s] \n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Writing the SDP program to inflation_sdp.mat\n", "Writing the inflation symmetries to inflation_symmetries.mat\n" ] } ], "source": [ "from inflation import InflationProblem, InflationSDP\n", "from scipy.io import savemat\n", "\n", "scenario = InflationProblem(dag={\"rho_AB\": [\"A\", \"B\"],\n", " \"rho_BC\": [\"B\", \"C\"],\n", " \"rho_AC\": [\"A\", \"C\"]},\n", " outcomes_per_party=(2, 2, 2),\n", " settings_per_party=(2, 2, 2),\n", " inflation_level_per_source=(3, 3, 3))\n", "sdp = InflationSDP(scenario, commuting=False, verbose=1)\n", "sdp.generate_relaxation(\"npa2\")\n", "\n", "mmnts = sdp.measurements\n", "A0, B0, C0, A1, B1, C1 = (1 - 2*mmnts[party][0][setting][0] # As correlators\n", " for setting in range(2) for party in range(3))\n", "sdp.set_objective(A1*B0*C0 + A0*B1*C0 + A0*B0*C1 - A1*B1*C1)\n", "\n", "sdp.write_to_file('inflation_sdp.mat')\n", "print(\"Writing the inflation symmetries to inflation_symmetries.mat\")\n", "savemat('inflation_symmetries.mat',\n", " {'inflation_symmetries': sdp.inflation_symmetries + 1} # because matlab indexing starts at 1\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Within a MATLAB script, after installing [RepLAB](https://github.com/replab/replab), we need to initialise it by running `replab_init`:\n", "\n", "```MATLAB\n", "run \"X:\\...\\replab-develop\\replab_init.m\"\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we load the SDP and symmetries from file, and create a Yalmip `sdpvar` for the moment matrix:\n", "\n", "```MATLAB\n", "load('inflation_sdp.mat');\n", "load('inflation_symmetries.mat');\n", "\n", "% Convert symmetries to cell for compatibility with RepLAB functions\n", "symmetries = num2cell(double(inflation_symmetries), 2);\n", "\n", "% Build momentmatrix as a Yalmip sdpvar\n", "IndexMatrix = double(momentmatrix); % int32 -> double\n", "vars = sdpvar(1, max(IndexMatrix(:)));\n", "Gamma = vars(IndexMatrix);\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Computing the block structure with RepLAB is as easy as running:\n", "\n", "```MATLAB\n", "GammaBlocks = replab.CommutantVar.fromSymSdpMatrix(Gamma, symmetries);\n", "```\n", "\n", "which generates a block diagonal moment matrix (in about 1.5 minutes):\n", "```\n", ">> disp(GammaBlocks)\n", "Commutant variable 1729x1729 (23 blocks, 6476 scalar variables)\n", " dim: 1729\n", " dimensions1: [1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 4, 4, 4, 4, 4, 4, 8]\n", " field: 'real'\n", " matrixType: 'symmetric'\n", "multiplicities: [70, 3, 3, 3, 14, 14, 14, 12, 12, 12, 12, 12, 12, 68, 68, 68, 4, 4, 4, 60, 60, 60, 36]\n", " types: 'RRRRRRRRRRRRRRRRRRRRRRR'\n", " U: 1729 x 1729 double\n", " blocks: [70, 3, 3, 3, 14, 14, 14, 12, 12, 12, 12, 12, 12, 68, 68, 68, 4, 4, 4, 60, 60, 60, 36]\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we solve the SDP with Yalmip:\n", "\n", "```MATLAB \n", "constraints = [GammaBlocks(1,1) == 1, GammaBlocks >= 0];\n", "mermin = [vars(objective(:, 1))] * objective(:, 2);\n", "optimize(constraints, mermin);\n", "```\n", "\n", "While solving the SDP without block diagonalisation takes about 4.5 minutes, *the block diagonalised SDP takes about 10 seconds, a 25x speedup!* \n", "\n", "While this is a significant speedup, it should be noted that this is notable only for problems with a large number of symmetries, such as the previous example which was an inflation with 3 copies of each source. For problems with moment matrices of similar dimension but with less symmetries, the speedup is present but not as significant. For example, the same problem but with (2,2,2) copies per source and NPA level 2 union local level 1 generates a moment matrix of 873 columns. Solving the original SDP on one computer takes around 2.3 minutes, while solving the block diagonalised SDP takes around 1.7 minutes, giving a speedup of around 1.3x. While still relevant, it is not as significant as the previous example." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.16" }, "vscode": { "interpreter": { "hash": "5185cb8988fc84c35117c94793cda6c5f0bb6718bc4f8ace0826abbce28c3e20" } } }, "nbformat": 4, "nbformat_minor": 4 }