From a3744756c34043736ea7d263d2d4f566366c326f Mon Sep 17 00:00:00 2001 From: Josh Combes Date: Sat, 15 Jun 2019 07:09:58 +1000 Subject: [PATCH] Operator tools (#140) This is a major module re-org plus adding new module to check plain old operators are unitary etc. * move the random_operators file to the module operator_tools * move the superoperator tools * fix the legend not showing and add padding to the title * superoperator_transformations * notebook update * type checking and docstring updates * add superoperator_transformations * tests * superoperator_validation * apply superop * proj superops * unused imports * add and remove files * validate operators * improvements * test superop transforms * test apply superop * test channel approx * first pass at operator_tools refactor * Spruce up permute_tesnor_factor and expand functionality. * Move project_state_matrix_to_physical out of tomography. * Fix import * Update names and imports for project_state_matrix * Remove nonsense code. * compose kraus + tests * improve * big update * format improvments * update * to fix the failed test add # NBVAL_IGNORE_OUTPUT * it should have been # NBVAL_RAISES_EXCEPTION * Typo fixing. * matplotlib >= 3.1.0 --- examples/random_operators.ipynb | 242 +++++-- examples/superoperator_tools.ipynb | 660 +++++++++++++++--- examples/tomography_state.ipynb | 6 +- .../benchmarking/operator_tools/__init__.py | 8 + .../operator_tools/apply_superoperator.py | 82 +++ .../operator_tools/channel_approximation.py | 52 ++ .../operator_tools/compose_superoperators.py | 42 ++ .../operator_tools/project_state_matrix.py | 52 ++ .../operator_tools/project_superoperators.py | 145 ++++ .../operator_tools/random_operators.py | 209 ++++++ .../superoperator_transformations.py | 376 ++++++++++ .../operator_tools/validate_operator.py | 171 +++++ .../operator_tools/validate_superoperator.py | 155 ++++ forest/benchmarking/plotting/state_process.py | 3 +- .../tests/test_apply_superoperator.py | 29 + .../tests/test_channel_approximation.py | 34 + .../tests/test_compose_superoperators.py | 32 + .../tests/test_project_state_matrix.py | 14 + .../tests/test_project_superoperators.py | 92 +++ .../tests/test_random_operators.py | 26 +- .../tests/test_state_tomography.py | 15 +- ... => test_superoperator_transformations.py} | 198 ------ .../tests/test_validate_operators.py | 77 ++ .../tests/test_validate_superoperator.py | 60 ++ forest/benchmarking/tomography.py | 57 +- requirements.txt | 2 +- 26 files changed, 2433 insertions(+), 406 deletions(-) create mode 100644 forest/benchmarking/operator_tools/__init__.py create mode 100644 forest/benchmarking/operator_tools/apply_superoperator.py create mode 100644 forest/benchmarking/operator_tools/channel_approximation.py create mode 100644 forest/benchmarking/operator_tools/compose_superoperators.py create mode 100644 forest/benchmarking/operator_tools/project_state_matrix.py create mode 100644 forest/benchmarking/operator_tools/project_superoperators.py create mode 100644 forest/benchmarking/operator_tools/random_operators.py create mode 100644 forest/benchmarking/operator_tools/superoperator_transformations.py create mode 100644 forest/benchmarking/operator_tools/validate_operator.py create mode 100644 forest/benchmarking/operator_tools/validate_superoperator.py create mode 100644 forest/benchmarking/tests/test_apply_superoperator.py create mode 100644 forest/benchmarking/tests/test_channel_approximation.py create mode 100644 forest/benchmarking/tests/test_compose_superoperators.py create mode 100644 forest/benchmarking/tests/test_project_state_matrix.py create mode 100644 forest/benchmarking/tests/test_project_superoperators.py rename forest/benchmarking/tests/{test_superoperator_tools.py => test_superoperator_transformations.py} (57%) create mode 100644 forest/benchmarking/tests/test_validate_operators.py create mode 100644 forest/benchmarking/tests/test_validate_superoperator.py diff --git a/examples/random_operators.ipynb b/examples/random_operators.ipynb index ce16698e..f2df678c 100644 --- a/examples/random_operators.ipynb +++ b/examples/random_operators.ipynb @@ -4,7 +4,16 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Random Operators: examples of random states and channels" + "# Random Operators: random quantum states and channels" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this notebook we explore a submodule of `operator_tools` called `random_operators`.\n", + "\n", + "In the context of forest benchmarking the primary use of random operators is to test the estimation routines. For example you might modify an existing state or process tomography routine (or develop a new method) and want to test that your modification works. One way to do that would be to test it on a bunch of random quantum states or channels. " ] }, { @@ -14,16 +23,27 @@ "outputs": [], "source": [ "import numpy as np\n", - "from sympy.combinatorics import Permutation\n", - "import forest.benchmarking.random_operators as rand_ops\n", - "import numpy.linalg" + "import forest.benchmarking.operator_tools.random_operators as rand_ops" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Complex Ginibre ensemble" + "## Random Operators" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Complex Ginibre ensemble\n", + "\n", + "This is a subroutine for other methods in the module. The complex Ginibre ensemble is a random matrix where the real and imaginary parts of each entry, $G(n,m)$, are drawn in an IID fashion from $\\mathcal N(0,1)$ e.g.\n", + "\n", + "$$G(n,m) = X(n,m) + i Y(n,m)$$\n", + "\n", + "where $X(n,m), Y(n,m)\\sim \\mathcal N(0,1)$. For our purpose we allow for non square matricies." ] }, { @@ -32,14 +52,20 @@ "metadata": {}, "outputs": [], "source": [ - "rand_ops.ginibre_matrix_complex(2,2)" + "gini_2by2 = rand_ops.ginibre_matrix_complex(2,2)\n", + "print(np.round(gini_2by2,3))\n", + "\n", + "print('Notice that the above matrix is not Hermitian.')\n", + "print('We can explicitly test if it is Hermitian: ', np.all(gini_2by2.T.conj()==gini_2by2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Haar random unitary" + "### Haar random unitary\n", + "\n", + "Here you simply specify the dimesion of the Hilbert space. " ] }, { @@ -49,16 +75,47 @@ "outputs": [], "source": [ "U = rand_ops.haar_rand_unitary(2)\n", + "print(U)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can test to see how close to unitary it is:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ "print(np.around(U.dot(np.transpose(np.conj(U))),decimals=15))\n", - "print(np.around(U.dot(np.transpose(np.conj(U))),decimals=16))\n", - "# only good to 16 decimal places..." + "print(np.around(U.dot(np.transpose(np.conj(U))),decimals=16))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Apparently it is only good to 16 decimal places." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Haar random pure state " + "## Random States" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Haar random pure state\n", + "\n", + "The simpliest random state is a state drawn from the Haar measure. It is a pure state, i.e. the purity is $P(\\rho)={\\rm Tr}[\\rho^2]=1$. These states are generated by applying a Haar random unitary to a fixed fuduial state, usually $|0\\ldots0\\rangle$." ] }, { @@ -67,15 +124,21 @@ "metadata": {}, "outputs": [], "source": [ - "psi = rand_ops.haar_rand_state(2)\n", - "print(psi)" + "psi2 = rand_ops.haar_rand_state(2)\n", + "print('The state vector is \\n', np.round(psi2,3))\n", + "print('It has shape', psi2.shape,'and purity P = ', np.real(np.round(np.trace(psi2@psi2.T.conj()@psi2@psi2.T.conj()),2)))\n", + "print('\\n')\n", + "print('Now lets look at a random pure state on two qubits.')\n", + "psi4 = rand_ops.haar_rand_state(4)\n", + "print('The state vector is \\n', np.round(psi4,3), )\n", + "print('It has shape', psi4.shape,'.')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Ginibre State (mixed state with rank K)" + "For fun lets plot the state in the Pauli representation." ] }, { @@ -84,11 +147,37 @@ "metadata": {}, "outputs": [], "source": [ - "# mixed single qubit state\n", - "print(np.around(rand_ops.ginibre_state_matrix(2,2),4))\n", - "print(\"----------------------\")\n", - "# mixed two qubit state\n", - "print(np.around(rand_ops.ginibre_state_matrix(4,2),4))" + "from forest.benchmarking.plotting.state_process import plot_pauli_rep_of_state\n", + "from forest.benchmarking.operator_tools.superoperator_transformations import computational2pauli_basis_matrix, vec \n", + "from forest.benchmarking.utils import n_qubit_pauli_basis\n", + "import matplotlib.pyplot as plt\n", + "\n", + "# change of basis\n", + "n_qubits = 2\n", + "pl_basis_twoq = n_qubit_pauli_basis(n_qubits)\n", + "c2p_twoq = computational2pauli_basis_matrix(2*n_qubits)\n", + "\n", + "# turn a state vector into a state matrix\n", + "rho = psi4@psi4.T.conj()\n", + "# convert the state to the Pauli representation which should be real\n", + "state = np.real(c2p_twoq@vec(rho))\n", + "\n", + "fig, ax = plt.subplots()\n", + "plot_pauli_rep_of_state(state.transpose(), ax, pl_basis_twoq.labels, 'Random Two qubit state in Pauli representation')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Ginibre State (mixed state with rank K)\n", + "\n", + "This function lets us generate mixed states with a specific rank.\n", + "\n", + "Specifically, given a Hilbert space dimension $D$ and a desired rank $K$, this function \n", + "gives a D by D positive semidefinite matrix of rank $K$ drawn from the Ginibre ensemble. \n", + " \n", + "For $D = K$ these are states drawn from the **Hilbert-Schmidt measure**." ] }, { @@ -97,15 +186,21 @@ "metadata": {}, "outputs": [], "source": [ - "# you cant have Rank > Hilbert space Dim\n", - "rand_ops.ginibre_state_matrix(2,3)" + "print('This is a mixed single qubit state drawn from Hilbert-Schmidt measure (as D=K)')\n", + "print(np.around(rand_ops.ginibre_state_matrix(2,2),3))\n", + "print(\"\\n\")\n", + "print('This is a mixed two qubit state with rank 2:')\n", + "print(np.around(rand_ops.ginibre_state_matrix(4,2),3))\n", + "evals, evec = np.linalg.eig(rand_ops.ginibre_state_matrix(4,2))\n", + "print('\\n')\n", + "print('Here are the eigenvalues:', np.round(evals,3),'. You can see only two are non zero.')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## State from Bures measure" + "### State from Bures measure" ] }, { @@ -114,14 +209,23 @@ "metadata": {}, "outputs": [], "source": [ - "rand_ops.bures_measure_state_matrix(2)" + "np.round(rand_ops.bures_measure_state_matrix(2),4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Uniform ensemble of CPTP maps" + "## Random quantum Channels" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Uniform ensemble of CPTP maps (BCSZ distribution)\n", + "\n", + "Given a Hilbert space dimension $D$ and a Kraus rank $K$, this function returns a $D^2 × D^2$ Choi matrix $J(Λ)$ of a channel drawn from the BCSZ distribution with Kraus rank $K$." ] }, { @@ -130,9 +234,9 @@ "metadata": {}, "outputs": [], "source": [ - "# random quantum channel on one qubit in Choi form\n", - "choi = rand_ops.rand_map_with_BCSZ_dist(2,2)\n", - "choi" + "rand_choi = rand_ops.rand_map_with_BCSZ_dist(2,2)\n", + "print('Here is a random quantum channel on one qubit in Choi form:')\n", + "np.round(rand_choi,3)" ] }, { @@ -141,7 +245,45 @@ "metadata": {}, "outputs": [], "source": [ - "choi.shape" + "rand_choi.shape" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To convert to different superoperator representations we import the module `operator_tools.superoperator_transformations`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import forest.benchmarking.operator_tools.superoperator_transformations as sot\n", + "print('We can convert this channel to Kraus form and enumerate the Kraus operators. We expect there to be two Kraus ops, consistent with the rank we specified.')\n", + "\n", + "for idx, kraus_op in enumerate(sot.choi2kraus(rand_choi)):\n", + " print('Kraus OP #'+str(1+idx)+' is: \\n', np.round(kraus_op,3))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rand_pl = sot.choi2pauli_liouville(rand_choi)\n", + "\n", + "from forest.benchmarking.plotting.state_process import plot_pauli_transfer_matrix\n", + "n_qubits = 1\n", + "pl_basis_oneq = n_qubit_pauli_basis(n_qubits)\n", + "c2p_oneq = computational2pauli_basis_matrix(2*n_qubits)\n", + "\n", + "f, (ax1) = plt.subplots(1, 1, figsize=(5.5, 4.2))\n", + "plot_pauli_transfer_matrix(np.real(rand_pl), ax1, pl_basis_oneq.labels, 'The Pauli transfer matrix of a random CPTP channel')\n", + "plt.show()" ] }, { @@ -158,26 +300,37 @@ "outputs": [], "source": [ "# pick a hilbert space dimension\n", - "D = 2\n", - "# pick a way you want to permute the operators\n", - "perm =[1,2,0]\n", - "# Note: the number of elements in the permutation determines\n", - "# the number of Hilbert spaces you are considering." + "D = 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "^^ here the Identity permutation is P = [0,1,2] which maps (a,b,c) to (a,b,c).\n", - "The permutation P = [1,2,0] maps (a,b,c) to (b,c,a)." + "Lets consider a tensor product of three Hilbert spaces: $$\\mathcal H_a \\otimes \\mathcal H_b\\otimes \\mathcal H_c.$$ \n", + "\n", + "Next we need to pick a way you want to permute the operators; specified by a permutation in [cycle notation](https://en.wikipedia.org/wiki/Permutation#Cycle_notation). \n", + "\n", + "For example the Identity permutation is $P = [0,1,2]$ which maps $(a,b,c)$ to $(a,b,c)$.\n", + "The permutation $P = [1,2,0]$ maps $(a,b,c)$ to $(b,c,a)$, so lets try that." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "perm =[1,2,0]\n", + "# Note: the number of elements in the permutation determines \n", + "# the number of Hilbert spaces you are considering." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Create the basis states in the Hilbert space" + "Create the basis states in the Hilbert space" ] }, { @@ -198,7 +351,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Create inital state and answer after applying the permutation [1,2,0]" + "Create inital state and answer after applying the permutation [1,2,0]" ] }, { @@ -207,15 +360,18 @@ "metadata": {}, "outputs": [], "source": [ - "inital_vector = np.kron(np.kron(states[0],states[0]),states[1]) # before permuting anything\n", - "perm_vector = np.kron(np.kron(states[0],states[1]),states[0]) # apply the permutation by hand" + "# before permuting anything\n", + "inital_vector = np.kron(np.kron(states[0],states[0]),states[1]) \n", + "\n", + "# apply the permutation by hand\n", + "perm_vector = np.kron(np.kron(states[0],states[1]),states[0]) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### create permutation operator" + "**create permutation operator**" ] }, { @@ -231,7 +387,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### check the permutation operator applied to the intial vector gives the correct answer" + "**check the permutation operator applied to the intial vector gives the correct answer**" ] }, { @@ -249,7 +405,7 @@ "metadata": {}, "outputs": [], "source": [ - "np.matmul(perm_vector.T,answer)" + "print('The inner product between the calculated and true answer is one', np.matmul(perm_vector.T,answer))" ] }, { @@ -261,10 +417,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "", - "name": "" - }, "language_info": { "name": "python", "pygments_lexer": "ipython3" diff --git a/examples/superoperator_tools.ipynb b/examples/superoperator_tools.ipynb index 13cff7af..9d965d75 100644 --- a/examples/superoperator_tools.ipynb +++ b/examples/superoperator_tools.ipynb @@ -4,7 +4,85 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Superoperator conversion: example conversion between representations" + "# Superoperator tools\n", + "\n", + "In this notebook we explore the submodules of `operator_tools` that enable easy manipulation of the various quantum channel representations.\n", + "\n", + "To summarize the functionality:\n", + "- vectorization and conversions between different repesentations of quantum channels\n", + "- apply quantum operations\n", + "- compose quantum operations\n", + "- validate that quantum channels are physical\n", + "- project unphysical channels to physical channels" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Breif motivation and introduction" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Perfect gates in **reversible classical computation** are described by permutation matricies, e.g. the [Toffoli gate](https://en.wikipedia.org/wiki/Toffoli_gate), while the input states are vectors. A noisy classical gate could be modeled as a perfect gate followed by a noise channel, e.g. [binary symmetric channel](https://en.wikipedia.org/wiki/Binary_symmetric_channel), on all the bits in the state vector.\n", + "\n", + "Perfect gates in **quantum computation** are described by unitary matricies and states are described by complex vectors, e.g.\n", + "$$\n", + "|\\psi\\rangle = U |\\psi_0\\rangle.\n", + "$$\n", + "\n", + "Modeling **noisy quantum computation** often makes use of [mixed states](https://en.wikipedia.org/wiki/Density_matrix) and quantum operations or quantum noise channels.\n", + "\n", + "Interestingly there are a number of ways to represent quantum noise channels, and depending on your task some can be more convenient than others. The simplest case to illustrate this point is to consider a mixed initial state $\\rho$ undergoing unitary evolution\n", + "\n", + "$$ \\rho' = U \\rho U^\\dagger.$$\n", + "\n", + "The fact that the unitary has to act on both sides of the inital state means it is a [*superoperator*](https://en.wikipedia.org/wiki/Superoperator), that is an object that can act on operators like the state matrix. \n", + "\n", + "\n", + "\n", + "It turns out using a special matrix multiiplication identity we can write this as\n", + "$$ |\\rho'\\rangle \\rangle = \\mathcal U |\\rho\\rangle\\rangle,\n", + "$$\n", + "where $\\mathcal U = U^*\\otimes U$ and $|\\rho\\rangle\\rangle = {\\rm vec}(\\rho)$. The nice thing about this is it looks like the pure state case. This is because the operator (the state) has become a vector and the superoperator (the left right action of $U$) has become an operator. \n", + "\n", + "\n", + "\n", + "**More information** \n", + "Below we will assume that you are already an expert in these topics. If you are unfamilar with these topics we recomend the following references\n", + "- chapter 8 of [Mike_N_Ike] which is on *Quantum noise and quantum operations*. \n", + "- chapter 3 of John Preskill's lecture notes [Physics 219/Computer Science 219](http://www.theory.caltech.edu/people/preskill/ph219/chap3_15.pdf)\n", + "- the file in `/docs/Superoperator representations.md` \n", + "- for an intuitive but advanced treatment see [GRAPTN]\n", + "\n", + "\n", + "\n", + "[Mike_N_Ike] Quantum Computation and Quantum Information \n", + " Michael A. Nielsen & Isaac L. Chuang \n", + " Cambridge: Cambridge University Press (2000) \n", + "\n", + "\n", + "[GRAPTN] Tensor networks and graphical calculus for open quantum systems \n", + " Christopher Wood et al. \n", + " Quant. Inf. Comp. 15, 0579-0811 (2015) \n", + " (no DOI) \n", + " https://arxiv.org/abs/1111.6950 " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Conversion between different descriptions of quantum channels\n", + "\n", + "We intentionally chose not to make quantum channels python objects with methods that would automatically transform between represtations. \n", + "\n", + "The functions to convert between different representations are called things like `kraus2chi`, `kraus2choi`, `pauli_liouville2choi` etc.\n", + "\n", + "This assumes the user does not do silly things like input a choi matrix to a function `chi2choi`." ] }, { @@ -14,7 +92,38 @@ "outputs": [], "source": [ "import numpy as np\n", - "from forest.benchmarking.superoperator_tools import vec, unvec" + "from pyquil.gate_matrices import I, X, Y, Z, H, CNOT" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define some channels" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def amplitude_damping_kraus(p):\n", + " Ad0 = np.asarray([[1, 0], [0, np.sqrt(1 - p)]])\n", + " Ad1 = np.asarray([[0, np.sqrt(p)], [0, 0]])\n", + " return [Ad0, Ad1]\n", + "\n", + "def bit_flip_kraus(p):\n", + " M0 = np.sqrt(1 - p) * I\n", + " M1 = np.sqrt(p) * X\n", + " return [M0, M1]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define some states" ] }, { @@ -23,21 +132,38 @@ "metadata": {}, "outputs": [], "source": [ - "I = np.asarray([[1, 0], [0, 1]])\n", - "X = np.asarray([[0, 1], [1, 0]])\n", - "Y = np.asarray([[0, -1.0j], [1.0j, 0]])\n", - "Z = np.asarray([[1, 0], [0, -1]])\n", - "H = np.asarray([[1, 1], [1, -1]]) / np.sqrt(2)\n", - "CNOT = np.array([[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]])\n", "one_state = np.asarray([[0,0],[0,1]])\n", - "zero_state = np.asarray([[1,0],[0,0]])" + "zero_state = np.asarray([[1,0],[0,0]])\n", + "rho_mixed = np.asarray([[0.9,0],[0,0.1]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# vec and unvec " + "### vec and unvec \n", + "\n", + "We can vectorize i.e. `vec` and unvec matricies.\n", + "\n", + "We chose a column stacking convention so that the matrix\n", + "$$\n", + "A = \\begin{pmatrix} 1 & 2\\\\ 3 & 4\\end{pmatrix}\n", + "$$\n", + "becomes\n", + "$$\n", + "|A\\rangle\\rangle = {\\rm vec}(A) = \\begin{pmatrix} 1\\\\ 3\\\\ 2\\\\ 4\\end{pmatrix}.\n", + "$$\n", + "\n", + "Let's check that" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from forest.benchmarking.operator_tools import vec, unvec" ] }, { @@ -59,14 +185,14 @@ "print(\" \")\n", "print(vec(A))\n", "print(\" \")\n", - "print(unvec(vec(A))==A)" + "print('Does the story check out? ', np.all(unvec(vec(A))==A))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# Apply Kraus operators to the excited state" + "### Kraus to $\\chi$ matrix (aka chi or process matrix)" ] }, { @@ -75,25 +201,121 @@ "metadata": {}, "outputs": [], "source": [ - "one_state = np.asarray([[0,0],[0,1]])\n", - "one = vec(one_state)\n", - "Ad0 = np.asarray([[1, 0], [0, np.sqrt(1-0.1)]])\n", - "Ad1 = np.asarray([[0, np.sqrt(0.1)], [0, 0]])\n", - "ko = [Ad0, Ad1]\n", - "print(ko)\n", - "# check kraus ops sum to ID\n", - "np.matmul(Ad0.transpose(), Ad0) + np.matmul(Ad1.transpose(), Ad1)" + "from forest.benchmarking.operator_tools import kraus2chi" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Lets do a unitary gate first, say the Hadamard" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print('The Kraus operator is:\\n', np.round(H,3))\n", + "print('\\n')\n", + "print('The Chi matrix is:\\n', kraus2chi(H))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now consider the Amplitude damping channel" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "AD_kraus = amplitude_damping_kraus(0.1)\n", + "\n", + "print('The Kraus operators are:\\n', np.round(AD_kraus,3))\n", + "print('\\n')\n", + "print('The Chi matrix is:\\n', np.round(kraus2chi(AD_kraus),3))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Kraus to Pauli Liouville aka the \"Pauli Transfer Matrix\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from forest.benchmarking.operator_tools import kraus2pauli_liouville" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Hpaulirep = kraus2pauli_liouville(H)\n", + "Hpaulirep" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Use Kraus operators to find out put of channel i.e.\n", + "We can visualize this using the tools from the plotting module." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from forest.benchmarking.plotting.state_process import plot_pauli_transfer_matrix\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "f, (ax1) = plt.subplots(1, 1, figsize=(5, 4.2))\n", "\n", - "$\\rho_{out} = A_0 \\rho A_0^\\dagger + A_1 \\rho A_1^\\dagger$\n", + "plot_pauli_transfer_matrix(Hpaulirep,ax=ax1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The above figure is a graphical representaiton of:\n", + " \n", + "(out operator) = H (in operator) H\n", + "\n", + "Z = H X H \n", + "-Y = H Y H \n", + "X = H Z H " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Evolving states using quantum channels\n", "\n", - "Except we only transpose as these matrices have real entries" + "In many superoperator representations evolution coresponds to mutiplying the vec'ed state by the superoperator. E.g." ] }, { @@ -102,23 +324,55 @@ "metadata": {}, "outputs": [], "source": [ - "rho_out = np.matmul(np.matmul(Ad0, one_state), Ad0.transpose()) \\\n", - " + np.matmul(np.matmul(Ad1, one_state), Ad1.transpose())\n", - "print(rho_out)" + "from forest.benchmarking.superoperator_tools import kraus2superop\n", + "\n", + "zero_state_vec = vec(zero_state)\n", + "\n", + "answer_vec = np.matmul(kraus2superop([H]), zero_state_vec)\n", + "\n", + "print('The vec\\'ed answer is', answer_vec)\n", + "print('\\n')\n", + "print('The unvec\\'ed answer is\\n', np.real(unvec(answer_vec)))\n", + "print('\\n')\n", + "print('Let\\'s compare it to the normal calculation\\n', H @ zero_state @ H)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Kraus to superoperator" + "For representations with this simple application there are no inbuilt functions in forest benchmarking. \n", + "\n", + "However applying a channel is more painful in the Choi and Kraus representation.\n", + "\n", + "Consider the amplitude damping channel where we need to perform the following calcualtion to find out put of channel \n", + "$\\rho_{out} = A_0 \\rho A_0^\\dagger + A_1 \\rho A_1^\\dagger.$\n", + "We provide helper functions to do these calculations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from forest.benchmarking.operator_tools import apply_kraus_ops_2_state, apply_choi_matrix_2_state, kraus2choi" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "apply_kraus_ops_2_state(AD_kraus, one_state)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Do the case of one Kraus operator i.e. a unitary." + "In the Choi representation we get the same answer:" ] }, { @@ -127,7 +381,18 @@ "metadata": {}, "outputs": [], "source": [ - "from forest.benchmarking.superoperator_tools import kraus2superop" + "AD_choi = kraus2choi(AD_kraus)\n", + "\n", + "apply_choi_matrix_2_state(AD_choi, one_state)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Compose quantum channels\n", + "\n", + "Composing channels is useful when describing larger circuits. In some representations e.g. in the superoperator or Liouville representation it is just matrix multiplication e.g." ] }, { @@ -136,10 +401,7 @@ "metadata": {}, "outputs": [], "source": [ - "zero_state_vec = vec(zero_state)\n", - "print(zero_state)\n", - "print(\"\")\n", - "print(zero_state_vec)" + "from forest.benchmarking.operator_tools import superop2kraus, kraus2superop" ] }, { @@ -148,15 +410,25 @@ "metadata": {}, "outputs": [], "source": [ - "# Hadamard * |0>\n", - "unvec(np.matmul(kraus2superop([H]), zero_state_vec))" + "H_super = kraus2superop(H)\n", + "\n", + "H_squared_super = H_super @ H_super\n", + "\n", + "print('Hadamard squared as a superoperator:\\n', np.round(H_squared_super,2))\n", + "\n", + "print('\\n As a Kraus operator:\\n', np.round(superop2kraus(H_squared_super),2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Do the case of amplitude damping i.e. two Kraus operators." + "Composing channels in the Kraus represntaion is more difficult. Consider composing two channels $\\mathcal A$ (with Kraus operators $[A_0, A_1]$) and $\\mathcal B$ (with Kraus operators $[B_0, B_1]$). The composition is \n", + "\n", + "$$\\begin{align}\n", + "\\mathcal B(\\mathcal A(\\rho)) & = \\sum_i \\sum_j B_j A_i \\rho A_i^\\dagger B_j^\\dagger \n", + "\\end{align}\n", + "$$" ] }, { @@ -165,11 +437,7 @@ "metadata": {}, "outputs": [], "source": [ - "# super op makes sense\n", - "one_state_vec = vec(one_state) \n", - "print(\"vec'ed version \\n\", np.matmul(kraus2superop(ko), one_state_vec))\n", - "print(\" \")\n", - "print(\"unvec'ed i.e. rho_out \\n\", unvec(np.matmul(kraus2superop(ko), one_state_vec)))" + "from forest.benchmarking.operator_tools import compose_channel_kraus, superop2kraus" ] }, { @@ -178,14 +446,16 @@ "metadata": {}, "outputs": [], "source": [ - "print(kraus2superop([X]))" + "BitFlip_kraus = bit_flip_kraus(0.2)\n", + "\n", + "kraus2superop(compose_channel_kraus(AD_kraus, BitFlip_kraus))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# Kraus to Choi" + "This is the same as if we do" ] }, { @@ -194,9 +464,36 @@ "metadata": {}, "outputs": [], "source": [ - "from scipy.linalg import expm\n", - "from forest.benchmarking.superoperator_tools import kraus2choi\n", - "from forest.benchmarking.utils import partial_trace" + "BitFlip_super = kraus2superop(BitFlip_kraus)\n", + "AD_super = kraus2superop(AD_kraus)\n", + "\n", + "AD_super @ BitFlip_super" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can also easily compose channels acting on independent spaces.\n", + "\n", + "Consider composing the same two channels as above, $\\mathcal A$ and $\\mathcal B$. However this time they act on different Hilbert spaces. With respect to the tensor product structure $H_2 \\otimes H_1$ the Kraus operators are $[A_0\\otimes I, A_1\\otimes I]$ and $[I \\otimes B_0, I \\otimes B_1]$.\n", + "\n", + "In this case the order of the operations commutes \n", + "$$\\begin{align}\n", + "\\mathcal A(\\mathcal B(\\rho))= \\mathcal B(\\mathcal A(\\rho)) & = \\sum_i \\sum_j A_i\\otimes B_j \\rho A_i^\\dagger\\otimes B_j^\\dagger \n", + "\\end{align}\n", + "$$\n", + "\n", + "In forest benchmarking you can specify the two channels without the Identity tensored on and it will take care of it for you:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from forest.benchmarking.operator_tools import tensor_channel_kraus" ] }, { @@ -205,8 +502,27 @@ "metadata": {}, "outputs": [], "source": [ - "choi = kraus2choi(ko)\n", - "print(choi)" + "np.round(tensor_channel_kraus(AD_kraus,BitFlip_kraus),3)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Validate quantum channels are physical\n", + "\n", + "When doing process tomography sometimes the estimates returned by various estimation methods can result in unphysical processes.\n", + "\n", + "The functions below can be used to check if the estimates are physical.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As a starting point, we might want to check if a process specified by Kraus operators is valid. \n", + "\n", + "Unless a process is unitary you need more than one Kraus operator to be a valid quantum operation." ] }, { @@ -215,9 +531,16 @@ "metadata": {}, "outputs": [], "source": [ - "prody = np.matmul(choi, np.kron(one_state.transpose(),I) )\n", - "print(partial_trace(prody, [1], [2,2]))\n", - "print(\"\\n which is the same as rho_out\")\n" + "from forest.benchmarking.operator_tools import kraus_operators_are_valid\n", + "\n", + "kraus_operators_are_valid(AD_kraus[0])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "However a full set is valid:" ] }, { @@ -226,8 +549,14 @@ "metadata": {}, "outputs": [], "source": [ - "choiCNOT = kraus2choi([CNOT])\n", - "print(choiCNOT)" + "kraus_operators_are_valid(AD_kraus)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can also validate other properties of quantum channels such as completely positivity and trace preservation. This is done on the **Choi** representation, so you many need to convert your quantum operation to the Choi representaiton first.\n" ] }, { @@ -236,52 +565,209 @@ "metadata": {}, "outputs": [], "source": [ - "two_q_id = np.kron(I,I)\n", + "from forest.benchmarking.operator_tools import (choi_is_unitary, \n", + " choi_is_unital, \n", + " choi_is_trace_preserving, \n", + " choi_is_completely_positive, \n", + " choi_is_cptp)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# amplitude damping is not unitary\n", + "print(choi_is_unitary(AD_choi),'\\n')\n", "\n", - "print(\"Lets see if the CNOT action is correct\")\n", - "print(\" \")\n", + "# amplitude damping is not unital\n", + "print(choi_is_unital(AD_choi))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# amplitude damping is trace preserving (TP)\n", + "print(choi_is_trace_preserving(AD_choi),'\\n')\n", "\n", - "print(\"one one state\")\n", - "input_state = np.kron(one_state, one_state)\n", - "prody = np.matmul(choiCNOT, np.kron(two_q_id, input_state.transpose()))\n", - "print(partial_trace(prody, [0,1], [2,2,2,2]))\n", + "# amplitude damping is completely positive (CP)\n", + "print(choi_is_completely_positive(AD_choi), '\\n')\n", "\n", - "print(\"one zero state\")\n", - "input_stated = np.kron(one_state, zero_state)\n", - "prody = np.matmul(choiCNOT, np.kron(two_q_id, input_stated.transpose()))\n", - "print(partial_trace(prody, [0,1], [2,2,2,2]))\n", + "# amplitude damping is CPTP\n", + "print(choi_is_cptp(AD_choi))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Project unphysical channels to physical channels\n", "\n", - "print(\"zero one state\")\n", - "input_stated = np.kron(zero_state,one_state)\n", - "prody = np.matmul(choiCNOT, np.kron(two_q_id, input_stated.transpose()))\n", - "print(partial_trace(prody, [0,1], [2,2,2,2]))\n", + "When doing process tomography often the estimates returned by maximum likelihood estimation or linear inversion methods can result in unphysical processes.\n", "\n", - "print(\"zero zero state\")\n", - "input_stated = np.kron(zero_state, zero_state)\n", - "prody = np.matmul(choiCNOT, np.kron(two_q_id, input_stated.transpose()))\n", - "print(partial_trace(prody, [0,1], [2,2,2,2]))\n", + "The functions below can be used to project the unphysical estimates back to physical estimates." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from forest.benchmarking.operator_tools.project_superoperators import (proj_choi_to_completely_positive,\n", + " proj_choi_to_trace_non_increasing,\n", + " proj_choi_to_trace_preserving,\n", + " proj_choi_to_physical)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "neg_Id_choi = -kraus2choi(I)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "proj_choi_to_completely_positive(neg_Id_choi)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "proj_choi_to_trace_non_increasing(neg_Id_choi)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "proj_choi_to_trace_preserving(neg_Id_choi)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "proj_choi_to_physical(neg_Id_choi)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Validate operators\n", + "\n", + "A lot of the work in validating the physicality of quantum channels comes down to validating properties of matricies:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from forest.benchmarking.operator_tools.validate_operator import (is_square_matrix, \n", + " is_identity_matrix, \n", + " is_idempotent_matrix, \n", + " is_unitary_matrix, \n", + " is_positive_semidefinite_matrix)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# a vector is not square\n", + "is_square_matrix(np.array([[1], [0]]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# NBVAL_RAISES_EXCEPTION\n", + "# the line above is for testing purposes, do not remove.\n", + "\n", + "# a tensor is not a matrix\n", + "\n", + "tensor = np.ones(8).reshape(2,2,2)\n", + "print(tensor)\n", + "\n", + "is_square_matrix(tensor)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "is_identity_matrix(X)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "projector_zero = np.array([[1, 0], [0, 0]])\n", "\n", - "print(\"\\n Story checks out.\")" + "is_idempotent_matrix(projector_zero)" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "is_unitary_matrix(AD_kraus[0])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "is_positive_semidefinite_matrix(I)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { - "kernelspec": { - "display_name": "Python 3", - "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.6.3" + "pygments_lexer": "ipython3" } }, "nbformat": 4, diff --git a/examples/tomography_state.ipynb b/examples/tomography_state.ipynb index b09f503d..a1f2df6c 100644 --- a/examples/tomography_state.ipynb +++ b/examples/tomography_state.ipynb @@ -265,10 +265,10 @@ "metadata": {}, "outputs": [], "source": [ - "from forest.benchmarking.tomography import project_density_matrix\n", + "from forest.benchmarking.operator_tools.project_state_matrix import project_state_matrix_to_physical\n", "rho_unphys = np.random.uniform(-1, 1, (2, 2)) \\\n", " * np.exp(1.j * np.random.uniform(-np.pi, np.pi, (2, 2)))\n", - "rho_phys = project_density_matrix(rho_unphys)\n", + "rho_phys = project_state_matrix_to_physical(rho_unphys)\n", "\n", "fig, (ax1, ax2) = plt.subplots(1, 2)\n", "hinton(rho_unphys, ax=ax1)\n", @@ -288,7 +288,7 @@ "# https://doi.org/10.1103/PhysRevLett.108.070502\n", "\n", "eigs = np.diag(np.array(list(reversed([3.0/5, 1.0/2, 7.0/20, 1.0/10, -11.0/20]))))\n", - "phys = project_density_matrix(eigs)\n", + "phys = project_state_matrix_to_physical(eigs)\n", "np.allclose(phys, np.diag([0, 0, 1.0/5, 7.0/20, 9.0/20]))" ] }, diff --git a/forest/benchmarking/operator_tools/__init__.py b/forest/benchmarking/operator_tools/__init__.py new file mode 100644 index 00000000..f576bcb4 --- /dev/null +++ b/forest/benchmarking/operator_tools/__init__.py @@ -0,0 +1,8 @@ +from .apply_superoperator import * +from .channel_approximation import * +from .compose_superoperators import * +from .project_superoperators import * +from .superoperator_transformations import * +from .random_operators import * +from .validate_operator import * +from .validate_superoperator import * diff --git a/forest/benchmarking/operator_tools/apply_superoperator.py b/forest/benchmarking/operator_tools/apply_superoperator.py new file mode 100644 index 00000000..b687b238 --- /dev/null +++ b/forest/benchmarking/operator_tools/apply_superoperator.py @@ -0,0 +1,82 @@ +"""A module containing tools for applying superoperators to states. + +We have arbitrarily decided to use a column stacking convention. + +For more information about the conventions used, look at the file in +/docs/Superoperator representations.md + +Further references include: + +[GRAPTN] Tensor networks and graphical calculus for open quantum systems + Wood et al. + Quant. Inf. Comp. 15, 0579-0811 (2015) + (no DOI) + https://arxiv.org/abs/1111.6950 + +[MATQO] On the Matrix Representation of Quantum Operations + Nambu et al., + arXiv: 0504091 (2005) + https://arxiv.org/abs/quant-ph/0504091 + +[DUAL] On duality between quantum maps and quantum states + Zyczkowski et al., + Open Syst. Inf. Dyn. 11, 3 (2004) + https://dx.doi.org/10.1023/B:OPSY.0000024753.05661.c2 + https://arxiv.org/abs/quant-ph/0401119 + +""" +from typing import Sequence +import numpy as np +from forest.benchmarking.utils import partial_trace + + +def apply_kraus_ops_2_state(kraus_ops: Sequence[np.ndarray], state: np.ndarray) -> np.ndarray: + r""" + Apply a quantum channel, specified by Kraus operators, to state. + + The Kraus operators need not be square. + + :param kraus_ops: A list or tuple of N Kraus operators, each operator is M by dim ndarray + :param state: A dim by dim ndarray which is the density matrix for the state + :return: M by M ndarray which is the density matrix for the state after the action of kraus_ops + """ + if isinstance(kraus_ops, np.ndarray): # handle input of single kraus op + if len(kraus_ops[0].shape) < 2: # first elem is not a matrix + kraus_ops = [kraus_ops] + + dim, _ = state.shape + rows, cols = kraus_ops[0].shape + + if dim != cols: + raise ValueError("Dimensions of state and Kraus operator are incompatible") + + new_state = np.zeros((rows, rows)) + for M in kraus_ops: + new_state += M @ state @ np.transpose(M.conj()) + + return new_state + + +def apply_choi_matrix_2_state(choi: np.ndarray, state: np.ndarray) -> np.ndarray: + r""" + Apply a quantum channel, specified by a Choi matrix (using the column stacking convention), + to a state. + + The Choi matrix is a dim**2 by dim**2 matrix and the state rho is a dim by dim matrix. The + output state is + + rho_{out} = Tr_{A_{in}}[(rho^T \otimes Id) Choi_matrix ], + + where T denotes transposition and Tr_{A_{in}} is the partial trace over input Hilbert space H_{ + A_{in}}; the Choi matrix representing a process mapping rho in H_{A_{in}} to rho_{out} + in H_{B_{out}} is regarded as an operator on the space H_{A_{in}} \otimes H_{B_{out}}. + + + :param choi: a dim**2 by dim**2 matrix + :param state: A dim by dim ndarray which is the density matrix for the state + :return: a dim by dim matrix. + """ + dim = int(np.sqrt(np.asarray(choi).shape[0])) + dims = [dim, dim] + tot_matrix = np.kron(state.transpose(), np.identity(dim)) @ choi + return partial_trace(tot_matrix, [1], dims) diff --git a/forest/benchmarking/operator_tools/channel_approximation.py b/forest/benchmarking/operator_tools/channel_approximation.py new file mode 100644 index 00000000..f3aa9cf8 --- /dev/null +++ b/forest/benchmarking/operator_tools/channel_approximation.py @@ -0,0 +1,52 @@ +"""A module containing tools for approximating channels. + +We have arbitrarily decided to use a column stacking convention. + +For more information about the conventions used, look at the file in +/docs/Superoperator representations.md + +Further references include: + +[GRAPTN] Tensor networks and graphical calculus for open quantum systems + Wood et al. + Quant. Inf. Comp. 15, 0579-0811 (2015) + (no DOI) + https://arxiv.org/abs/1111.6950 + +[MATQO] On the Matrix Representation of Quantum Operations + Nambu et al., + arXiv: 0504091 (2005) + https://arxiv.org/abs/quant-ph/0504091 + +[DUAL] On duality between quantum maps and quantum states + Zyczkowski et al., + Open Syst. Inf. Dyn. 11, 3 (2004) + https://dx.doi.org/10.1023/B:OPSY.0000024753.05661.c2 + https://arxiv.org/abs/quant-ph/0401119 + +""" +import numpy as np + + +def pauli_twirl_chi_matrix(chi_matrix: np.ndarray) -> np.ndarray: + r""" + Implements a Pauli twirl of a chi matrix (aka a process matrix). + + See the folloiwng reference for more details + + [SPICC] Scalable protocol for identification of correctable codes + Silva et al., + PRA 78, 012347 2008 + http://dx.doi.org/10.1103/PhysRevA.78.012347 + https://arxiv.org/abs/0710.1900 + + Note: Pauli twirling a quantum channel can give rise to a channel that is less noisy; use with + care. + + :param chi_matrix: a dim**2 by dim**2 chi or process matrix + :return: dim**2 by dim**2 chi or process matrix + """ + return np.diag(chi_matrix.diagonal()) + + +# TODO: Honest approximations for Channels that act on one or MORE qubits. diff --git a/forest/benchmarking/operator_tools/compose_superoperators.py b/forest/benchmarking/operator_tools/compose_superoperators.py new file mode 100644 index 00000000..d4ba2ecc --- /dev/null +++ b/forest/benchmarking/operator_tools/compose_superoperators.py @@ -0,0 +1,42 @@ +"""A module containing tools for composing superoperators. +""" +from typing import Sequence +import numpy as np + + +def tensor_channel_kraus(k2: Sequence[np.ndarray], + k1: Sequence[np.ndarray]) -> Sequence[np.ndarray]: + r""" + Given the Kraus representaion for two channels, $\mathcal E$ and $\mathcal F$, acting on + different systems this function returns the Kraus operators representing the composition of + these independent channels. + + Suppose $\mathcal E$ and $\mathcal F$ both have one Kraus operator K_1 = X and K_2 = H, + that is they are unitary. Then, with respect to the tensor product structure + + $H_2 \otimes H_1$ + + of the individual systems this function returns $K_{\rm tot} = H \otimes X$. + + :param k1: The list of Kraus operators on the first system. + :param k2: The list of Kraus operators on the second system. + :return: A list of tensored Kraus operators. + """ + # TODO: make this function work for an arbitrary number of Kraus operators + return [np.kron(k2l, k1j) for k1j in k1 for k2l in k2] + + +def compose_channel_kraus(k2: Sequence[np.ndarray], + k1: Sequence[np.ndarray]) -> Sequence[np.ndarray]: + r""" + Given two channels, K_1 and K_2, acting on the same system in the Kraus representation this + function return the Kraus operators representing the composition of the channels. + + It is assumed that K_1 is applied first then K_2 is applied. + + :param k2: The list of Kraus operators that are applied second. + :param k1: The list of Kraus operators that are applied first. + :return: A combinatorially generated list of composed Kraus operators. + """ + # TODO: make this function work for an arbitrary number of Kraus operators + return [np.dot(k2l, k1j) for k1j in k1 for k2l in k2] diff --git a/forest/benchmarking/operator_tools/project_state_matrix.py b/forest/benchmarking/operator_tools/project_state_matrix.py new file mode 100644 index 00000000..132ae76a --- /dev/null +++ b/forest/benchmarking/operator_tools/project_state_matrix.py @@ -0,0 +1,52 @@ +import numpy as np +import functools +from scipy.linalg import eigh + + +def project_state_matrix_to_physical(rho: np.ndarray) -> np.ndarray: + """ + Project a possibly unphysical estimated density matrix to the closest (with respect to the + 2-norm) positive semi-definite matrix with trace 1, that is a valid quantum state. + + This is the so called "wizard" method. It is described in the following reference: + + [MLEWIZ] Efficient Method for Computing the Maximum-Likelihood Quantum State from + Measurements with Additive Gaussian Noise + Smolin et al., + Phys. Rev. Lett. 108, 070502 (2012) + https://doi.org/10.1103/PhysRevLett.108.070502 + https://arxiv.org/abs/1106.5458 + + :param rho: the density (state) matrix with shape (N, N) + :return rho_projected: The closest positive semi-definite trace 1 matrix to rho. + """ + # Rescale to trace 1 if the matrix is not already + rho_impure = rho / np.trace(rho) + + dimension = rho_impure.shape[0] # the dimension of the Hilbert space + [eigvals, eigvecs] = eigh(rho_impure) + + # If matrix is already trace one PSD, we are done + if np.min(eigvals) >= 0: + return rho_impure + + # Otherwise, continue finding closest trace one, PSD matrix + eigvals = list(eigvals) + eigvals.reverse() + eigvals_new = [0.0] * len(eigvals) + + i = dimension + accumulator = 0.0 # Accumulator + while eigvals[i - 1] + accumulator / float(i) < 0: + accumulator += eigvals[i - 1] + i -= 1 + for j in range(i): + eigvals_new[j] = eigvals[j] + accumulator / float(i) + eigvals_new.reverse() + + # Reconstruct the matrix + rho_projected = functools.reduce(np.dot, (eigvecs, + np.diag(eigvals_new), + np.conj(eigvecs.T))) + + return rho_projected diff --git a/forest/benchmarking/operator_tools/project_superoperators.py b/forest/benchmarking/operator_tools/project_superoperators.py new file mode 100644 index 00000000..71fc6922 --- /dev/null +++ b/forest/benchmarking/operator_tools/project_superoperators.py @@ -0,0 +1,145 @@ +"""A module containing tools for projecting superoperators to CP, TNI, TP, and physical. + +We have arbitrarily decided to use a column stacking convention. + +A good reference for these methods is: + +[PGD] Maximum-likelihood quantum process tomography via projected gradient descent + Knee et al., + Phys. Rev. A 98, 062336 (2018) + https://dx.doi.org/10.1103/PhysRevA.98.062336 + https://arxiv.org/abs/1803.10062 +""" +import numpy as np +from forest.benchmarking.utils import partial_trace +from forest.benchmarking.operator_tools.superoperator_transformations import vec + + +def proj_choi_to_completely_positive(choi: np.ndarray) -> np.ndarray: + """ + Projects the Choi representation of a process into the nearest Choi matrix in the space of + completely positive maps. + + Equation 8 of [PGD] + + [PGD] Maximum-likelihood quantum process tomography via projected gradient descent + Knee et al., + Phys. Rev. A 98, 062336 (2018) + https://dx.doi.org/10.1103/PhysRevA.98.062336 + https://arxiv.org/abs/1803.10062 + + :param choi: Choi representation of a process + :return: closest Choi matrix in the space of completely positive maps + """ + hermitian = (choi + choi.conj().T) / 2 # enforce Hermiticity + evals, v = np.linalg.eigh(hermitian) + evals[evals < 0] = 0 # enforce completely positive by removing negative eigenvalues + diag = np.diag(evals) + return v @ diag @ v.conj().T + + +def proj_choi_to_trace_non_increasing(choi: np.ndarray) -> np.ndarray: + """ + Projects the Choi matrix of a process into the space of trace non-increasing maps. + + Equation 33 of [PGD] + + :param choi: Choi representation of a process + :return: Choi representation of the projected trace non-increasing process + """ + dim = int(np.sqrt(choi.shape[0])) + + # trace out the output Hilbert space + pt = partial_trace(choi, dims=[dim, dim], keep=[0]) + + hermitian = (pt + pt.conj().T) / 2 # enforce Hermiticity + d, v = np.linalg.eigh(hermitian) + d[d > 1] = 1 # enforce trace preserving + D = np.diag(d) + projection = v @ D @ v.conj().T + + trace_increasing_part = np.kron((pt - projection) / dim, np.eye(dim)) + + return choi - trace_increasing_part + + +def proj_choi_to_trace_preserving(choi: np.ndarray) -> np.ndarray: + """ + Projects the Choi representation of a process to the closest processes in the space of trace + preserving maps. + + Equation 12 of [PGD], but without vecing the Choi matrix. See choi_is_trace_preserving for + comparison. + + :param choi: Choi representation of a process + :return: Choi representation of the projected trace preserving process + """ + dim = int(np.sqrt(choi.shape[0])) + + # trace out the output Hilbert space, keep the input space at index 0 + pt = partial_trace(choi, dims=[dim, dim], keep=[0]) + # isolate the part the violates the condition we want, namely pt = Id + diff = pt - np.eye(dim) + # we want to subtract off the violation from the larger operator, so 'invert' the partial_trace + subtract = np.kron(diff/dim, np.eye(dim)) + return choi - subtract + + +def proj_choi_to_physical(choi: np.ndarray, make_trace_preserving: bool = True) -> np.ndarray: + """ + Projects the given Choi matrix into the subspace of Completetly Positive and either + Trace Perserving (TP) or Trace-Non-Increasing maps. + + Uses Dykstra's algorithm with the stopping criterion presented in: + + [DYKALG] Dykstra’s algorithm and robust stopping criteria + Birgin et al., + (Springer US, Boston, MA, 2009), pp. 828–833, ISBN 978-0-387-74759-0. + https://doi.org/10.1007/978-0-387-74759-0_143 + + This method is suggested in [PGD] + + :param choi: the Choi representation estimate of a quantum process. + :param make_trace_preserving: default true, projects the estimate to a trace-preserving + process. If false the output process may only be trace non-increasing + :return: The Choi representation of the Completely Positive, Trace Preserving (CPTP) or Trace + Non-Increasing map that is closest to the given state. + """ + old_CP_change = np.zeros_like(choi) + old_TP_change = np.zeros_like(choi) + last_CP_projection = np.zeros_like(choi) + last_state = choi + + while True: + # Dykstra's algorithm + pre_CP = last_state - old_CP_change + CP_projection = proj_choi_to_completely_positive(pre_CP) + new_CP_change = CP_projection - pre_CP + + pre_TP = CP_projection - old_TP_change + if make_trace_preserving: + new_state = proj_choi_to_trace_preserving(pre_TP) + else: + new_state = proj_choi_to_trace_non_increasing(pre_TP) + new_TP_change = new_state - pre_TP + + CP_change_change = new_CP_change - old_CP_change + TP_change_change = new_TP_change - old_TP_change + state_change = new_state - last_state + + # stopping criterion + # norm(mat) is the frobenius norm + # norm(mat)**2 is thus equivalent to the dot product vec(mat) dot vec(mat) + if np.linalg.norm(CP_change_change) ** 2 + np.linalg.norm(TP_change_change) ** 2 \ + + 2 * abs(np.dot(vec(old_TP_change).conj().T, vec(state_change))) \ + + 2 * abs(np.dot(vec(old_CP_change).conj().T, + vec(CP_projection - last_CP_projection))) < 1e-4: + break + + # store results from this iteration + old_CP_change = new_CP_change + old_TP_change = new_TP_change + last_CP_projection = CP_projection + last_state = new_state + + return new_state diff --git a/forest/benchmarking/operator_tools/random_operators.py b/forest/benchmarking/operator_tools/random_operators.py new file mode 100644 index 00000000..5f64082a --- /dev/null +++ b/forest/benchmarking/operator_tools/random_operators.py @@ -0,0 +1,209 @@ +"""A module for generating random quantum states and processes. + +Pseudocode for many of these routines can be found in the appendix of the paper: + +[BAYES] Practical Bayesian Tomography + Granade et al., + New Journal of Physics 18, 033024 (2016) + https://dx.doi.org/10.1088/1367-2630/18/3/033024 + https://arxiv.org/abs/1509.03770 +""" +from typing import Optional, List, Union + +import numpy as np +from numpy import linalg as la +from scipy.linalg import sqrtm +from sympy.combinatorics import Permutation +from numpy.random import RandomState +from forest.benchmarking.utils import partial_trace + + +def ginibre_matrix_complex(dim: int, k: int, rs: Optional[RandomState] = None) -> np.ndarray: + r""" + Given a scalars dim and k, returns a dim by k matrix, drawn from the complex Ginibre + ensemble, i.e. each element is distributed ~ [N(0, 1) + i · N(0, 1)]. Here X ~ N(0,1) + denotes a normally distributed random variable. + + [IM] Induced measures in the space of mixed quantum states + Zyczkowski et al., + J. Phys A: Math. and Gen. 34, 7111 (2001) + https://doi.org/10.1088/0305-4470/34/35/335 + https://arxiv.org/abs/quant-ph/0012101 + + :param dim: Hilbert space dimension. + :param k: Ultimately becomes the rank of a state. + :param rs: Optional random state. + :return: Returns a dim by k matrix, drawn from the Ginibre ensemble. + """ + if rs is None: + rs = np.random + return rs.randn(dim, k) + 1j * rs.randn(dim, k) + + +def haar_rand_unitary(dim: int, rs=None) -> np.ndarray: + """ + Given a Hilbert space dimension dim this function + returns a unitary operator U ∈ C^(dim by dim) drawn from the Haar measure. + + The error is of order 10^-16. + + [MEZ] How to generate random matrices from the classical compact groups + Mezzadri + Notices of the American Mathematical Society 54, 592 (2007). + http://www.ams.org/notices/200705/fea-mezzadri-web.pdf + https://arxiv.org/abs/math-ph/0609050 + + :param dim: Hilbert space dimension (scalar). + :param rs: Optional random state + :return: Returns a dim by dim unitary operator U drawn from the Haar measure. + """ + if rs is None: + rs = np.random + Z = ginibre_matrix_complex(dim=dim, k=dim, rs=rs) # /np.sqrt(2) + Q, R = np.linalg.qr(Z) + diag = np.diagonal(R) + lamb = np.diag(diag) / np.absolute(diag) + return np.matmul(Q, lamb) + + +def haar_rand_state(dim: int) -> np.ndarray: + """ + Given a Hilbert space dimension dim this function returns a vector + representing a random pure state operator drawn from the Haar measure. + + :param dim: Hilbert space dimension. + :return: Returns a dim by 1 vector drawn from the Haar measure. + + """ + unitary = haar_rand_unitary(dim) + fiducial_vec = np.zeros((dim, 1)) + fiducial_vec[0] = 1 + return np.matmul(unitary, fiducial_vec) + + +def ginibre_state_matrix(dim: int, rank: int) -> np.ndarray: + """ + Given a Hilbert space dimension dim and a desired rank K, returns a dim by dim positive + semidefinite matrix of rank K drawn from the Ginibre ensemble. For dim = K these are states + drawn from the Hilbert-Schmidt measure. + + See reference [IM] for more details. + + :param dim: Hilbert space dimension. + :param rank: The rank of a state. + :return: Returns a dim by rank matrix, drawn from the Ginibre ensemble. + """ + if rank > dim: + raise ValueError("The rank of the state matrix cannot exceed the dimension.") + A = ginibre_matrix_complex(dim, rank) + M = A.dot(np.transpose(np.conjugate(A))) + return M / np.trace(M) + + +def bures_measure_state_matrix(dim: int) -> np.ndarray: + """ + Given a Hilbert space dimension dim, returns a dim by dim positive semidefinite matrix drawn + from the Bures measure. + + [OSZ] Random Bures mixed states and the distribution of their purity + Osipov et al., + J. Phys. A: Math. Theor. 43, 055302 (2010). + https://doi.org/10.1088/1751-8113/43/5/055302 + https://arxiv.org/abs/0909.5094 + + :param dim: Hilbert space dimension. + :return: Returns a dim by dim matrix, drawn from the Bures measure. + """ + A = ginibre_matrix_complex(dim, dim) + U = haar_rand_unitary(dim) + Udag = np.transpose(np.conjugate(U)) + Id = np.eye(dim) + M = A.dot(np.transpose(np.conjugate(A))) + P = (Id + U).dot(M).dot(Id + Udag) + return P / np.trace(P) + + +def rand_map_with_BCSZ_dist(dim: int, kraus_rank: int) -> np.ndarray: + """ + Given a Hilbert space dimension dim and a Kraus rank K, returns a $dim^2 by dim^2$ Choi + matrix $J(Λ)$ of a channel drawn from the BCSZ distribution with Kraus rank $K$. + + [RQO] Random quantum operations, + Bruzda et al., + Physics Letters A 373, 320 (2009). + https://doi.org/10.1016/j.physleta.2008.11.043 + https://arxiv.org/abs/0804.2361 + + :param dim: Hilbert space dimension. + :param kraus_rank: The number of Kraus operators in the operator sum description of the channel. + :return: dim^2 by dim^2 Choi matrix, drawn from the BCSZ distribution with Kraus rank K. + """ + # TODO: this ^^ is CPTP, might want a flag that allows for just CP quantum operations. + X = ginibre_matrix_complex(dim ** 2, kraus_rank) + rho = X @ X.conj().T + rho_red = partial_trace(rho, [0], [dim, dim]) + # Note that Eqn. 8 of [RQO] uses a *row* stacking convention so in that case we would write + # Q = np.kron(np.eye(D), sqrtm(la.inv(rho_red))) + # But as we use column stacking we need: + Q = np.kron(sqrtm(la.inv(rho_red)), np.eye(dim)) + Z = Q @ rho @ Q + return Z + + +def permute_tensor_factors(dims: Union[int, List[int]], perm: List[int]) -> np.ndarray: + r""" + Return a permutation matrix that appropriately swaps spaces of the given dimension(s). + + Given a Hilbert space dimension dim and an list representing the permutation perm of the + tensor product Hilbert spaces, returns a $dim^len(perm)$ by $dim^len(perm)$ permutation matrix. + + E.g. 1) Suppose dims=2 and perm=[0,1] + Returns the identity operator on two qubits + + 2) Suppose dims=2 and perm=[1,0] + Returns the SWAP operator on two qubits which + maps A_1 \otimes A_2 --> A_2 \otimes A_1. + + 3) Suppose dims=[2, 4] and perm=[1,0] + Returns the SWAP operator on three qubits which + maps A_1 \otimes (A_2 \otimes A_3) -> (A_2 \otimes A_3) + + See: Equations 5.11, 5.12, and 5.13 in + + [SCOTT] Optimizing quantum process tomography with unitary 2-designs + A. J. Scott, + J. Phys. A 41, 055308 (2008) + https://dx.doi.org/10.1088/1751-8113/41/5/055308 + https://arxiv.org/abs/0711.1017 + + This function is used in tests for other functions. However, it can also be useful when + thinking about higher moment (N>2) integrals over the Haar measure. + + :param dims: The dimension of each Hilbert space factor given in order of the pre-permuted + factorization of the total space. If an int is specified then each space in the + permutation is assumed to have this dimension. + :param perm: A list representing the permutation of the tensor factors. + :return: a permutation matrix that permutes the factors of the given dimension. + """ + if isinstance(dims, int): + dim_list = [dims for _ in range(2 * len(perm))] + total_dim = dims ** len(perm) + else: + assert len(dims) == len(perm), "Please specify the dimension of each factor to be permuted." + dim_list = [dim for _ in range(2) for dim in dims] + total_dim = np.prod(dims) + + transpositions = Permutation(perm).transpositions() + + # start with identity + perm_matrix = np.eye(total_dim, total_dim) + + # reshape for easy implementation of transpositions + perm_matrix = np.reshape(perm_matrix, dim_list) + + # build up the permutation one transposition at a time + for swap in transpositions: + perm_matrix = np.swapaxes(perm_matrix, swap[0], swap[1]) + + # reshape to act on total_dim space + return np.reshape(perm_matrix, [total_dim, total_dim]) diff --git a/forest/benchmarking/operator_tools/superoperator_transformations.py b/forest/benchmarking/operator_tools/superoperator_transformations.py new file mode 100644 index 00000000..512a69cc --- /dev/null +++ b/forest/benchmarking/operator_tools/superoperator_transformations.py @@ -0,0 +1,376 @@ +"""A module containing tools for converting between different representations of superoperators. + +We have arbitrarily decided to use a column stacking convention. + +For more information about the conventions used, look at the file in +/docs/Superoperator representations.md + +Further references include: + +[GRAPTN] Tensor networks and graphical calculus for open quantum systems + Wood et al. + Quant. Inf. Comp. 15, 0579-0811 (2015) + (no DOI) + https://arxiv.org/abs/1111.6950 + +[MATQO] On the Matrix Representation of Quantum Operations + Nambu et al., + arXiv: 0504091 (2005) + https://arxiv.org/abs/quant-ph/0504091 + +[DUAL] On duality between quantum maps and quantum states + Zyczkowski et al., + Open Syst. Inf. Dyn. 11, 3 (2004) + https://dx.doi.org/10.1023/B:OPSY.0000024753.05661.c2 + https://arxiv.org/abs/quant-ph/0401119 + +""" +from typing import Sequence, Tuple, List +import numpy as np +from forest.benchmarking.utils import n_qubit_pauli_basis + + +def vec(matrix: np.ndarray) -> np.ndarray: + """ + Vectorize, i.e. "vec", a matrix by column stacking. + + For example the 2 by 2 matrix A + + A = [[a, b] + [c, d]] becomes |A>> := vec(A) = (a, c, b, d)^T , + + where |A>> denotes the vec'ed version of A and T denotes transpose. + + :param matrix: A N (rows) by M (columns) numpy array. + :return: Returns a column vector with N by M rows. + """ + return np.asarray(matrix).T.reshape((-1, 1)) + + +def unvec(vector: np.ndarray, shape: Tuple[int, int] = None) -> np.ndarray: + """ + Take a column vector and turn it into a matrix. + + By default, the unvec'ed matrix is assumed to be square. Specifying shape = [N, M] will + produce a N by M matrix where N is the number of rows and M is the number of columns. + + Consider |A>> := vec(A) = (a, c, b, d)^T. `unvec(|A>>)` should return + + A = [[a, b] + [c, d]]. + + :param vector: A (N*M) by 1 numpy array. + :param shape: The shape of the output matrix; by default, the matrix is assumed to be square. + :return: Returns a N by M matrix. + """ + vector = np.asarray(vector) + if shape is None: + dim = int(np.sqrt(vector.size)) + shape = dim, dim + matrix = vector.reshape(*shape).T + return matrix + + +def kraus2chi(kraus_ops: Sequence[np.ndarray]) -> np.ndarray: + r""" + Convert a set of Kraus operators (representing a channel) to + a chi matrix which is also known as a process matrix. + + :param kraus_ops: A list or tuple of N Kraus operators + :return: Returns a dim**2 by dim**2 matrix. + """ + if isinstance(kraus_ops, np.ndarray): # handle input of single kraus op + if len(kraus_ops[0].shape) < 2: # first elem is not a matrix + kraus_ops = [kraus_ops] + + dim = np.asarray(kraus_ops[0]).shape[0] # kraus op is dim by dim matrix + c_vecs = [computational2pauli_basis_matrix(dim) @ vec(kraus) for kraus in kraus_ops] + chi_mat = sum([c_vec @ c_vec.conj().T for c_vec in c_vecs]) + return chi_mat + + +def kraus2superop(kraus_ops: Sequence[np.ndarray]) -> np.ndarray: + r""" + Convert a set of Kraus operators (representing a channel) to + a superoperator using the column stacking convention. + + Suppose the N Kraus operators M_i are dim by dim matrices. Then the + the superoperator is a (dim**2) \otimes (dim**2) matrix. Using the relation + column stacking relation, + vec(ABC) = (C^T\otimes A) vec(B), we can show + + super_operator = \sum_i^N ( M_i^\dagger )^T \otimes M_i + = \sum_i^N M_i^* \otimes M_i + + where A^* is the complex conjugate of a matrix A, A^T is the transpose, + and A^\dagger is the complex conjugate and transpose. + + Note: This function can also convert non-square Kraus operators to a superoperator, + these frequently arise in quantum measurement theory and quantum error correction. In that + situation consider a single Kraus operator that is M by N then the superoperator will be a + M**2 by N**2 matrix. + + :param kraus_ops: A tuple of N Kraus operators + :return: Returns a dim**2 by dim**2 matrix. + """ + if isinstance(kraus_ops, np.ndarray): # handle input of single kraus op + if len(kraus_ops[0].shape) < 2: # first elem is not a matrix + kraus_ops = [kraus_ops] + + rows, cols = np.asarray(kraus_ops[0]).shape + + # Standard case of square Kraus operators is if rows==cols. + # When representing a partial projection, e.g. a single measurement operator + # M_i = Id \otimes np.ndarray: + """ + Convert a set of Kraus operators (representing a channel) to + a Pauli-Liouville matrix (aka Pauli Transfer matrix). + + :param kraus_ops: A list of Kraus operators + :return: Returns dim**2 by dim**2 Pauli-Liouville matrix + """ + return superop2pauli_liouville(kraus2superop(kraus_ops)) + + +def kraus2choi(kraus_ops: Sequence[np.ndarray]) -> np.ndarray: + r""" + Convert a set of Kraus operators (representing a channel) to + a Choi matrix using the column stacking convention. + + Suppose the N Kraus operators M_i are dim by dim matrices. Then the + the Choi matrix is a dim**2 by dim**2 matrix + + choi_matrix = \sum_i^N |M_i>> (|M_i>>)^\dagger + = \sum_i^N |M_i>> <> = vec(M_i) + + :param kraus_ops: A list of N Kraus operators + :return: Returns a dim**2 by dim**2 matrix. + """ + if isinstance(kraus_ops, np.ndarray): # handle input of single kraus op + if len(kraus_ops[0].shape) < 2: # first elem is not a matrix + kraus_ops = [kraus_ops] + + return sum([vec(op) @ vec(op).conj().T for op in kraus_ops]) + + +def chi2pauli_liouville(chi_matrix: np.ndarray) -> np.ndarray: + r""" + Converts a chi matrix (aka a process matrix) to the Pauli Liouville representation. + + :param chi_matrix: a dim**2 by dim**2 process matrix + :return: dim**2 by dim**2 Pauli-Liouville matrix + """ + return choi2pauli_liouville(chi2choi(chi_matrix)) + + +def chi2kraus(chi_matrix: np.ndarray) -> List[np.ndarray]: + """ + Converts a chi matrix into a list of Kraus operators. (operators with small norm may be + excluded) + + :param chi_matrix: a dim**2 by dim**2 process matrix + :return: list of Kraus operators + """ + return pauli_liouville2kraus(chi2pauli_liouville(chi_matrix)) + + +def chi2superop(chi_matrix: np.ndarray) -> np.ndarray: + """ + Converts a chi matrix into a superoperator. + + :param chi_matrix: a dim**2 by dim**2 process matrix + :return: a dim**2 by dim**2 superoperator matrix + """ + return pauli_liouville2superop(chi2pauli_liouville(chi_matrix)) + + +def chi2choi(chi_matrix: np.ndarray) -> np.ndarray: + """ + Converts a chi matrix into a Choi matrix. + + :param chi_matrix: a dim**2 by dim**2 process matrix + :return: a dim**2 by dim**2 Choi matrix + """ + dim = int(np.sqrt(np.asarray(chi_matrix).shape[0])) + p2c = pauli2computational_basis_matrix(dim) + return p2c @ chi_matrix @ p2c.conj().T + + +def superop2kraus(superop: np.ndarray) -> List[np.ndarray]: + """ + Converts a superoperator into a list of Kraus operators. (operators with small norm may be excluded) + + :param superop: a dim**2 by dim**2 superoperator + :return: list of Kraus operators + """ + return choi2kraus(superop2choi(superop)) + + +def superop2chi(superop: np.ndarray) -> np.ndarray: + """ + Converts a superoperator into a list of Kraus operators. (operators with small norm may be excluded) + + :param superop: a dim**2 by dim**2 superoperator + :return: a dim**2 by dim**2 process matrix + """ + return kraus2chi(superop2kraus(superop)) + + +def superop2pauli_liouville(superop: np.ndarray) -> np.ndarray: + """ + Converts a superoperator into a pauli_liouville matrix. This is achieved by a linear change of basis. + + :param superop: a dim**2 by dim**2 superoperator + :return: dim**2 by dim**2 Pauli-Liouville matrix + """ + dim = int(np.sqrt(np.asarray(superop).shape[0])) + c2p_basis_transform = computational2pauli_basis_matrix(dim) + return c2p_basis_transform @ superop @ c2p_basis_transform.conj().T * dim + + +def superop2choi(superop: np.ndarray) -> np.ndarray: + """ + Convert a superoperator into a choi matrix. The operation acts equivalently to choi2superop, as it is a bijection. + + :param superop: a dim**2 by dim**2 superoperator + :return: dim**2 by dim**2 choi matrix + """ + dim = int(np.sqrt(np.asarray(superop).shape[0])) + return np.reshape(superop, [dim] * 4).swapaxes(0, 3).reshape([dim ** 2, dim ** 2]) + + +def pauli_liouville2kraus(pl_matrix: np.ndarray) -> List[np.ndarray]: + """ + Converts a pauli_liouville matrix into a list of Kraus operators. (operators with small norm may be excluded) + + :param pl_matrix: a dim**2 by dim**2 pauli_liouville matrix + :return: list of Kraus operators + """ + return choi2kraus(pauli_liouville2choi(pl_matrix)) + + +def pauli_liouville2chi(pl_matrix: np.ndarray) -> np.ndarray: + """ + Converts a pauli_liouville matrix into a chi matrix. (operators with small norm may be excluded) + + :param pl_matrix: a dim**2 by dim**2 pauli_liouville matrix + :return: a dim**2 by dim**2 process matrix + """ + return kraus2chi(pauli_liouville2kraus(pl_matrix)) + + +def pauli_liouville2superop(pl_matrix: np.ndarray) -> np.ndarray: + """ + Converts a pauli_liouville matrix into a superoperator. This is achieved by a linear change of basis. + + :param pl_matrix: a dim**2 by dim**2 Pauli-Liouville matrix + :return: dim**2 by dim**2 superoperator + """ + dim = int(np.sqrt(np.asarray(pl_matrix).shape[0])) + p2c_basis_transform = pauli2computational_basis_matrix(dim) + return p2c_basis_transform @ pl_matrix @ p2c_basis_transform.conj().T / dim + + +def pauli_liouville2choi(pl_matrix: np.ndarray) -> np.ndarray: + """ + Convert a Pauli-Liouville matrix into a choi matrix. + + :param pl_matrix: a dim**2 by dim**2 Pauli-Liouville matrix + :return: dim**2 by dim**2 choi matrix + """ + return superop2choi(pauli_liouville2superop(pl_matrix)) + + +def choi2kraus(choi: np.ndarray, tol: float = 1e-9) -> List[np.ndarray]: + """ + Converts a Choi matrix into a list of Kraus operators. (operators with small norm may be + excluded) + + :param choi: a dim**2 by dim**2 choi matrix + :param tol: optional threshold parameter for eigenvalues/kraus ops to be discarded + :return: list of Kraus operators + """ + eigvals, v = np.linalg.eigh(choi) + return [np.lib.scimath.sqrt(eigval) * unvec(np.array([evec]).T) for eigval, evec in + zip(eigvals, v.T) if abs(eigval) > tol] + + +def choi2chi(choi: np.ndarray) -> np.ndarray: + """ + Converts a Choi matrix into a chi matrix. (operators with small norm may be excluded) + :param choi: a dim**2 by dim**2 choi matrix + :return: a dim**2 by dim**2 process matrix + """ + return kraus2chi(choi2kraus(choi)) + + +def choi2superop(choi: np.ndarray) -> np.ndarray: + """ + Convert a choi matrix into a superoperator. The operation acts equivalently to superop2choi, as it is a bijection. + + :param choi: a dim**2 by dim**2 choi matrix + :return: dim**2 by dim**2 superoperator + """ + dim = int(np.sqrt(np.asarray(choi).shape[0])) + return np.reshape(choi, [dim] * 4).swapaxes(0, 3).reshape([dim ** 2, dim ** 2]) + + +def choi2pauli_liouville(choi: np.ndarray) -> np.ndarray: + """ + Convert a choi matrix into a Pauli-Liouville matrix. + + :param choi: a dim**2 by dim**2 choi matrix + :return: dim**2 by dim**2 Pauli-Liouville matrix + """ + return superop2pauli_liouville(choi2superop(choi)) + + +def pauli2computational_basis_matrix(dim) -> np.ndarray: + """ + Produces a basis transform matrix that converts from a pauli basis to the computational basis. + p2c_transform = sum_{k=1}^{dim**2} | sigma_k >> > + + :param dim: dimension of the hilbert space on which the operators act. + :return: A dim**2 by dim**2 basis transform matrix + """ + n_qubits = int(np.log2(dim)) + + conversion_mat = np.zeros((dim ** 2, dim ** 2), dtype=complex) + + for i, pauli in enumerate(n_qubit_pauli_basis(n_qubits)): + pauli_label = np.zeros((dim ** 2, 1)) + pauli_label[i] = 1. + pauli_mat = pauli[1] + conversion_mat += np.kron(vec(pauli_mat), pauli_label.T) + + return conversion_mat + + +def computational2pauli_basis_matrix(dim) -> np.ndarray: + """ + Produces a basis transform matrix that converts from a computational basis to a pauli basis. + Conjugate transpose of pauli2computational_basis_matrix with an extra dimensional factor. + c2p_transform = sum_{k=1}^{dim**2} | k > << sigma_k | + For example, + vec(sigma_z) = | sigma_z >> = [1, 0, 0, -1].T in the computational basis + c2p * | sigma_z >> = [0, 0, 0, 1].T + + :param dim: dimension of the hilbert space on which the operators act. + :return: A dim**2 by dim**2 basis transform matrix + """ + return pauli2computational_basis_matrix(dim).conj().T / dim diff --git a/forest/benchmarking/operator_tools/validate_operator.py b/forest/benchmarking/operator_tools/validate_operator.py new file mode 100644 index 00000000..ec2bd740 --- /dev/null +++ b/forest/benchmarking/operator_tools/validate_operator.py @@ -0,0 +1,171 @@ +"""A module allowing one to check properties of operators or matrices. +""" +import numpy as np + + +def is_square_matrix(matrix: np.ndarray) -> bool: + """ + Checks if a matrix is square. + + :param matrix: a M by N matrix. + :param rtol: The relative tolerance parameter in np.allclose + :param atol: The absolute tolerance parameter in np.allclose + :return: True if the matrix is square; False otherwise. + """ + if len(matrix.shape) != 2: + raise ValueError("The object is not a matrix.") + rows, cols = matrix.shape + return rows == cols + + +def is_symmetric_matrix(matrix: np.ndarray, rtol: float = 1e-05, atol: float = 1e-08) -> bool: + """ + Checks if a square matrix is symmetric. That is + + A is symmetric iff $A = A ^T$, + + where $T$ denotes transpose. + + :param matrix: a M by M matrix. + :param rtol: The relative tolerance parameter in np.allclose + :param atol: The absolute tolerance parameter in np.allclose + :return: True if the matrix is symmetric; False otherwise. + """ + if not is_square_matrix(matrix): + raise ValueError("The matrix is not square.") + return np.allclose(matrix, matrix.T, rtol=rtol, atol=atol) + + +def is_identity_matrix(matrix: np.ndarray, rtol: float = 1e-05, atol: float = 1e-08) -> bool: + """ + Checks if a square matrix is the identity matrix. + + :param matrix: a M by M matrix. + :param rtol: The relative tolerance parameter in np.allclose + :param atol: The absolute tolerance parameter in np.allclose + :return: True if the matrix is the identity matrix; False otherwise. + """ + if not is_square_matrix(matrix): + raise ValueError("The matrix is not square.") + Id = np.eye(len(matrix)) + return np.allclose(matrix, Id, rtol=rtol, atol=atol) + + +def is_idempotent_matrix(matrix: np.ndarray, rtol: float = 1e-05, atol: float = 1e-08) -> bool: + """ + Checks if a square matrix A is idempotent. That is + + A is idempotent iff $A^2 = A.$ + + :param matrix: a M by M matrix. + :param rtol: The relative tolerance parameter in np.allclose + :param atol: The absolute tolerance parameter in np.allclose + :return: True if the matrix is the idempotent; False otherwise. + """ + if not is_square_matrix(matrix): + raise ValueError("The matrix is not square.") + return np.allclose(matrix, matrix @ matrix, rtol=rtol, atol=atol) + + +def is_normal_matrix(matrix: np.ndarray, rtol: float = 1e-05, atol: float = 1e-08) -> bool: + r""" + Checks if a square matrix A is normal. That is + + A is normal iff $A^{\dagger} A = A A^{\dagger}$, + + where $\dagger$ denotes conjugate transpose. + + :param matrix: a M by M matrix. + :param rtol: The relative tolerance parameter in np.allclose + :param atol: The absolute tolerance parameter in np.allclose + :return: True if the matrix is normal; False otherwise. + """ + if not is_square_matrix(matrix): + raise ValueError("The matrix is not square.") + AB = matrix.T.conj() @ matrix + BA = matrix @ matrix.T.conj() + return np.allclose(AB, BA, rtol=rtol, atol=atol) + + +def is_hermitian_matrix(matrix: np.ndarray, rtol: float = 1e-05, atol: float = 1e-08) -> bool: + r""" + Checks if a square matrix A is Hermitian. That is + + A is Hermitian iff $A = A^{\dagger}$, + + where $\dagger$ denotes conjugate transpose. + + :param matrix: a M by M matrix. + :param rtol: The relative tolerance parameter in np.allclose + :param atol: The absolute tolerance parameter in np.allclose + :return: True if the matrix is Hermitian; False otherwise. + """ + if not is_square_matrix(matrix): + raise ValueError("The matrix is not square.") + return np.allclose(matrix, matrix.T.conj(), rtol=rtol, atol=atol) + + +def is_unitary_matrix(matrix: np.ndarray, rtol: float = 1e-05, atol: float = 1e-08) -> bool: + r""" + Checks if a square matrix A is unitary. That is + + A is unitary iff $A^{\dagger} A = A A^{\dagger}$ = Id, + + where $\dagger$ denotes conjugate transpose and Id denotes the identity. + + :param matrix: a M by M matrix. + :param rtol: The relative tolerance parameter in np.allclose + :param atol: The absolute tolerance parameter in np.allclose + :return: True if the matrix is normal; False otherwise. + """ + if not is_square_matrix(matrix): + raise ValueError("The matrix is not square.") + AB = matrix.T.conj() @ matrix + BA = matrix @ matrix.T.conj() + Id = np.eye(len(matrix)) + return np.allclose(AB, Id, rtol=rtol, atol=atol) and np.allclose(BA, Id, rtol=rtol, atol=atol) + + +def is_positive_definite_matrix(matrix: np.ndarray, + rtol: float = 1e-05, + atol: float = 1e-08) -> bool: + r""" + Checks if a square Hermitian matrix A is positive definite. That is + + A is positive definite iff eig(A) > 0. + + In this numerical implementation we check if each eigenvalue obeys eig(A) > -|atol|, + the strict condition can be recoved by setting `atol = 0`. + + :param matrix: a M by M Hermitian matrix. + :param rtol: The relative tolerance parameter in np.allclose + :param atol: The absolute tolerance parameter in np.allclose + :return: True if the matrix is normal; False otherwise. + """ + if not is_hermitian_matrix(matrix, rtol, atol): + raise ValueError("The matrix is not Hermitian.") + evals, _ = np.linalg.eigh(matrix) + return all(x > -abs(atol) for x in evals) + + +def is_positive_semidefinite_matrix(matrix: np.ndarray, + rtol: float = 1e-05, + atol: float = 1e-08) -> bool: + r""" + Checks if a square Hermitian matrix A is positive semi-definite. That is + + A is positive semi-definite iff eig(A) >= 0. + + In this numerical implementation we check if each eigenvalue obeys + + eig(A) >= -|atol|. + + :param matrix: a M by M Hermitian matrix. + :param rtol: The relative tolerance parameter in np.allclose + :param atol: The absolute tolerance parameter in np.allclose + :return: True if the matrix is normal; False otherwise. + """ + if not is_hermitian_matrix(matrix, rtol, atol): + raise ValueError("The matrix is not Hermitian.") + evals, _ = np.linalg.eigh(matrix) + return all(x >= -abs(atol) for x in evals) diff --git a/forest/benchmarking/operator_tools/validate_superoperator.py b/forest/benchmarking/operator_tools/validate_superoperator.py new file mode 100644 index 00000000..00c8ce6e --- /dev/null +++ b/forest/benchmarking/operator_tools/validate_superoperator.py @@ -0,0 +1,155 @@ +"""A module allowing one to check if superoperators or channels are physical. + +We have arbitrarily decided to use a column stacking convention. + +For more information about the conventions used, look at the file in +/docs/Superoperator representations.md + +Further references include: + +[GRAPTN] Tensor networks and graphical calculus for open quantum systems + Wood et al. + Quant. Inf. Comp. 15, 0579-0811 (2015) + (no DOI) + https://arxiv.org/abs/1111.6950 + +[MATQO] On the Matrix Representation of Quantum Operations + Nambu et al., + arXiv: 0504091 (2005) + https://arxiv.org/abs/quant-ph/0504091 + +[DUAL] On duality between quantum maps and quantum states + Zyczkowski et al., + Open Syst. Inf. Dyn. 11, 3 (2004) + https://dx.doi.org/10.1023/B:OPSY.0000024753.05661.c2 + https://arxiv.org/abs/quant-ph/0401119 + +""" +from typing import Sequence +import numpy as np +from forest.benchmarking.utils import partial_trace +from forest.benchmarking.operator_tools.superoperator_transformations import choi2kraus +from forest.benchmarking.operator_tools.validate_operator import ( + is_positive_semidefinite_matrix, is_identity_matrix, is_hermitian_matrix) +from forest.benchmarking.operator_tools.apply_superoperator import apply_choi_matrix_2_state + + +# ================================================================================================== +# Check physicality of Channels +# ================================================================================================== +def kraus_operators_are_valid(kraus_ops: Sequence[np.ndarray], + rtol: float = 1e-05, + atol: float = 1e-08) -> bool: + """ + Checks if a set of Kraus operators are valid. + + :param kraus_ops: A list of Kraus operators + :param rtol: The relative tolerance parameter in np.allclose + :param atol: The absolute tolerance parameter in np.allclose + :return: True if the Kraus operators are valid with the given tolerance; False otherwise. + """ + if isinstance(kraus_ops, np.ndarray): # handle input of single kraus op + if len(kraus_ops[0].shape) < 2: # first elem is not a matrix + kraus_ops = [kraus_ops] + rows, _ = np.asarray(kraus_ops[0]).shape + # Standard case of square Kraus operators is if rows==cols. For non-square Kraus ops it is + # required that sum_i M_i^\dagger M_i = np.eye(rows,rows). + POVM = [np.transpose(op).conjugate().dot(op) for op in kraus_ops] + all_POVM_elements_are_PSD = all(is_positive_semidefinite_matrix(elem) for elem in POVM) + POVM_elements_sum_to_Id = is_identity_matrix(sum(POVM), rtol, atol) + return all_POVM_elements_are_PSD and POVM_elements_sum_to_Id + + +def choi_is_hermitian_preserving(choi: np.ndarray, rtol: float = 1e-05, + atol: float = 1e-08) -> bool: + """ + Checks if a quantum process, specified by a Choi matrix, is hermitian-preserving. + + :param choi: a dim**2 by dim**2 Choi matrix + :param rtol: The relative tolerance parameter in np.allclose + :param atol: The absolute tolerance parameter in np.allclose + :return: Returns True if the quantum channel is hermitian preserving with the given tolerance; + False otherwise. + """ + # Equation 3.31 of [GRAPTN] + return is_hermitian_matrix(choi, rtol, atol) + + +def choi_is_trace_preserving(choi: np.ndarray, rtol: float = 1e-05, atol: float = 1e-08) -> bool: + """ + Checks if a quantum process, specified by a Choi matrix, is trace-preserving (TP). + + :param choi: A dim**2 by dim**2 Choi matrix + :param rtol: The relative tolerance parameter in np.allclose + :param atol: The absolute tolerance parameter in np.allclose + :return: Returns True if the quantum channel is trace-preserving with the given tolerance; + False otherwise. + """ + rows, cols = choi.shape + dim = int(np.sqrt(rows)) + # the choi matrix acts on the Hilbert space H_{in} \otimes H_{out}. + # We want to trace out H_{out} and so keep the H_{in} space at index 0. + keep = [0] + id_iff_tp = partial_trace(choi, keep, [dim, dim]) + # Equation 3.33 of [GRAPTN] + return is_identity_matrix(id_iff_tp, rtol, atol) + + +def choi_is_completely_positive(choi: np.ndarray, rtol: float = 1e-05, atol: float = 1e-08) -> bool: + """ + Checks if a quantum process, specified by a Choi matrix, is completely positive (CP). + + :param choi: A dim**2 by dim**2 Choi matrix + :param rtol: The relative tolerance parameter in np.allclose + :param atol: The absolute tolerance parameter in np.allclose + :return: Returns True if the quantum channel is completely positive with the given tolerance; + False otherwise. + """ + # Equation 3.35 of [GRAPTN] + return is_positive_semidefinite_matrix(choi, rtol, atol) + + +def choi_is_cptp(choi: np.ndarray, rtol: float = 1e-05, atol: float = 1e-08) -> bool: + """ + Checks if a quantum process, specified by a Choi matrix, is completely positive and + trace-preserving (CPTP). + + :param choi: A dim**2 by dim**2 Choi matrix + :param rtol: The relative tolerance parameter in np.allclose + :param atol: The absolute tolerance parameter in np.allclose + :return: Returns True if the quantum channel is CPTP with the given tolerance; False otherwise. + """ + choi_is_TP = choi_is_trace_preserving(choi, rtol, atol) + choi_is_CP = choi_is_completely_positive(choi, rtol, atol) + return choi_is_CP and choi_is_TP + + +def choi_is_unital(choi: np.ndarray, rtol: float = 1e-05, atol: float = 1e-08) -> bool: + """ + Checks if a quantum process, specified by a Choi matrix, is unital. + + A process is unital iff it maps the identity to itself. + + :param choi: A dim**2 by dim**2 Choi matrix + :param rtol: The relative tolerance parameter in np.allclose + :param atol: The absolute tolerance parameter in np.allclose + :return: Returns True if the quantum channel is unital with the given tolerance; False + otherwise. + """ + rows, cols = choi.shape + dim = int(np.sqrt(rows)) + id_iff_unital = apply_choi_matrix_2_state(choi, np.identity(dim)) + return is_identity_matrix(id_iff_unital, rtol, atol) + + +def choi_is_unitary(choi: np.ndarray, limit: float = 1e-09) -> bool: + """ + Checks if a quantum process, specified by a Choi matrix, is unitary. + + :param choi: A dim**2 by dim**2 Choi matrix + :param limit: A tolerance parameter to discard Kraus operators with small norm. + :return: Returns True if the quantum channel is unitary with the given tolerance; False + otherwise. + """ + kraus_ops = choi2kraus(choi, tol=limit) + return len(kraus_ops) == 1 diff --git a/forest/benchmarking/plotting/state_process.py b/forest/benchmarking/plotting/state_process.py index 3cb173fb..ff5725fd 100644 --- a/forest/benchmarking/plotting/state_process.py +++ b/forest/benchmarking/plotting/state_process.py @@ -110,11 +110,12 @@ def plot_pauli_transfer_matrix(ptransfermatrix, ax, labels=None, title='', fonts ticklabs = cb.ax.get_yticklabels() cb.ax.set_yticklabels(ticklabs, ha='right') cb.ax.yaxis.set_tick_params(pad=35) + cb.draw_all() ax.set_xticks(range(dim_squared)) ax.set_xlabel("Input Pauli Operator", fontsize=fontsizes) ax.set_yticks(range(dim_squared)) ax.set_ylabel("Output Pauli Operator", fontsize=fontsizes) - ax.set_title(title, fontsize= int(np.floor(1.2*fontsizes))) + ax.set_title(title, fontsize= int(np.floor(1.2*fontsizes)), pad=15) ax.set_xticklabels(labels, rotation=45, fontsize=int(np.floor(0.7*fontsizes))) ax.set_yticklabels(labels, fontsize=int(np.floor(0.7*fontsizes))) ax.grid(False) diff --git a/forest/benchmarking/tests/test_apply_superoperator.py b/forest/benchmarking/tests/test_apply_superoperator.py new file mode 100644 index 00000000..430c72a9 --- /dev/null +++ b/forest/benchmarking/tests/test_apply_superoperator.py @@ -0,0 +1,29 @@ +import numpy as np +from forest.benchmarking.tests.test_superoperator_transformations import ( + amplitude_damping_kraus, amplitude_damping_choi, ONE_STATE, ZERO_STATE, rho_out) +from forest.benchmarking.operator_tools.apply_superoperator import (apply_kraus_ops_2_state, + apply_choi_matrix_2_state) + + +def test_apply_kraus_ops_2_state(): + AD_kraus = amplitude_damping_kraus(0.1) + # rho_out was calculated by hand + assert np.allclose(rho_out, apply_kraus_ops_2_state(AD_kraus, ONE_STATE)) + + +def test_apply_non_square_kraus_ops_2_state(): + Id = np.eye(2) + bra_zero = np.asarray([[1], [0]]) + bra_one = np.asarray([[0], [1]]) + state_one = np.kron(Id / 2, ONE_STATE) + state_zero = np.kron(Id / 2, ZERO_STATE) + Kraus1 = np.kron(Id, bra_one.transpose()) + Kraus0 = np.kron(Id, bra_zero.transpose()) + assert np.allclose(apply_kraus_ops_2_state(Kraus0, state_zero), Id / 2) + assert np.allclose(apply_kraus_ops_2_state(Kraus1, state_one), Id / 2) + assert np.allclose(apply_kraus_ops_2_state(Kraus0, state_one), 0) + + +def test_apply_choi_matrix_2_state(): + choi = amplitude_damping_choi(0.1) + assert np.allclose(rho_out, apply_choi_matrix_2_state(choi, ONE_STATE)) diff --git a/forest/benchmarking/tests/test_channel_approximation.py b/forest/benchmarking/tests/test_channel_approximation.py new file mode 100644 index 00000000..ac5f2e79 --- /dev/null +++ b/forest/benchmarking/tests/test_channel_approximation.py @@ -0,0 +1,34 @@ +import numpy as np +from forest.benchmarking.operator_tools.channel_approximation import pauli_twirl_chi_matrix +from forest.benchmarking.tests.test_superoperator_transformations import (one_q_pauli_channel_chi, + amplitude_damping_chi) + + +# Pauli twirled Amplitude damping channel +def analytical_pauli_twirl_of_AD_chi(p): + # see equation 7 of https://arxiv.org/pdf/1701.03708.pdf + poly1 = (2 + 2 * np.sqrt(1 - p) - p) / 4 + poly2 = p / 4 + poly3 = (2 - 2 * np.sqrt(1 - p) - p) / 4 + pp_chi = np.asarray([[poly1, 0, 0, 0], + [0, poly2, 0, 0], + [0, 0, poly2, 0], + [0, 0, 0, poly3]]) + return pp_chi + + +def test_pauli_twirl_of_pauli_channel(): + # diagonal channel so should not change anything + px = np.random.rand() + py = np.random.rand() + pz = np.random.rand() + pauli_chan_chi_matrix = one_q_pauli_channel_chi(px, py, pz) + pauli_twirled_chi_matrix = pauli_twirl_chi_matrix(pauli_chan_chi_matrix) + assert np.allclose(pauli_chan_chi_matrix, pauli_twirled_chi_matrix) + + +def test_pauli_twirl_of_amp_damp(): + p = np.random.rand() + ana_chi = analytical_pauli_twirl_of_AD_chi(p) + num_chi = pauli_twirl_chi_matrix(amplitude_damping_chi(p)) + assert np.allclose(ana_chi, num_chi) diff --git a/forest/benchmarking/tests/test_compose_superoperators.py b/forest/benchmarking/tests/test_compose_superoperators.py new file mode 100644 index 00000000..0f53aa4a --- /dev/null +++ b/forest/benchmarking/tests/test_compose_superoperators.py @@ -0,0 +1,32 @@ +import numpy as np +from pyquil.gate_matrices import I, X, H +from forest.benchmarking.tests.test_superoperator_transformations import amplitude_damping_kraus +from forest.benchmarking.operator_tools.superoperator_transformations import kraus2superop +from forest.benchmarking.operator_tools.compose_superoperators import (compose_channel_kraus, + tensor_channel_kraus) + + +def bit_flip_kraus(p): + M0 = np.sqrt(1 - p) * I + M1 = np.sqrt(p) * X + return [M0, M1] + + +AD_kraus = amplitude_damping_kraus(0.1) +BitFlip_kraus = bit_flip_kraus(0.2) +BitFlip_super = kraus2superop(BitFlip_kraus) +AD_super = kraus2superop(AD_kraus) + + +def test_compose_channel_kraus(): + function_output = kraus2superop(compose_channel_kraus(AD_kraus, BitFlip_kraus)) + independent_answer = AD_super @ BitFlip_super + assert np.allclose(function_output, independent_answer) + + +def test_tensor_channel_kraus(): + function_output = tensor_channel_kraus([X], [H]) + independent_answer = np.kron(X,H) + wrong_answer = np.kron(H,X) + assert np.allclose(function_output, independent_answer) + assert not np.allclose(function_output, wrong_answer) diff --git a/forest/benchmarking/tests/test_project_state_matrix.py b/forest/benchmarking/tests/test_project_state_matrix.py new file mode 100644 index 00000000..5d6a7273 --- /dev/null +++ b/forest/benchmarking/tests/test_project_state_matrix.py @@ -0,0 +1,14 @@ +import numpy as np +from forest.benchmarking.operator_tools.project_state_matrix import project_state_matrix_to_physical + + +def test_project_state_matrix(): + """ + Test the wizard method. Example from fig 1 of maximum likelihood minimum effort + https://doi.org/10.1103/PhysRevLett.108.070502 + + :return: + """ + eigs = np.diag(np.array(list(reversed([3.0 / 5, 1.0 / 2, 7.0 / 20, 1.0 / 10, -11.0 / 20])))) + phys = project_state_matrix_to_physical(eigs) + assert np.allclose(phys, np.diag([0, 0, 1.0 / 5, 7.0 / 20, 9.0 / 20])) diff --git a/forest/benchmarking/tests/test_project_superoperators.py b/forest/benchmarking/tests/test_project_superoperators.py new file mode 100644 index 00000000..a4d6b809 --- /dev/null +++ b/forest/benchmarking/tests/test_project_superoperators.py @@ -0,0 +1,92 @@ +import numpy as np +from pyquil.gate_matrices import X, Y, Z, H +from forest.benchmarking.operator_tools.superoperator_transformations import kraus2choi +from forest.benchmarking.operator_tools.validate_superoperator import * +from forest.benchmarking.operator_tools.project_superoperators import * + + +def test_proj_to_cp(): + state = np.eye(2) + assert np.allclose(state, proj_choi_to_completely_positive(state)) + + state = np.array([[1.5, 0], [0, 10]]) + assert choi_is_completely_positive(state) + assert np.allclose(state, proj_choi_to_completely_positive(state)) + + state = -Z + cp_state = proj_choi_to_completely_positive(state) + target = np.array([[0, 0], [0, 1.]]) + assert choi_is_completely_positive(cp_state) + assert np.allclose(target, cp_state) + + state = X + cp_state = proj_choi_to_completely_positive(state) + target = np.array([[.5, .5], [.5, .5]]) + assert choi_is_completely_positive(cp_state) + assert np.allclose(target, cp_state) + + state = Y + cp_state = proj_choi_to_completely_positive(state) + target = np.array([[.5, -.5j], [.5j, .5]]) + assert choi_is_completely_positive(cp_state) + assert np.allclose(target, cp_state) + + choi = kraus2choi(np.kron(X, Z)) + assert choi_is_completely_positive(choi) + assert np.allclose(choi, proj_choi_to_completely_positive(choi)) + + +def test_proj_to_tp(): + # Identity process is trace preserving, so no change + choi = kraus2choi(np.eye(2)) + assert np.allclose(choi, proj_choi_to_trace_preserving(choi)) + + # Bit flip process is trace preserving, so no change + choi = kraus2choi(X) + assert np.allclose(choi, proj_choi_to_trace_preserving(choi)) + + # start with a non-trace-preserving choi. + choi = kraus2choi(X - np.eye(2)*.01) + assert choi_is_trace_preserving(proj_choi_to_trace_preserving(choi)) + + +def test_proj_to_tni(): + # Bit flip process is trace preserving, so no change + choi = kraus2choi(X) + assert np.allclose(choi, proj_choi_to_trace_non_increasing(choi)) + + choi = np.array( + [[0., 0., 0., 0.], [0., 1.01, 1.01, 0.], [0., 1., 1., 0.], [0., 0., 0., 0.]]) + assert choi_is_trace_preserving(proj_choi_to_trace_non_increasing(choi)) + + # start with a non-trace-preserving choi. + choi = kraus2choi(np.kron(X - np.eye(2) * .01, np.eye(2))) + choi_tni = proj_choi_to_trace_non_increasing(choi) + plusplus = np.array([[1, 1, 1, 1]]).T / 2 + rho_pp = plusplus @ plusplus.T + output = apply_choi_matrix_2_state(choi_tni, rho_pp) + assert 0 < np.trace(output) <= 1 + + +def test_proj_to_cptp(): + # Identity process is cptp, so no change + choi = kraus2choi(np.eye(2)) + assert np.allclose(choi, proj_choi_to_physical(choi)) + + # Bit flip process is cptp, so no change + choi = kraus2choi(X) + assert np.allclose(choi, proj_choi_to_physical(choi)) + + # Small perturbation shouldn't change too much + choi = np.array([[1.001, 0., 0., .99], [0., 0., 0., 0.], [0., 0., 0., 0.], + [1.004, 0., 0., 1.01]]) + assert np.allclose(choi, proj_choi_to_physical(choi), atol=1e-2) + + # Ensure final product is cptp with arbitrary perturbation + choi = np.array([[1.1, 0.2, -0.4, .9], [.5, 0., 0., 0.], [0., 0., 0., 0.], + [1.4, 0., 0., .8]]) + physical_choi = proj_choi_to_physical(choi) + assert choi_is_trace_preserving(physical_choi) + assert choi_is_completely_positive(physical_choi, atol=1e-1) + assert choi_is_cptp(physical_choi, atol=1e-1) + diff --git a/forest/benchmarking/tests/test_random_operators.py b/forest/benchmarking/tests/test_random_operators.py index c1b6bbcf..93091d3c 100644 --- a/forest/benchmarking/tests/test_random_operators.py +++ b/forest/benchmarking/tests/test_random_operators.py @@ -1,3 +1,4 @@ +import pytest import numpy.random numpy.random.seed(1) # seed random number generation for all calls to rand_ops @@ -5,9 +6,9 @@ from sympy.combinatorics import Permutation from numpy import linalg as la import forest.benchmarking.distance_measures as dm -import forest.benchmarking.random_operators as rand_ops -from forest.benchmarking.superoperator_tools import choi_is_trace_preserving, \ - choi_is_completely_positive +import forest.benchmarking.operator_tools.random_operators as rand_ops +from forest.benchmarking.operator_tools.validate_superoperator import ( + choi_is_trace_preserving,choi_is_completely_positive) D2_SWAP = np.array([[1, 0, 0, 0], [0, 0, 1, 0], @@ -29,6 +30,10 @@ # ================================================================================================= # Test: Permute tensor factors # ================================================================================================= +def test_permute_tensor_factor_id(): + np.testing.assert_array_equal(rand_ops.permute_tensor_factors(2, [0, 1, 2]), np.eye(8)) + + def test_permute_tensor_factor_SWAP(): # test the Dimension two and three SWAP operators D2 = 2 @@ -63,6 +68,17 @@ def test_permute_tensor_factor_three_qubits(): assert np.dot(fn_ans3.T, by_hand_ans3) == 1.0 +def test_permute_tensor_factor_diff_size_spaces(): + dims = [2, 4, 8] + perm = [2, 1, 0] + + alt_dims = 2 + alt_perm = [3, 4, 5, 1, 2, 0] + + np.testing.assert_array_equal(rand_ops.permute_tensor_factors(dims, perm), + rand_ops.permute_tensor_factors(alt_dims, alt_perm)) + + def test_permute_tensor_factor_four_qubit_permutation_operator(): # Test the generation of the permutation operator from the symbolic list # given to us by SymPy and if that is the same thing from Andrew Scott's papers. @@ -168,7 +184,7 @@ def test_random_unitaries_are_unitary(): assert np.allclose(avg_trace3, 3.0) assert np.allclose(avg_det3, 1.0) - +# ~ 11 sec; passed 2019/06/11 def test_random_unitaries_first_moment(): # the first moment should be proportional to P_21/D N_avg = 50000 @@ -203,6 +219,8 @@ def test_random_unitaries_first_moment(): # for dimensions 2 and 3 +# ~ 12 sec; passed 2019/06/11 +@pytest.mark.slow def test_random_unitaries_second_moment(): # the second moment should be proportional to # diff --git a/forest/benchmarking/tests/test_state_tomography.py b/forest/benchmarking/tests/test_state_tomography.py index 959ecaad..8849f02c 100644 --- a/forest/benchmarking/tests/test_state_tomography.py +++ b/forest/benchmarking/tests/test_state_tomography.py @@ -6,8 +6,7 @@ from forest.benchmarking.compilation import basic_compile from forest.benchmarking.random_operators import haar_rand_unitary from forest.benchmarking.tomography import generate_state_tomography_experiment, _R, \ - iterative_mle_state_estimate, project_density_matrix, estimate_variance, \ - linear_inv_state_estimate + iterative_mle_state_estimate, estimate_variance, linear_inv_state_estimate from pyquil.api import ForestConnection, QuantumComputer, QVM from pyquil.api._compiler import _extract_attribute_dictionary_from_program from pyquil.api._qac import AbstractCompiler @@ -263,18 +262,6 @@ def test_hedged_two_qubit(two_q_tomo_fixture): np.testing.assert_allclose(rho_true, rho_est, atol=0.02) -def test_project_density_matrix(): - """ - Test the wizard method. Example from fig 1 of maximum likelihood minimum effort - https://doi.org/10.1103/PhysRevLett.108.070502 - - :return: - """ - eigs = np.diag(np.array(list(reversed([3.0 / 5, 1.0 / 2, 7.0 / 20, 1.0 / 10, -11.0 / 20])))) - phys = project_density_matrix(eigs) - assert np.allclose(phys, np.diag([0, 0, 1.0 / 5, 7.0 / 20, 9.0 / 20])) - - def test_variance_bootstrap(two_q_tomo_fixture): qubits = [0, 1] results, rho_true = two_q_tomo_fixture diff --git a/forest/benchmarking/tests/test_superoperator_tools.py b/forest/benchmarking/tests/test_superoperator_transformations.py similarity index 57% rename from forest/benchmarking/tests/test_superoperator_tools.py rename to forest/benchmarking/tests/test_superoperator_transformations.py index 4e59a316..65757516 100644 --- a/forest/benchmarking/tests/test_superoperator_tools.py +++ b/forest/benchmarking/tests/test_superoperator_transformations.py @@ -1,7 +1,6 @@ import numpy as np from pyquil.gate_matrices import X, Y, Z, H from forest.benchmarking.superoperator_tools import * -import forest.benchmarking.random_operators as rand_ops # Test philosophy: @@ -115,10 +114,6 @@ def analytical_pauli_twirl_of_AD_chi(p): np.matmul(np.matmul(AdKrausOps[1], ONE_STATE), AdKrausOps[1].transpose().conj()) -# =================================================================================================== -# Test superoperator conversion tools -# =================================================================================================== - def test_vec(): A = np.asarray([[1, 2], [3, 4]]) B = np.asarray([[1, 2, 5], [3, 4, 6]]) @@ -290,196 +285,3 @@ def test_choi_pl_bijectivity(): h_superop = kraus2superop(H) assert np.allclose(choi2superop(choi2superop(h_choi)), h_choi) assert np.allclose(superop2choi(superop2choi(h_superop)), h_superop) - - -# =================================================================================================== -# Test channel and superoperator approximation tools -# =================================================================================================== - - -def test_pauli_twirl_of_pauli_channel(): - # diagonal channel so should not change anything - px = np.random.rand() - py = np.random.rand() - pz = np.random.rand() - pauli_chan_chi_matrix = one_q_pauli_channel_chi(px, py, pz) - pauli_twirled_chi_matrix = pauli_twirl_chi_matrix(pauli_chan_chi_matrix) - assert np.allclose(pauli_chan_chi_matrix, pauli_twirled_chi_matrix) - - -def test_pauli_twirl_of_amp_damp(): - p = np.random.rand() - ana_chi = analytical_pauli_twirl_of_AD_chi(p) - num_chi = pauli_twirl_chi_matrix(amplitude_damping_chi(p)) - assert np.allclose(ana_chi, num_chi) - - -# ================================================================================================== -# Test apply channel -# ================================================================================================== - - -def test_apply_kraus_ops_2_state(): - AD_kraus = amplitude_damping_kraus(0.1) - assert np.allclose(rho_out, apply_kraus_ops_2_state(AD_kraus, ONE_STATE)) - - -def test_apply_non_square_kraus_ops_2_state(): - Id = np.eye(2) - bra_zero = np.asarray([[1], [0]]) - bra_one = np.asarray([[0], [1]]) - state_one = np.kron(Id / 2, ONE_STATE) - state_zero = np.kron(Id / 2, ZERO_STATE) - Kraus1 = np.kron(Id, bra_one.transpose()) - Kraus0 = np.kron(Id, bra_zero.transpose()) - assert np.allclose(apply_kraus_ops_2_state(Kraus0, state_zero), Id / 2) - assert np.allclose(apply_kraus_ops_2_state(Kraus1, state_one), Id / 2) - assert np.allclose(apply_kraus_ops_2_state(Kraus0, state_one), 0) - - -def test_apply_choi_matrix_2_state(): - choi = amplitude_damping_choi(0.1) - assert np.allclose(rho_out, apply_choi_matrix_2_state(choi, ONE_STATE)) - - -# ================================================================================================== -# Test physicality of Channels -# ================================================================================================== -def test_kraus_operators_are_valid(): - assert kraus_operators_are_valid(amplitude_damping_kraus(np.random.rand())) - assert kraus_operators_are_valid(H) - assert not kraus_operators_are_valid(AdKrausOps[0]) - -def test_choi_is_hermitian_preserving(): - D = 2 - K = 2 - choi = rand_ops.rand_map_with_BCSZ_dist(D, K) - assert choi_is_hermitian_preserving(choi) - - -def test_choi_is_trace_preserving(): - D = 2 - K = 2 - choi = rand_ops.rand_map_with_BCSZ_dist(D, K) - assert choi_is_trace_preserving(choi) - - -def test_choi_is_completely_positive(): - D = 2 - K = 2 - choi = rand_ops.rand_map_with_BCSZ_dist(D, K) - assert choi_is_completely_positive(choi) - D = 3 - K = 2 - choi = rand_ops.rand_map_with_BCSZ_dist(D, K) - assert choi_is_completely_positive(choi) - - -def test_choi_is_unital(): - px = np.random.rand() - py = np.random.rand() - pz = np.random.rand() - choi = chi2choi(one_q_pauli_channel_chi(px, py, pz)) - assert choi_is_unital(choi) - assert choi_is_unital(HADChoi) - assert not choi_is_unital(amplitude_damping_choi(0.1)) - - -def test_choi_is_unitary(): - px = np.random.rand() - py = np.random.rand() - pz = np.random.rand() - choi = chi2choi(one_q_pauli_channel_chi(px, py, pz)) - assert not choi_is_unitary(choi) - assert choi_is_unital(choi) - assert choi_is_unitary(HADChoi) - assert not choi_is_unitary(amplitude_damping_choi(0.1)) - -# ================================================================================================== -# Test project Channels to CP, TNI, TP, and physical -# ================================================================================================== - -def test_proj_to_cp(): - state = np.eye(2) - assert np.allclose(state, proj_choi_to_completely_positive(state)) - - state = np.array([[1.5, 0], [0, 10]]) - assert choi_is_completely_positive(state) - assert np.allclose(state, proj_choi_to_completely_positive(state)) - - state = -Z - cp_state = proj_choi_to_completely_positive(state) - target = np.array([[0, 0], [0, 1.]]) - assert choi_is_completely_positive(cp_state) - assert np.allclose(target, cp_state) - - state = X - cp_state = proj_choi_to_completely_positive(state) - target = np.array([[.5, .5], [.5, .5]]) - assert choi_is_completely_positive(cp_state) - assert np.allclose(target, cp_state) - - state = Y - cp_state = proj_choi_to_completely_positive(state) - target = np.array([[.5, -.5j], [.5j, .5]]) - assert choi_is_completely_positive(cp_state) - assert np.allclose(target, cp_state) - - choi = kraus2choi(np.kron(X, Z)) - assert choi_is_completely_positive(choi) - assert np.allclose(choi, proj_choi_to_completely_positive(choi)) - - -def test_proj_to_tp(): - # Identity process is trace preserving, so no change - choi = kraus2choi(np.eye(2)) - assert np.allclose(choi, proj_choi_to_trace_preserving(choi)) - - # Bit flip process is trace preserving, so no change - choi = kraus2choi(X) - assert np.allclose(choi, proj_choi_to_trace_preserving(choi)) - - # start with a non-trace-preserving choi. - choi = kraus2choi(X - np.eye(2)*.01) - assert choi_is_trace_preserving(proj_choi_to_trace_preserving(choi)) - - -def test_proj_to_tni(): - # Bit flip process is trace preserving, so no change - choi = kraus2choi(X) - assert np.allclose(choi, proj_choi_to_trace_non_increasing(choi)) - - choi = np.array( - [[0., 0., 0., 0.], [0., 1.01, 1.01, 0.], [0., 1., 1., 0.], [0., 0., 0., 0.]]) - assert choi_is_trace_preserving(proj_choi_to_trace_non_increasing(choi)) - - # start with a non-trace-preserving choi. - choi = kraus2choi(np.kron(X - np.eye(2) * .01, np.eye(2))) - choi_tni = proj_choi_to_trace_non_increasing(choi) - plusplus = np.array([[1, 1, 1, 1]]).T / 2 - rho_pp = plusplus @ plusplus.T - output = apply_choi_matrix_2_state(choi_tni, rho_pp) - assert 0 < np.trace(output) <= 1 - - -def test_proj_to_cptp(): - # Identity process is cptp, so no change - choi = kraus2choi(np.eye(2)) - assert np.allclose(choi, proj_choi_to_physical(choi)) - - # Bit flip process is cptp, so no change - choi = kraus2choi(X) - assert np.allclose(choi, proj_choi_to_physical(choi)) - - # Small perturbation shouldn't change too much - choi = np.array([[1.001, 0., 0., .99], [0., 0., 0., 0.], [0., 0., 0., 0.], - [1.004, 0., 0., 1.01]]) - assert np.allclose(choi, proj_choi_to_physical(choi), atol=.01) - - # Ensure final product is cptp with arbitrary perturbation - choi = np.array([[1.1, 0.2, -0.4, .9], [.5, 0., 0., 0.], [0., 0., 0., 0.], - [1.4, 0., 0., .8]]) - physical_choi = proj_choi_to_physical(choi) - assert choi_is_trace_preserving(physical_choi) - assert choi_is_completely_positive(physical_choi, limit=1e-1) - diff --git a/forest/benchmarking/tests/test_validate_operators.py b/forest/benchmarking/tests/test_validate_operators.py new file mode 100644 index 00000000..837fbac4 --- /dev/null +++ b/forest/benchmarking/tests/test_validate_operators.py @@ -0,0 +1,77 @@ +import pytest +import numpy as np +from pyquil.gate_matrices import X, Y, Z, H +from forest.benchmarking.operator_tools.random_operators import haar_rand_unitary +from forest.benchmarking.operator_tools.validate_operator import * + + +# Matrix below is from https://en.wikipedia.org/wiki/Normal_matrix +# it is normal but NOT unitary or Hermitian +NORMAL = np.array([[1, 1, 0], [0, 1, 1], [1, 0, 1]]) + +# not symmetric +SIGMA_MINUS = (X + 1j*Y)/2 # np.array([[0, 1], [0, 0]]) + +# idempotent +PROJ_ZERO = np.array([[1, 0], [0, 0]]) + + +def test_is_square_matrix(): + assert is_square_matrix(np.eye(3)) + with pytest.raises(ValueError): + is_square_matrix(np.ndarray(shape=(2, 2, 2))) + assert not is_square_matrix(np.array([[1, 0]])) + + +def test_is_symmetric_matrix(): + assert is_symmetric_matrix(X) + assert not is_symmetric_matrix(SIGMA_MINUS) + with pytest.raises(ValueError): + is_symmetric_matrix(np.ndarray(shape=(2, 2, 2))) + with pytest.raises(ValueError): + is_symmetric_matrix(np.array([[1, 0]])) + + +def test_is_identity_matrix(): + assert not is_identity_matrix(Z) + assert is_identity_matrix(np.eye(3)) + with pytest.raises(ValueError): + is_identity_matrix(np.ndarray(shape=(2, 2, 2))) + with pytest.raises(ValueError): + is_identity_matrix(np.array([[1, 0]])) + + +def test_is_idempotent_matrix(): + assert not is_idempotent_matrix(SIGMA_MINUS) + assert is_idempotent_matrix(PROJ_ZERO) + assert is_idempotent_matrix(np.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]])) + + +def test_is_normal_matrix(): + assert is_normal_matrix(NORMAL) + assert not is_normal_matrix(SIGMA_MINUS) + + +def test_is_hermitian_matrix(): + assert not is_hermitian_matrix(NORMAL) + assert is_hermitian_matrix(X) + assert is_hermitian_matrix(Y) + + +def test_is_unitary_matrix(): + assert not is_unitary_matrix(NORMAL) + assert is_unitary_matrix(Y) + assert is_unitary_matrix(haar_rand_unitary(4)) + + +def test_is_positive_definite_matrix(): + # not atol is = 1e-08 + assert not is_positive_definite_matrix(np.array([[-1e-08, 0], [0, 0.1]])) + assert is_positive_definite_matrix(np.array([[0.5e-08, 0], [0, 0.1]])) + + +def test_is_positive_semidefinite_matrix(): + # not atol is = 1e-08 + assert not is_positive_semidefinite_matrix(np.array([[-1e-07, 0], [0, 0.1]])) + assert is_positive_semidefinite_matrix(np.array([[-1e-08, 0], [0, 0.1]])) + assert is_positive_semidefinite_matrix(np.array([[0.5e-08, 0], [0, 0.1]])) diff --git a/forest/benchmarking/tests/test_validate_superoperator.py b/forest/benchmarking/tests/test_validate_superoperator.py new file mode 100644 index 00000000..e209fe78 --- /dev/null +++ b/forest/benchmarking/tests/test_validate_superoperator.py @@ -0,0 +1,60 @@ +from pyquil.gate_matrices import H +import forest.benchmarking.operator_tools.random_operators as rand_ops +from forest.benchmarking.operator_tools.superoperator_transformations import chi2choi +from forest.benchmarking.tests.test_superoperator_transformations import ( + amplitude_damping_kraus, AdKrausOps, HADChoi, amplitude_damping_choi, one_q_pauli_channel_chi) +from forest.benchmarking.operator_tools.validate_superoperator import * + + +def test_kraus_operators_are_valid(): + assert kraus_operators_are_valid(amplitude_damping_kraus(np.random.rand())) + assert kraus_operators_are_valid(H) + assert not kraus_operators_are_valid(AdKrausOps[0]) + + +def test_choi_is_hermitian_preserving(): + D = 2 + K = 2 + choi = rand_ops.rand_map_with_BCSZ_dist(D, K) + assert choi_is_hermitian_preserving(choi) + + +def test_choi_is_trace_preserving(): + D = 2 + K = 2 + choi = rand_ops.rand_map_with_BCSZ_dist(D, K) + assert choi_is_trace_preserving(choi) + + +def test_choi_is_completely_positive(): + D = 2 + K = 2 + choi = rand_ops.rand_map_with_BCSZ_dist(D, K) + assert choi_is_completely_positive(choi) + D = 3 + K = 2 + choi = rand_ops.rand_map_with_BCSZ_dist(D, K) + assert choi_is_completely_positive(choi) + + +def test_choi_is_unital(): + px = np.random.rand() + py = np.random.rand() + pz = np.random.rand() + norm = px + py + pz + choi = chi2choi(one_q_pauli_channel_chi(px/norm, py/norm, pz/norm)) + assert choi_is_unital(choi) + assert choi_is_unital(HADChoi) + assert not choi_is_unital(amplitude_damping_choi(0.1)) + + +def test_choi_is_unitary(): + px = np.random.rand() + py = np.random.rand() + pz = np.random.rand() + norm = px + py + pz + choi = chi2choi(one_q_pauli_channel_chi(px/norm, py/norm, pz/norm)) + assert not choi_is_unitary(choi) + assert choi_is_unital(choi) + assert choi_is_unitary(HADChoi) + assert not choi_is_unitary(amplitude_damping_choi(0.1)) diff --git a/forest/benchmarking/tomography.py b/forest/benchmarking/tomography.py index cb06bf12..856ab0f3 100644 --- a/forest/benchmarking/tomography.py +++ b/forest/benchmarking/tomography.py @@ -5,7 +5,7 @@ import warnings import numpy as np -from scipy.linalg import logm, pinv, eigh +from scipy.linalg import logm, pinv from pyquil import Program from pyquil.unitary_tools import lifted_pauli as pauli2matrix, lifted_state_operator as state2matrix @@ -13,6 +13,7 @@ import forest.benchmarking.distance_measures as dm from forest.benchmarking.utils import all_traceless_pauli_terms from forest.benchmarking.superoperator_tools import vec, unvec, proj_choi_to_physical +from forest.benchmarking.operator_tools.project_state_matrix import project_state_matrix_to_physical from forest.benchmarking.observable_estimation import ExperimentSetting, ObservablesExperiment, \ ExperimentResult, SIC0, SIC1, SIC2, SIC3, plusX, minusX, plusY, minusY, plusZ, minusZ, \ TensorProductState, zeros_state @@ -352,56 +353,6 @@ def state_log_likelihood(state: np.ndarray, results: Iterator[ExperimentResult], return ll -def project_density_matrix(rho) -> np.ndarray: - """ - Project a possibly unphysical estimated density matrix to the closest (with respect to the - 2-norm) positive semi-definite matrix with trace 1, that is a valid quantum state. - - This is the so called "wizard" method. It is described in the following reference: - - [MLEWIZ] Efficient Method for Computing the Maximum-Likelihood Quantum State from - Measurements with Additive Gaussian Noise - Smolin et al., - Phys. Rev. Lett. 108, 070502 (2012) - https://doi.org/10.1103/PhysRevLett.108.070502 - https://arxiv.org/abs/1106.5458 - - :param rho: Numpy array containing the density matrix with dimension (N, N) - :return rho_projected: The closest positive semi-definite trace 1 matrix to rho. - """ - - # Rescale to trace 1 if the matrix is not already - rho_impure = rho / np.trace(rho) - - dimension = rho_impure.shape[0] # the dimension of the Hilbert space - [eigvals, eigvecs] = eigh(rho_impure) - - # If matrix is already trace one PSD, we are done - if np.min(eigvals) >= 0: - return rho_impure - - # Otherwise, continue finding closest trace one, PSD matrix - eigvals = list(eigvals) - eigvals.reverse() - eigvals_new = [0.0] * len(eigvals) - - i = dimension - accumulator = 0.0 # Accumulator - while eigvals[i - 1] + accumulator / float(i) < 0: - accumulator += eigvals[i - 1] - i -= 1 - for j in range(i): - eigvals_new[j] = eigvals[j] + accumulator / float(i) - eigvals_new.reverse() - - # Reconstruct the matrix - rho_projected = functools.reduce(np.dot, (eigvecs, - np.diag(eigvals_new), - np.conj(eigvecs.T))) - - return rho_projected - - def _resample_expectations_with_beta(results, prior_counts=1): """Resample expectation values by constructing a beta distribution and sampling from it. @@ -455,7 +406,7 @@ def estimate_variance(results: List[ExperimentResult], functional is measured. Not applicable if functional is ``dm.purity``. :param n_resamples: The number of times to resample. :param project_to_physical: Whether to project the estimated state to a physical one - with :py:func:`project_density_matrix`. + with :py:func:`project_state_matrix_to_physical`. """ if functional != dm.purity: if target_state is None: @@ -468,7 +419,7 @@ def estimate_variance(results: List[ExperimentResult], rho = tomo_estimator(resampled_results, qubits) if project_to_physical: - rho = project_density_matrix(rho) + rho = project_state_matrix_to_physical(rho) # Calculate functional of the state if functional == dm.purity: diff --git a/requirements.txt b/requirements.txt index 906b2298..e0f400ee 100644 --- a/requirements.txt +++ b/requirements.txt @@ -11,7 +11,7 @@ cvxpy>=1.0.0 dataclasses; python_version < "3.7" tqdm gitpython -matplotlib +matplotlib>=3.1.0 # test dependencies flake8