diff --git a/notebooks/Matrix to PauliwordOp.ipynb b/notebooks/Matrix to PauliwordOp.ipynb
index 74aa5a55..217dbd1e 100644
--- a/notebooks/Matrix to PauliwordOp.ipynb	
+++ b/notebooks/Matrix to PauliwordOp.ipynb	
@@ -410,7 +410,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.13"
+   "version": "3.9.7"
   },
   "vscode": {
    "interpreter": {
diff --git a/notebooks/hilbert schmidt inner product.ipynb b/notebooks/hilbert schmidt inner product.ipynb
new file mode 100644
index 00000000..34358ed4
--- /dev/null
+++ b/notebooks/hilbert schmidt inner product.ipynb	
@@ -0,0 +1,1176 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "1ff590e9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from symmer.symplectic import PauliwordOp\n",
+    "import numpy as np\n",
+    "from scipy.sparse import rand"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "91670474",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(4, 4)"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "n_qubits = 2\n",
+    "\n",
+    "D = 2**n_qubits\n",
+    "\n",
+    "A = rand(D, D, density=0.5, format='csr')\n",
+    "B = rand(D, D, density=0.5, format='csr')\n",
+    "A.shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "9ccf6e66",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# np.intersect1d(A.nonzero(), B.nonzero())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0a32b2ea",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "02ec37e0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "inds = np.array(list(zip(A.nonzero()[0], A.nonzero()[1])))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "6bd24219",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "matrix([[0.20998285, 0.02277784, 0.97400745, 0.7810161 , 0.1003445 ,\n",
+       "         0.52807743, 0.51582736, 0.81581773]])"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "A[A.nonzero()]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "83b71154",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[0.        , 0.20998285, 0.02277784, 0.        ],\n",
+       "       [0.        , 0.97400745, 0.7810161 , 0.        ],\n",
+       "       [0.        , 0.1003445 , 0.52807743, 0.51582736],\n",
+       "       [0.81581773, 0.        , 0.        , 0.        ]])"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "A.toarray()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "9694d227",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.3607448291650869"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "A_terms = np.array(list(zip(A.nonzero()[0], A.nonzero()[1])))\n",
+    "B_terms = np.array(list(zip(B.nonzero()[0], B.nonzero()[1])))\n",
+    "# keep = []\n",
+    "# for a in A_terms:\n",
+    "#     for b in B_terms:\n",
+    "#         if np.all(a==b):\n",
+    "#             keep.append(a)\n",
+    "# survived = np.array(keep)\n",
+    "\n",
+    "aset = set([tuple(x) for x in A_terms])\n",
+    "bset = set([tuple(x) for x in B_terms])\n",
+    "survived = np.array([x for x in aset & bset])\n",
+    "\n",
+    "# np.einsum('i,i->', A[survived[:,0]], A[survived[:,1]], B[survived[:,0]], B[survived[:,1]]\n",
+    "\n",
+    "# np.einsum('i,i->', A[survived[:,0], survived[:,1]], B[survived[:,0], survived[:,1]])\n",
+    "\n",
+    "A_survied = np.array(A[survived[:,0], survived[:,1]])[0]\n",
+    "B_survied = np.array(B[survived[:,0], survived[:,1]])[0]\n",
+    "np.einsum('i,i->', A_survied, B_survied)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "05ff8942",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.3607448291650868"
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.trace(A.toarray().T@B.toarray())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "id": "dc606d7c",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([0.02277784, 0.7810161 , 0.1003445 , 0.52807743])"
+      ]
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.array(A[survived[:,0], survived[:,1]])[0]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5aeb0341",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bfad087b",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "886a9e75",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cf659d7b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "i_terms = np.intersect1d(A.nonzero()[0], B.nonzero()[0])\n",
+    "j_terms = np.intersect1d(A.nonzero()[1], B.nonzero()[1])\n",
+    "\n",
+    "# i_terms = np.union1d(A.nonzero()[0], B.nonzero()[0])\n",
+    "# j_terms = np.union1d(A.nonzero()[1], B.nonzero()[1])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7c977b9f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# A[i_terms, j_terms]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5dc7a47d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "B.todense()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fa8b0891",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "A.todense()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fb9b0ba9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.einsum('ij,ij->',A[i_terms, j_terms], B[i_terms, j_terms])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "71894deb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.einsum('ij,ji->',A.T.toarray(), B.toarray())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "db888d1d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.trace(A.T.toarray() @ B.toarray())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "43af26f8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.einsum('ij,ij->',A.toarray(), B.toarray())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "aee23d44",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "i_terms= A.nonzero()[0]\n",
+    "j_terms = A.nonzero()[1]\n",
+    "\n",
+    "minA = A[*i_terms.tolist(), *j_terms.tolist()]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "13c13bb6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "W = A[i_terms.tolist()]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "59c99f1e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "W.shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "aaf6f2bf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "i_terms"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ae0acaa9",
+   "metadata": {},
+   "source": [
+    "$$Tr(A^{\\dagger} B) = \\sum_{ij}^{2^{n}} A_{ij} B_{ij}$$\n",
+    "\n",
+    "- only need to sum over nonzero $ij$ pairs of indices!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 131,
+   "id": "d302bdd9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def hilbert_schmidt_inner_product(mat1, mat2):\n",
+    "    \n",
+    "    ## find nonzero ij indices of each matrix\n",
+    "    mat1_ij = set([tuple(x) for x in np.array(list(zip(*mat1.nonzero())))])\n",
+    "    mat2_ij = set([tuple(x) for x in np.array(list(zip(*mat2.nonzero())))])\n",
+    "    \n",
+    "    ## find common ij indices between these that are both nonzero\n",
+    "    common_ij = np.array(list(mat1_ij & mat2_ij))\n",
+    "    \n",
+    "    ## select common i,j indices from both (now will be 1D array)\n",
+    "    mat1_survied = np.array(mat1[common_ij[:,0], common_ij[:,1]])[0]\n",
+    "    mat2_survied = np.array(mat2[common_ij[:,0], common_ij[:,1]])[0]\n",
+    "    \n",
+    "    ## multiply the nonzero ij common elements (1D dot product!)\n",
+    "    trace = np.dot(mat1_survied.conj(),mat2_survied)\n",
+    "    return trace"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 124,
+   "id": "8ca78926",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n_qubits = 12\n",
+    "\n",
+    "D = 2**n_qubits\n",
+    "\n",
+    "A = rand(D, D, density=0.001, format='csr')\n",
+    "B = rand(D, D, density=0.001, format='csr')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 125,
+   "id": "be80da65",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# n_qubits = 8\n",
+    "\n",
+    "# D = 2**n_qubits\n",
+    "\n",
+    "# A = rand(D, D, density=0.1, format='csr')\n",
+    "# B = rand(D, D, density=0.1, format='csr')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 126,
+   "id": "a9edc7e6",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "49.2 ms ± 3.13 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit hilbert_schmidt_inner_product(A,B)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 127,
+   "id": "ffe87d5f",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "53.9 ms ± 2.74 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit np.einsum('ij,ij->', A.conj().todense(), B.todense())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 128,
+   "id": "d1365a3f",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "67.7 ms ± 1.24 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit np.einsum('ij,ij->', A.conj().todense(), B.todense(), optimize=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 129,
+   "id": "25c5f5c6",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "1.48 ms ± 32 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit sum((A.conj().T@B).diagonal())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 134,
+   "id": "f45f7ddd",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "5.217258722681082"
+      ]
+     },
+     "execution_count": 134,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "sum((A.conj().T@B).diagonal())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 135,
+   "id": "91fbd650",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "5.217258722681082"
+      ]
+     },
+     "execution_count": 135,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "hilbert_schmidt_inner_product(A,B)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 136,
+   "id": "b9fc7638",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "5.217258722681082"
+      ]
+     },
+     "execution_count": 136,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.einsum('ij,ij->', A.conj().todense(), B.todense())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ec6692a7",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 218,
+   "id": "88aa7825",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n_qubits = 14\n",
+    "\n",
+    "D = 2**n_qubits\n",
+    "\n",
+    "A = rand(D, D, density=0.001, format='csr')\n",
+    "B = rand(D, D, density=0.001, format='csr')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "22a81ade",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 223,
+   "id": "68c81d95",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "A_dense = A.toarray()\n",
+    "B_dense = B.toarray()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 224,
+   "id": "0fb55f44",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "202 ms ± 15.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit np.einsum('ij, ij', A_dense, B_dense)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 226,
+   "id": "748be2ff",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "63.86005862182765"
+      ]
+     },
+     "execution_count": 226,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "(A.multiply(B)).sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 227,
+   "id": "e7f2e219",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "63.86005862182765"
+      ]
+     },
+     "execution_count": 227,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.einsum('ij, ij', A_dense, B_dense)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 217,
+   "id": "789080eb",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "1.32 s ± 42.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit np.sum(A_dense * B_dense)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 225,
+   "id": "c6795827",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "4.53 ms ± 196 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit (A.conj().multiply(B)).sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 222,
+   "id": "6924ccf3",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "51.3 ms ± 745 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit sum((A.conj().T @ B).diagonal())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 210,
+   "id": "d37fc043",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "1.42 s ± 205 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
+     ]
+    }
+   ],
+   "source": [
+    "%timeit (A_dense.conj() * (B_dense)).sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 212,
+   "id": "f391ce4c",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<1x16384 sparse matrix of type '<class 'numpy.float64'>'\n",
+       "\twith 16384 stored elements in Compressed Sparse Row format>"
+      ]
+     },
+     "execution_count": 212,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 196,
+   "id": "cb5c6ac0",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "10559.513739928236"
+      ]
+     },
+     "execution_count": 196,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 170,
+   "id": "563d17dc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %timeit sum((A.conj().T @ B).diagonal())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 150,
+   "id": "782887a6",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "matrix([[0.52829966, 0.31149272, 0.        ],\n",
+       "        [0.89450939, 0.59275915, 0.18832217],\n",
+       "        [0.69545409, 0.45883224, 0.11505853]])"
+      ]
+     },
+     "execution_count": 150,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# B = rand(3, 3, density=0.9, format='csr')\n",
+    "# B.todense()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 158,
+   "id": "82fcbeac",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sum((A.conj().T @ B).diagonal())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 159,
+   "id": "05e40e0b",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "1048435.4757661638"
+      ]
+     },
+     "execution_count": 159,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "(A.conj().multiply(B)).sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 166,
+   "id": "b62f7af3",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "1048435.4757661638"
+      ]
+     },
+     "execution_count": 166,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "AB = A.conj().T @ B\n",
+    "np.trace(AB.toarray())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 165,
+   "id": "d235721a",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "1048435.4757661638"
+      ]
+     },
+     "execution_count": 165,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.trace(AB.toarray())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "38ea09fc",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 94,
+   "id": "0086742d",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "matrix([[0.],\n",
+       "        [0.],\n",
+       "        [0.],\n",
+       "        ...,\n",
+       "        [0.],\n",
+       "        [0.],\n",
+       "        [0.]])"
+      ]
+     },
+     "execution_count": 94,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "d = (A.multiply(B.T)).sum(axis=1)\n",
+    "d"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e674071a",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "32d6b5ef",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 47,
+   "id": "7c15ca65",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "matrix([[0.31685082, 0.39924373, 0.        , 0.93879449, 0.42521169,\n",
+       "         0.23145935, 0.17405214, 0.69974732, 0.        , 0.00576625,\n",
+       "         0.        , 0.        ]])"
+      ]
+     },
+     "execution_count": 47,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 42,
+   "id": "37500ed3",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "matrix([[0.02587247, 0.84722346, 0.85563884, ..., 0.55659597, 0.35041813,\n",
+       "         0.52862651]])"
+      ]
+     },
+     "execution_count": 42,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "A[A.nonzero()]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 45,
+   "id": "c09dbb53",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(4, 4)"
+      ]
+     },
+     "execution_count": 45,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "n_qubits = 2\n",
+    "\n",
+    "D = 2**n_qubits\n",
+    "\n",
+    "A = rand(D, D, density=0.5, format='csr')\n",
+    "B = rand(D, D, density=0.5, format='csr')\n",
+    "A.shape"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 59,
+   "id": "c24ef62e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "non_Z = (abs(A) + abs(B)).nonzero()\n",
+    "survived = np.array(list(zip(*non_Z)))\n",
+    "\n",
+    "A_survied = np.array(A[survived[:,0], survived[:,1]])[0]\n",
+    "B_survied = np.array(B[survived[:,0], survived[:,1]])[0]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 64,
+   "id": "439e8df1",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "matrix([[1.28363658, 0.39924373, 0.42043598, 0.93879449],\n",
+       "        [0.50186747, 0.27474963, 0.        , 0.88260442],\n",
+       "        [0.69974732, 0.        , 0.12702344, 0.00576625],\n",
+       "        [0.        , 0.22478913, 0.2690172 , 0.        ]])"
+      ]
+     },
+     "execution_count": 64,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "(A+B).todense()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 69,
+   "id": "97ae477c",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "matrix([[0.31685082, 0.39924373, 0.93879449, 0.42521169, 0.23145935,\n",
+       "         0.17405214, 0.69974732, 0.00576625]])"
+      ]
+     },
+     "execution_count": 69,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "A[A.nonzero()]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 49,
+   "id": "4c640e29",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(array([0, 0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3], dtype=int32),\n",
+       " array([0, 1, 2, 3, 0, 1, 3, 0, 2, 3, 1, 2], dtype=int32))"
+      ]
+     },
+     "execution_count": 49,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "non_Z"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 51,
+   "id": "6238e009",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "matrix([[0.31685082, 0.39924373, 0.        , 0.93879449],\n",
+       "        [0.42521169, 0.23145935, 0.        , 0.17405214],\n",
+       "        [0.69974732, 0.        , 0.        , 0.00576625],\n",
+       "        [0.        , 0.        , 0.        , 0.        ]])"
+      ]
+     },
+     "execution_count": 51,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "A.todense()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 52,
+   "id": "03b2ccac",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "matrix([[0.96678575, 0.        , 0.42043598, 0.        ],\n",
+       "        [0.07665578, 0.04329028, 0.        , 0.70855227],\n",
+       "        [0.        , 0.        , 0.12702344, 0.        ],\n",
+       "        [0.        , 0.22478913, 0.2690172 , 0.        ]])"
+      ]
+     },
+     "execution_count": 52,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "B.todense()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 54,
+   "id": "dff14f1a",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "matrix([[0.31685082, 0.39924373, 0.        , 0.93879449],\n",
+       "        [0.42521169, 0.23145935, 0.        , 0.17405214],\n",
+       "        [0.69974732, 0.        , 0.        , 0.00576625],\n",
+       "        [0.        , 0.        , 0.        , 0.        ]])"
+      ]
+     },
+     "execution_count": 54,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "A.todense()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dab4d104",
+   "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
+}