From 15348d530b0fdb84a2ed9928f4ecea85a4aa57cf Mon Sep 17 00:00:00 2001 From: AlexisRalli Date: Mon, 24 Oct 2022 09:53:34 +0100 Subject: [PATCH] projector notebook added and projector code slightly modified to improve speed --- notebooks/PauliwordOp projector.ipynb | 625 ++++++++++++++++++++++++++ symmer/symplectic/base.py | 22 +- 2 files changed, 637 insertions(+), 10 deletions(-) create mode 100644 notebooks/PauliwordOp projector.ipynb diff --git a/notebooks/PauliwordOp projector.ipynb b/notebooks/PauliwordOp projector.ipynb new file mode 100644 index 00000000..272e944e --- /dev/null +++ b/notebooks/PauliwordOp projector.ipynb @@ -0,0 +1,625 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "9923bd2a", + "metadata": {}, + "source": [ + "# Notes on Pauli projectors" + ] + }, + { + "cell_type": "markdown", + "id": "2872f745", + "metadata": {}, + "source": [ + "A Pauli projector on a **single qubit** can be defined as:\n", + "\n", + "\n", + "$$A = \\frac{1}{2}\\big( I \\pm \\sigma \\big)$$\n", + "\n", + "where $\\sigma \\in \\{X,Y,Z\\}$. The $+$ term above projects onto the $+1$ eigenvector and $-$ projects onto the $-1$ eigenvector" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ce0f610f", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "b2766f25", + "metadata": {}, + "source": [ + "Next consider projection of an $n$-qubits (requires taking products of single qubit projectors)\n", + "\n", + "$$B = \\prod_{i=0}^{n-1} \\frac{1}{2}\\Bigg[ \\bigg(\\bigotimes_{j=0}^{n-1} I \\bigg) \\pm \\bigg( \\bigotimes_{j=0}^{i-1} I \\otimes \\sigma_{i} \\otimes \\bigotimes_{j=i+1}^{n-1} I \\bigg) \\Bigg]$$\n", + "\n", + "$$= \\prod_{i=0}^{n-1} \\Bigg( \\frac{1}{2} \\big( II \\text{...}I_{i}\\text{...}I \\pm II\\text{...} \\sigma_{i} \\text{...} I\\big) \\Bigg)$$\n", + "\n", + "- first term is all identity string\n", + "- second term is all identities bar one unique single Pauli matrix\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8ec1d18b", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "ccccdda7", + "metadata": {}, + "source": [ + "In the computational basis $\\sigma = Z$ and we get:\n", + "\n", + "\n", + "$$C = \\prod_{i=0}^{n-1} \\Bigg( \\frac{1}{2} \\big( II \\text{...} I \\pm II\\text{...} Z_{i} \\text{...} I\\big) \\Bigg)$$\n", + "\n", + "This product result in a binary expansion with all combinations of $Z$ on each qubit:\n", + "\n", + "\n", + "$$C = \\frac{1}{2^{n}} \\sum_{i=0}^{2^{n}-1} \\pm P$$\n", + "\n", + "Which bitstring is chosen to project onto determines signs in the expansion...\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e96e265f", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "63af8040", + "metadata": {}, + "source": [ + "E.g. projection onto the state: $|10 \\rangle = |2 \\rangle$\n", + "\n", + "$$|2 \\rangle \\langle 2| = \\frac{1}{2^{2}}\\big( II - ZI \\big) \\big( II + IZ \\big)$$\n", + "$$ = \\frac{1}{4} \\big( II + IZ - ZI - ZZ \\big)$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9b18fe77", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "965d1c6d", + "metadata": {}, + "source": [ + "Let the set $\\mathcal{W}$ be all $2^{n}$ the binary bitstrings of size $n$, then we have: for projecting onto a particular computational basis state $| b \\rangle$:\n", + "\n", + "$$|b \\rangle \\langle b| = \\frac{1}{2^{n}} \\sum_{k \\in \\mathcal{W}}\\Bigg[ (-1)^{(b\\cdot k) \\: mod \\: 2}\\bigg( \\bigotimes_{i=0}^{n-1} \\sigma_{i}^{k[i]} \\bigg) \\Bigg]$$\n", + "\n", + "where:\n", + "- $k[i]$ is the $i$-th bit of the bitstring $| k \\rangle$ this determines each $\\sigma \\in \\{I, Z \\}$:\n", + " - $\\sigma_{i}^{0} = I_{i}$\n", + " - $\\sigma_{i}^{1} = Z_{i}$\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "e8e333e6", + "metadata": {}, + "source": [ + "This gives a **FAST** method to generate projectors in the symplectic picture" + ] + }, + { + "cell_type": "markdown", + "id": "bda1ca61", + "metadata": {}, + "source": [ + "To project onto the $X$ and $Y$ basis, the same process as above will work except one needs to replace the $Z_{i}$ operators with $X_{i}$ or $Y_{i}$\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a86745b4", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "09588529", + "metadata": {}, + "source": [ + "In general for a bitstring $| b^{'} \\rangle \\in \\{0,1,+,-, +i, -i \\}^{\\otimes n}$\n", + "\n", + "$$|b' \\rangle \\langle b'| = \\frac{1}{2^{n}} \\sum_{k \\in \\mathcal{W}}\\Bigg[ (-1)^{(b''\\cdot k) \\: mod \\: 2}\\bigg( \\bigotimes_{j=0}^{n-1} \\sigma_{j}^{k[j]} \\bigg) \\Bigg]$$\n", + "\n", + "where:\n", + "- $k[j]$ is the $j$-th bit of the bitstring $| k \\rangle$ this determines each $\\sigma \\in \\{I, X, Y, Z \\}$:\n", + " - $\\sigma_{j}^{0} = I_{i}$\n", + " - $\\sigma_{j}^{1} = Z_{i}$\n", + " - $\\sigma_{j}^{-} = X_{i}$\n", + " - $\\sigma_{j}^{+} = X_{i}$\n", + " - $\\sigma_{j}^{-i} = Y_{i}$\n", + " - $\\sigma_{j}^{+i} = Y_{i}$\n", + "- and $b''$ is the mapping of $b'$ to a $\\{0,1\\}^{\\otimes n}$ state\n", + " - $- \\mapsto 1$\n", + " - $+ \\mapsto 0$\n", + " - $-i \\mapsto 1$\n", + " - $+i \\mapsto 0$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a99dad2d", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "93db18af", + "metadata": {}, + "source": [ + "This is what `get_PauliwordOp_projector` does" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "9cae4f42", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + " 0.125+0.000j IIII +\n", + "-0.125+0.000j IIIZ +\n", + " 0.125+0.000j IXII +\n", + "-0.125+0.000j IXIZ +\n", + " 0.125+0.000j ZIII +\n", + "-0.125+0.000j ZIIZ +\n", + " 0.125+0.000j ZXII +\n", + "-0.125+0.000j ZXIZ" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from symmer.symplectic import get_PauliwordOp_projector\n", + "\n", + "# projector onto |0 + > ⊗ I ⊗ |1 >\n", + "Proj = get_PauliwordOp_projector('0+I1') # <--- can use 0,1, +,-, *, % for Z, X, Y basis and I to leave qubit unchanged\n", + "Proj" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "80e0b881", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "a9363223", + "metadata": {}, + "source": [ + "# FAST(er) decomposition of sparse matrix into Pauli operators" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ae3c406d", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "2eb20d05", + "metadata": {}, + "source": [ + "The projector $|i \\rangle \\langle j |$ defines an all zero matrix with a single one in the $i$-th row and $j$-th column.\n", + "\n", + "\n", + "Therefore a sparse matrix can be decomposed in an efficent manner (better than the $4^{n}$ trace approach that is)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1b0143ea", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "f3c0aedb", + "metadata": {}, + "source": [ + "Given the projector onto the zero state: \n", + "\n", + "$$|00...0 \\rangle \\langle 00...0| = \\frac{1}{2^{n}} \\sum_{k \\in \\mathcal{W}}\\Bigg[ \\bigg( \\bigotimes_{i=0}^{n-1} \\sigma_{i}^{k[i]} \\bigg) \\Bigg]$$\n", + "\n", + "$$= \\frac{1}{2^{n}} \\Big(II\\text{...}I + II\\text{...}Z +\\text{...} + ZZ\\text{...}Z \\Big) $$\n", + "\n", + "\n", + "- note this has $2^{n}$ terms!\n", + "\n", + "The **Heisenberg picture** lets us convert it into an arbitrary projector of $|i \\rangle \\langle j |$, but multiplying the appropriate operators on the left and right.\n", + "\n", + "Each multiplication takes $2^{n}$ (of which there are two)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73486e61", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "378efae3", + "metadata": {}, + "source": [ + "What is the cost of decomposing an arbitrary matrix $M$ using this strategy?\n", + "\n", + "Given the all zero projector.\n", + "Find all non-zero row and column pairs $i,j$ in the matrix $M$.\n", + "\n", + "1. for each $i,j$ convert the all zero projector to $|i \\rangle \\langle j |$ \n", + "2. multiply the new projector by $M[i,j]$\n", + "3. add to overall decomposition\n" + ] + }, + { + "cell_type": "markdown", + "id": "074b864e", + "metadata": {}, + "source": [ + "\n", + "How does this scale?\n", + "- step 1. requires two multiplications of a single operator over a $2^{n}$ size term\n", + "- step 2. is a multiplication of $2^{n}$ coefficients by a single number\n", + "- step 3. addition to tracking term (at worst the size of this object will become $4^{n})$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "82341cf6", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "69c08908", + "metadata": {}, + "source": [ + "\n", + "What is the computation cost?\n", + "\n", + "$$2^{N} \\times 2 \\times S_{M}$$\n", + "- two multiplications for every element in M (denoted $S$ for sparsity)\n", + "\n", + "$$ = \\mathcal{O}(S_{M} \\cdot 2^{n+1})$$\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0efb3a39", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "17c43561", + "metadata": {}, + "source": [ + "WHEREAS standard method does:\n", + "1. Build each $4^{n}$ pauli matrices\n", + "2. Multiply the matrix $M$ by each of the matrix for of each Pauli operators\n", + "3. take trace\n", + "4. each operation gives coeff weight\n", + "\n", + "A dense verion would require:\n", + "$$4^{n} \\times \\underbrace{\\big[ (2^{n})^{2} \\big]}_{\\text{mat mult}} \\times \\underbrace{2^{n}}_{\\text{trace}} = 2^{5n}$$\n", + "\n", + "Whereas a sparse form would need:\n", + "\n", + "$$4^{n} \\times \\underbrace{\\big[ S_{M}\\times S_{P}\\times 2^{n} \\big]}_{\\text{mat mul}} \\times \\underbrace{2^{n}}_{\\text{trace}} = 2^{4n} (S_{M}\\times S_{P})$$\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "642ce39a", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "7ea8035c", + "metadata": {}, + "source": [ + "Ignoring the sparsity of the each $4^{n}$ Pauli matrix, we find at best $\\mathcal{O}(2^{4n}\\cdot S_{M})$ \n", + "\n", + "The new approach improves upon this by:\n", + "\n", + "$$\\frac{(2^{4n}\\cdot S_{M})}{(S_{M} \\cdot 2^{n+1})} = 2^{3n-1}$$\n", + "\n", + "which is an **exponential speedup!**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a4bb8654", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "a108491e", + "metadata": {}, + "source": [ + "In the dense version:\n", + "\n", + "$$\\frac{(2^{5n})}{(S_{M} \\cdot 2^{n+1})} = \\frac{1}{S_{M}} (2^{4n-1})$$\n", + "\n", + "with a completely dense matrix:\n", + "\n", + "$$\\frac{(2^{5n})}{(2^{2n} \\cdot 2^{n+1})} = \\frac{(2^{5n})}{(2^{3n+1})} = 2^{2n-1}$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "babedbe1", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "6f80d638", + "metadata": {}, + "source": [ + "# Matrix to Pauli Code" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "b6770ca4", + "metadata": {}, + "outputs": [], + "source": [ + "from symmer.symplectic import PauliwordOp\n", + "import numpy as np\n", + "from scipy.sparse import rand" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "id": "ad4fa97b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(64, 64)" + ] + }, + "execution_count": 38, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "n_qubits = 6\n", + "\n", + "D = 2**n_qubits\n", + "x = rand(D, D, density=0.01, format='csr')\n", + "x.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "id": "51041632", + "metadata": {}, + "outputs": [], + "source": [ + "# n_qubits = 12\n", + "\n", + "# D = 2**n_qubits\n", + "# x = rand(D, D, density=0.00001, format='csr')\n", + "# x.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "id": "07fa9f59", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(41,)" + ] + }, + "execution_count": 40, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mat = x.toarray()\n", + "mat.nonzero()[0].shape" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "id": "5db9f421", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "31 ms ± 641 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" + ] + }, + { + "data": { + "text/plain": [ + "288" + ] + }, + "execution_count": 41, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# strategy == projector or full_basis\n", + "\n", + "%timeit Pop = PauliwordOp.from_matrix(mat, strategy='projector')\n", + "Pop.n_terms" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "id": "d7bac7bf", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.allclose(Pop.to_sparse_matrix.toarray(),\n", + " mat)" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "id": "aaa37721", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "1048576" + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# no chance:\n", + "4**n_qubits * (2**n_qubits)**2" + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "id": "59b62ce9", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2.67 s ± 9.08 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" + ] + } + ], + "source": [ + "%timeit Pop2 = PauliwordOp.from_matrix(mat, strategy='full_basis')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "483fe57f", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a679ac3b", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.7" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/symmer/symplectic/base.py b/symmer/symplectic/base.py index 77b868ff..5e1827e7 100644 --- a/symmer/symplectic/base.py +++ b/symmer/symplectic/base.py @@ -199,7 +199,7 @@ def _from_matrix_projector(cls, base_projector = get_PauliwordOp_projector('0'*n_qubits) P_out = cls.empty(n_qubits) for i,j in zip(row, col): - left = np.binary_repr(i, width=n_qubits) + left = np.binary_repr(i, width=n_qubits) right = np.binary_repr(j, width=n_qubits) left_op = cls.from_list([left.replace('0', 'I').replace('1', 'X')]) right_op = cls.from_list([right.replace('0', 'I').replace('1', 'X')]) @@ -1442,9 +1442,9 @@ def get_PauliwordOp_projector(projector: Union[str, List[str], np.array]) -> "Pa else: projector = np.asarray(projector) basis_dict = {'I':1, - '0':1, '1':-1, - '+':1, '-':-1, - '*':1, '%':-1} + '0':0, '1':1, + '+':0, '-':1, + '*':0, '%':1} assert len(projector.shape) == 1, 'projector can only be defined over a single string or single list of strings (each a single letter)' assert set(projector).issubset(list(basis_dict.keys())), 'unknown qubit state (must be I,X,Y,Z basis)' @@ -1462,13 +1462,15 @@ def get_PauliwordOp_projector(projector: Union[str, List[str], np.array]) -> "Pa dtype=object))[ ::-1])) > 0).astype(int) - # assign a sign only to 'active positions' (0 in binary not relevent) - sign_from_binary = binary_vec * state_sign + # # assign a sign only to 'active positions' (0 in binary not relevent) + # sign_from_binary = binary_vec * state_sign + # + # # need to turn 0s in matrix to 1s before taking product across rows + # sign_from_binary = sign_from_binary + (sign_from_binary + 1) % 2 + # + # sign = np.product(sign_from_binary, axis=1) - # need to turn 0s in matrix to 1s before taking product across rows - sign_from_binary = sign_from_binary + (sign_from_binary + 1) % 2 - - sign = np.product(sign_from_binary, axis=1) + sign = (-1)**((binary_vec@state_sign.T)%2) coeff = 1 / 2 ** (N_qubits_fixed) * np.ones(2 ** N_qubits_fixed) sym_arr = np.zeros((coeff.shape[0], 2 * N_qubits))