asserts the completeness condition: b0@b0.T.conj() + b1@b1.T.conj() = I\n", + " *is the Moore-Penrose pseudo-inverse of \n", + " \"\"\"\n", + " mid = int(start + (end-start)/2)\n", + " # computing \n", + " M,D,V = get_diagonalization(povm, start, end)\n", + " # computing the null space of \n", + " P = np.sign(D.round(5))\n", + " Pc = np.eye(len(M))-P\n", + " Q = V@Pc@V.T.conj()\n", + " # computing \n", + " #D_inv = np.zeros_like(D)\n", + " #for i in range(len(M)):\n", + " # if D[i,i].round(5) == 0: continue\n", + " # else: D_inv[i,i] = np.real(1/D[i,i])\n", + " D_inv = np.linalg.pinv(D)\n", + " M_psinv = V@D_inv@V.T.conj()\n", + " # computing and \n", + " M0,_,_ = get_diagonalization(povm, start, mid)\n", + " M1,_,_ = get_diagonalization(povm, mid, end)\n", + " # computing and \n", + " b0 = M0@M_psinv + Q/np.sqrt(2)\n", + " b1 = M1@M_psinv + Q/np.sqrt(2)\n", + " return b0, b1\n", + "\n", + "def closest_unitary(A):\n", + " \"\"\" Calculates the unitary matrix U that is closest with respect to the\n", + " operator norm distance to the general matrix A.\n", + " \"\"\"\n", + " V, _, Wh = scipy.linalg.svd(A)\n", + " U = np.matrix(V.dot(Wh))\n", + " return U\n", + "\n", + "def closest_unitary_schur(A):\n", + " T, Z = scipy.linalg.schur(A, output='complex')\n", + " return Z @ np.diag(np.diag(T)/abs(np.diag(T))) @ Z.T.conj()\n", + "\n", + "def extend_to_unitary(b0, b1):\n", + " \"\"\" Creates a coupling unitary between the system and ancilla.\n", + "\n", + " The condition for unitary is: <0|U|0> = b0 and <1|U|0> = b1,\n", + " whereby the ancilla is projected onto states |0> and |1>.\n", + "\n", + " A two-column matrix A, with its upper left block given by b0, \n", + " and bottom left block by b1, is extended to unitary U by appending \n", + " a basis of the null space of A.\n", + " \"\"\" \n", + " A = np.concatenate((b0,b1))\n", + " u, _, _ = scipy.linalg.svd(A)\n", + " y = u[:, len(A[0]):]\n", + " U = np.hstack((A, y))\n", + " # verify U is close to unitary\n", + " assert np.allclose(U.T.conj()@U, np.eye(len(A)),atol=1e-03), \"Failed to construct U\"\n", + " return closest_unitary(U)\n", + "\n", + "def theoretical_probs(povm, state, binary_tree=False):\n", + " \"\"\" Precalculates the expected probabilities for a pure state\n", + " for the given POVM.\n", + " \"\"\" \n", + " rho = state.T.conj()@state # density matrix of the state\n", + " probs = [np.trace(povm[i]@rho).round(2).real for i in range(len(povm))] \n", + " \n", + " l = len(probs)\n", + " n = int(np.ceil(np.log2(l)))\n", + " prob_counts = {}\n", + " for i in range(l):\n", + " if binary_tree: key = np.binary_repr(i, n)[::-1]\n", + " else: key = np.binary_repr(i, n)\n", + " value = probs[i]\n", + " prob_counts[key] = value\n", + " \n", + " return prob_counts\n", + "\n", + "# Class definitions\n", + "class POVM:\n", + " \"\"\"Base class that holds an arbitrary POVM as a list of POVM elements.\n", + " \"\"\"\n", + " def __init__(self, povm):\n", + " \"\"\"\n", + " Constructor asserts that the given POVM is valid.\n", + " \"\"\"\n", + " self.povm = povm\n", + " self.N = len(povm)\n", + " self.depth = int(np.ceil(np.log2(self.N))) # required depth of the binary tree\n", + " self.povm_dim = len(povm[0]) # dimension of the POVM operators\n", + " self.n_qubits = int(np.log2(self.povm_dim)) # number of system qubits\n", + " assert self.is_valid()\n", + " def is_valid(self):\n", + " \"\"\"Verifies the hermiticity, positivity of the POVM and that\n", + " the POVM resolves the identity.\n", + " Returns:\n", + " True: if all conditions are satisfied\n", + " Raises:\n", + " Assertion Error: if one of conditions is not satisfied\n", + " \"\"\"\n", + " for E in self.povm:\n", + " assert np.allclose(E.conj().T, E), \"Some POVM elements are not hermitian\"\n", + " assert np.all(np.linalg.eigvals(E).round(3) >= 0), \"Some POVM elements are not positive semi-definite\"\n", + " assert np.allclose(sum(self.povm), np.eye(self.povm_dim)), \"POVM does not resolve the identity\"\n", + " return True\n", + "\n", + "class BinaryTreePOVM(POVM):\n", + " \"\"\"Class which implements the binary tree approach as described in https://arxiv.org/abs/0712.2665 \n", + " to contruct a POVM measurement tree. \n", + " \"\"\"\n", + " def __init__(self, povm, dual=False):\n", + " \"\"\"Creates a binary stree structure with BinaryMeasurementNode objects\n", + " stored in dictionary. The keys in the dictionary are the bitstrings \n", + " corresponding to the states of the classical register at the point when the\n", + " corresponding node has been \"reached\".\n", + "\n", + " Args:\n", + " povm: list of POVM elements\n", + " \"\"\"\n", + " super().__init__(povm)\n", + " # pad with zero operators if necessary\n", + " while np.log2(self.N)-self.depth != 0:\n", + " self.povm.append(np.zeros_like(self.povm[0]))\n", + " self.N += 1\n", + " \n", + " self.nodes = {}\n", + " self.dual = dual\n", + "\n", + " self.create_binary_tree(key=\"0\", start=0, end=self.N)\n", + " self.qc = self.construct_measurement_circuit()\n", + " def create_binary_tree(self, key, start, end):\n", + " \"\"\"Recursive method to build the measurement tree.\n", + " Terminates when the fine-grain level corresponding to the single POVM\n", + " elements is reached.\n", + "\n", + " and are the first and (last-1) indices of POVM elements\n", + " which were grouped together to obtain a cumulative coarse-grain operator. \n", + " The range [start, end) corresponds to the possible outcomes which \"sit\" in \n", + " the branches below.\n", + " \"\"\"\n", + " if start >= (end-1):\n", + " return\n", + " new_node = BinaryMeasurementNode(self.povm, key=key, start=start, end=end)\n", + " self.nodes[key] = new_node\n", + " mid = int(start + (end-start)/2)\n", + " self.create_binary_tree(new_node.left, start=start, end=mid)\n", + " self.create_binary_tree(new_node.right, start=mid, end=end) \n", + " def construct_measurement_circuit(self):\n", + " \"\"\"Contructs a quantum circuit for a given POVM by sequentially appending\n", + " coupling unitaries and measurements conditioned on the state of the\n", + " classical register . The method uses BFS traversal of the precomputed \n", + " binary measurement tree, i.e. the measurement nodes are visited in level-order.\n", + " \n", + " * Traversal terminates when the fine-grain level was reached.\n", + " * Ancilla qubit is reset before each level.\n", + " * The root node has the key \"0\".\n", + "\n", + " * The instruction is applied to the entire classical register ,\n", + " whereby the value is the key of the corresponding node - padded with zeros \n", + " from right to the length of the register - and interpreted as an integer.\n", + " \n", + " Example:\n", + " At the first level the two nodes have keys:\n", + " left = \"00\" and right = \"01\"\n", + " If the is 3 bits long, then the left/right unitary is applied if \n", + " the state of is int(\"000\",2) = 0 / int(\"010\",2) = 2\n", + " \"\"\"\n", + " system = QuantumRegister(self.n_qubits, name='system')\n", + " anc = QuantumRegister(1,name='ancilla')\n", + " cr = ClassicalRegister(self.depth)\n", + " qc = QuantumCircuit(system, anc, cr, name=\"measurement-circuit\")\n", + "\n", + " root = self.nodes[\"0\"]\n", + " U_gate = UnitaryGate(root.U, label=root.key)\n", + " qc.append(U_gate, system[:] + anc[:])\n", + " if self.dual: # works for depth = 2\n", + " qc.x(anc)\n", + " qc.measure(anc, cr[0])\n", + " if self.depth == 1: return qc\n", + " qc.x(anc).c_if(cr[0],1)\n", + " \n", + " current_level = [self.nodes[\"00\"],self.nodes[\"10\"]]\n", + "\n", + " for i in range(1,self.depth):\n", + " next_level = []\n", + " for node in current_level:\n", + " U_gate = UnitaryGate(node.U, label=node.key)\n", + " cr_state = int(node.key[:-1],2)\n", + " with qc.if_test((cr, cr_state)):\n", + " qc.append(U_gate, system[:]+anc[:])\n", + " if node.left in self.nodes: next_level.append(self.nodes[node.left])\n", + " if node.right in self.nodes: next_level.append(self.nodes[node.right])\n", + " current_level = next_level\n", + " # dual condition must be checked here for larger systems\n", + " qc.measure(anc, cr[i])\n", + " if i == self.depth-1: continue\n", + " qc.x(anc).c_if(cr[i],1) # instead of resetting apply conditional X gate\n", + " \n", + " return qc\n", + "\n", + "class BinaryMeasurementNode(POVM):\n", + " \"\"\"A BinaryMeasurementNode object is a node in the BinaryTreePOVM.\n", + " It contains:\n", + " 1. Its in the dictionary.\n", + " 2. and : the first and (last-1) indices of the accumulated \n", + " POVM elements, corresponding this node.\n", + " 3. Coupling unitary .\n", + " 4. Keys and of the two children nodes.\n", + " 5. Attributes of the POVM class: , , .\n", + " 6. Its level in the binary tree, where level of the root node is 0.\n", + " \"\"\"\n", + " def __init__(self, povm, key, start, end):\n", + " super().__init__(povm)\n", + " self.key = key\n", + " self.level = len(self.key)-1\n", + " self.start = start\n", + " self.end = end\n", + " self.left = \"0\" + self.key\n", + " self.right = \"1\" + self.key\n", + " b0, b1 = get_next_level_binary_kraus_ops(self.povm, self.start, self.end)\n", + " self.U = extend_to_unitary(b0, b1)\n", + " def __str__(self):\n", + " line1 = 'Node with the key {} at level {}\\n'.format(self.key, self.level)\n", + " line2 = 'Cumulative operator = [{},{})'.format(self.start, self.end)\n", + " line3 = 'left = {}, right = {}\\n'.format(self.left, self.right)\n", + " line4 = 'U = \\n{}\\n'.format(self.U)\n", + " return line1+line2+line3+line4" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Implementing Naimark's extension circuit" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [], + "source": [ + "from numpy.linalg import svd\n", + "from qiskit.visualization import plot_histogram\n", + "\n", + "def check_for_rank_one(povm):\n", + " \"\"\"\n", + " function to check if a povm is a rank-1 povm\n", + " \"\"\"\n", + " rank_one = True\n", + " for p in povm:\n", + " if np.linalg.matrix_rank(p)!=1:\n", + " rank_one = False\n", + " return rank_one\n", + " else:\n", + " continue\n", + " return rank_one\n", + " \n", + "# %%\n", + "def compute_rank_one_unitary(povm, atol=1e-13, rtol=0):\n", + " \"\"\"\n", + " This function computes the unitary that rotates the system to the Hilbert space of the ancilla\n", + " Input: POVM ---> a list of the elements of POVM\n", + " Output: Unitary matrix\n", + " \"\"\"\n", + " \n", + " # check if povm is a rank-1 povm:\n", + " assert check_for_rank_one(povm), \"This is not a rank-1 povm\"\n", + " new_povm = []\n", + " for p in povm:\n", + " if np.log2(len(povm))%2==0: #still under investigation\n", + " w, v = np.linalg.eig(p)\n", + " else: \n", + " w, v = np.linalg.eigh(p) #note the that the eigenvenvector is computer for hermitian eigh\n", + " for eigenvalue, engenvector in zip(w,v):\n", + " if np.isclose(np.abs(eigenvalue), 0):\n", + " continue\n", + " else:\n", + " new_p = np.sqrt(eigenvalue)*engenvector\n", + " new_povm.append(new_p)\n", + " v = np.vstack(new_povm) # arrange the povm elements to form a matrix of dimension: Row X Col*len(povm)\n", + " v = np.atleast_2d(v.T) # convert to 2d matrix\n", + " \n", + " u, s, vh = svd(v) # apply svd\n", + " tol = max(atol, rtol * s[0])\n", + " nnz = (s >= tol).sum()\n", + " ns = vh[nnz:] # missing rows of v\n", + " \n", + " # add the missing rows of v to v\n", + " V = np.vstack((v, ns)) \n", + " \n", + " \n", + " # make the unitary a square matrix of dimension N=2^n where n = int(np.ceil(np.log(V.shape[0])))\n", + " n = int(np.ceil(np.log2(V.shape[0])))\n", + " N = 2**n # dimension of system and ancilla Hilber space\n", + " r,c = V.shape \n", + " \n", + " U = np.eye(N, dtype=complex) # initialize Unitary matrix to the identity. Ensure it is complex\n", + " U[:r,:c] = V[:r,:c] # assign all the elements of V to the corresponding elements of U\n", + " \n", + " U = U.conj().T # Transpose the unitary so that the rows are the povm\n", + " \n", + " # check for unitarity of U\n", + " assert np.allclose(U.T.conj()@U, np.eye(N),atol=1e-13), \"Failed to construct U\"\n", + " \n", + " return U\n", + " \n", + "# %%\n", + "# Using the original unitary generator\n", + "\n", + "def compute_full_rank_unitary(povm, atol=1e-13, rtol=0):\n", + " \"\"\"\n", + " This function computes the unitary that rotates the system to the Hilbert space of the ancilla\n", + " Input: POVM ---> a list of the elements of POVM\n", + " Output: Unitary matrix\n", + " \"\"\"\n", + " \n", + " \n", + " # Here square root of the POVM elements were used as a replacement for the vector that form the povm\n", + " povm = [sqrtm(M)for M in povm]\n", + " \n", + " v = np.hstack(povm) # arrange the povm elements to form a matrix of dimension: Row X Col*len(povm)\n", + " v = np.atleast_2d(v) # convert to 2d matrix\n", + " u, s, vh = svd(v) # apply svd\n", + " tol = max(atol, rtol * s[0])\n", + " nnz = (s >= tol).sum()\n", + " ns = vh[nnz:] # missing rows of v\n", + " \n", + " # add the missing rows of v to v\n", + " V = np.vstack((v, ns)) \n", + " \n", + " \n", + " # make the unitary a square matrix of dimension N=2^n where n = int(np.ceil(np.log(V.shape[0])))\n", + " n = int(np.ceil(np.log2(V.shape[0])))\n", + " N = 2**n # dimension of system and ancilla Hilber space\n", + " r,c = V.shape \n", + " \n", + " U = np.eye(N, dtype=complex) # initialize Unitary matrix to the identity. Ensure it is complex\n", + " U[:r,:c] = V[:r,:c] # assign all the elements of V to the corresponding elements of U\n", + " \n", + " U = U.conj().T # Transpose the unitary so that the rows are the povm\n", + " \n", + " # check for unitarity of U\n", + " assert np.allclose(U.T.conj()@U, np.eye(N),atol=1e-07), \"Failed to construct U\"\n", + " \n", + " return U\n", + "\n", + "# %%\n", + "def rank_one_circuit(povm, state, U):\n", + " \n", + " # Define the quantum and classical registers\n", + " dim_system = state.shape[1] # dimension of state\n", + " num_system_qubit = int(np.ceil(np.log2(dim_system))) # total number of qubits for system\n", + "\n", + " system_reg = QuantumRegister(num_system_qubit, name='system') # system register\n", + "\n", + " N = U.shape[0] # Dimension of the unitary to be applied to system and ancilla\n", + "\n", + " num_ancilla_qubit = int(np.ceil(np.log2(N))) - num_system_qubit # total number of qubits for system\n", + "\n", + " ancilla_reg = QuantumRegister(num_ancilla_qubit, name='ancilla') # ancilla register\n", + "\n", + " U_gate = UnitaryGate(U, label='U') # unitary gate to be applied between system and ancilla\n", + " \n", + " \n", + " # create the quantum circuit for the system and ancilla\n", + " qc = QuantumCircuit(system_reg, ancilla_reg, name='circuit')\n", + " qc.initialize(state[0],system_reg)\n", + "\n", + " # reset ancilla to zero\n", + " qc.reset(ancilla_reg)\n", + "\n", + " # append the unitary gate\n", + " qc.append(U_gate, range(system_reg.size + ancilla_reg.size))\n", + "\n", + " # measure only the ancilliary qubits\n", + " qc.measure_all()\n", + " \n", + " return qc\n", + "\n", + "# %%\n", + "def full_rank_circuit(povm, state, U):\n", + " \n", + " # Define the quantum and classical registers\n", + " dim_system = state.shape[1] # dimension of state\n", + " num_system_qubit = int(np.ceil(np.log2(dim_system))) # total number of qubits for system\n", + "\n", + " system_reg = QuantumRegister(num_system_qubit, name='system') # system register\n", + "\n", + " N = U.shape[0] # Dimension of the unitary to be applied to system and ancilla\n", + "\n", + " num_ancilla_qubit = int(np.ceil(np.log2(N))) - num_system_qubit # total number of qubits for system\n", + "\n", + " ancilla_reg = QuantumRegister(num_ancilla_qubit, name='ancilla') # ancilla register\n", + "\n", + " classical_reg = ClassicalRegister(num_ancilla_qubit, name='measure') # classical register\n", + "\n", + " U_gate = UnitaryGate(U, label='U') # unitary gate to be applied between system and ancilla\n", + " \n", + " \n", + " # create the quantum circuit for the system and ancilla\n", + " qc = QuantumCircuit(system_reg, ancilla_reg, classical_reg, name='circuit')\n", + " qc.initialize(state[0],system_reg)\n", + "\n", + " # reset ancilla to zero\n", + " qc.reset(ancilla_reg)\n", + "\n", + " # append the unitary gate\n", + " qc.append(U_gate, range(system_reg.size + ancilla_reg.size))\n", + "\n", + " # measure only the ancilliary qubits\n", + " qc.measure(ancilla_reg, classical_reg)\n", + " \n", + " return qc\n", + "\n", + "# %%\n", + "def construct_quantum_circuit(povm, state):\n", + " \n", + " # compute unitary matrix\n", + " if check_for_rank_one(povm):\n", + " U = compute_rank_one_unitary(povm)\n", + " qc = rank_one_circuit(povm, state, U)\n", + " else:\n", + " U = compute_full_rank_unitary(povm)\n", + " qc = full_rank_circuit(povm, state, U)\n", + " \n", + " return qc\n", + "\n", + "# %%\n", + "def draw_circuit(qc, idle_wires=True):\n", + " \"\"\"\n", + " This functions draws the naimark extension quantum circuit\n", + " \"\"\"\n", + " \n", + " return qc.draw(output='mpl', idle_wires=idle_wires)\n", + "\n", + "def naimark_plot(povm, state, counts, names, save=False, file_name=\"povm_output.pdf\"): \n", + " theory_count = theoretical_probs(povm, state)\n", + " count_list = [theory_count]\n", + " #legend_list = [\"Theoretical result\"]\n", + " legend_list = ['Theory (t)']\n", + " \n", + " for count, name in zip(counts, names):\n", + " count_list.append(count)\n", + " legend_list.append(name)\n", + " \n", + " fig = plot_histogram(count_list, legend=legend_list, bar_labels=True)\n", + " ax = fig.axes[0]\n", + " title = f\"Tetrad POVM with Naimark's extension approach\"\n", + " ax.set_title(title, fontsize=12)\n", + " ax.set_ylabel('Probabilities', fontsize=12)\n", + " ax.set_xlabel('POVM elements', fontsize=12)\n", + " \n", + " x_lables = []\n", + " for i in range(len(povm)):\n", + " x_lables.append(\"M\"+str(i))\n", + " ax.set_xticks(range(len(povm)))\n", + " ax.set_xticklabels(x_lables)\n", + " \n", + " ax.legend(fontsize=12)\n", + " fig.tight_layout()\n", + " \n", + " if save:\n", + " fig.savefig(file_name)\n", + " \n", + " \n", + " return fig" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Qiskit version and copyright " + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + " Version Information
" + ], + "text/plain": [ + "
Qiskit Software Version qiskit-terra0.22.3 qiskit-aer0.11.2 qiskit-ibmq-provider0.19.2 qiskit0.39.4 System information Python version 3.8.8 Python compiler MSC v.1916 64 bit (AMD64) Python build default, Apr 13 2021 15:08:03 OS Windows CPUs 4 Memory (Gb) 7.857898712158203 Wed Jan 11 12:10:31 2023 W. Europe Standard Time " + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + " " + ], + "text/plain": [ + "This code is a part of Qiskit
© Copyright IBM 2017, 2023.
This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals." + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import qiskit.tools.jupyter\n", + "%qiskit_version_table\n", + "%qiskit_copyright" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3.8.8 ('qiskit_env': venv)", + "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.8.8" + }, + "vscode": { + "interpreter": { + "hash": "d40343b2b8fe07cae983373af2e631c3feb542de4b1fdf9434c5290f2e9f60fc" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}