diff --git a/docs/tomography.rst b/docs/tomography.rst index 4f46e50a..999c5af3 100644 --- a/docs/tomography.rst +++ b/docs/tomography.rst @@ -77,6 +77,14 @@ Finally, we analyze our data with one of the analysis routines:: [-0.1869-0.0077j -0.1794-0.0188j -0.169 -0.0169j 0.2202-0.j ]] Purity = (0.6889520199999999+4.597017211338539e-17j) +Debugger +~~~~~~~~~ + +The above steps can be automated to create a basic debugger that can be used to +peek into the state of a program running on a qc. This can be done using the +tomographize function:: + + rho = tomographize(qc, program, qubits, pauli_num=10, t_type="compressed_sensing") API Reference ------------- @@ -96,4 +104,5 @@ API Reference iterative_mle_state project_density_matrix estimate_variance + tomographize diff --git a/examples/tomography_debugger.ipynb b/examples/tomography_debugger.ipynb new file mode 100644 index 00000000..19a6320a --- /dev/null +++ b/examples/tomography_debugger.ipynb @@ -0,0 +1,235 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Tomography Debugger\n", + "State tomography involves measuring a quantum state repeatedly in the bases given by `itertools.product(['X', 'Y', 'Z], repeat=n_qubits)`. From these measurements, we can reconstruct a density matrix $\\rho$ using a varaiety of methods described in forest.benchmarking.tomography under the heading \"state tomography\". This is all done automaticly in using the forest.benchmarking.tomography.tomographize function allowing it to be use effectivly as a quantum debugger." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import time\n", + "\n", + "from pyquil import Program, get_qc\n", + "from pyquil.gates import *\n", + "from forest.benchmarking.tomography import tomographize" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Construct a state with a `Program`\n", + "We'll construct a two-qubit graph state by Hadamarding all qubits and then applying a controlled-Z operation across edges of our graph. In the two-qubit case, there's only one edge. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "qubits = [0, 1]\n", + "\n", + "program = Program()\n", + "for qubit in qubits:\n", + " program += H(qubit)\n", + "program += CZ(qubits[0], qubits[1])\n", + "program += RY(-np.pi/2, qubits[0])\n", + "program += X(qubits[1])\n", + "program += CNOT(qubits[1], qubits[0])\n", + "\n", + "print(program)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run the tomography debugger and print output" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "qc = get_qc('%dq-qvm' % len(qubits))\n", + "\n", + "\n", + "start_linear = time.time()\n", + "m = 10\n", + "rho_linear = tomographize(qc, program, qubits, pauli_num=10, t_type=\"linear_inv\")\n", + "end_linear = time.time() - start_linear\n", + "\n", + "print(\"Linear tomography took %gs\" % end_linear)\n", + "print(\"Recovered density matrix:\\n\")\n", + "print(rho_linear)\n", + "\n", + "start_compressed = time.time()\n", + "rho_compressed = tomographize(qc, program, qubits, pauli_num=10, t_type=\"compressed_sensing\")\n", + "end_compressed = time.time() - start_compressed\n", + "print(\"Compressed tomography took %gs\" % end_compressed)\n", + "print(\"Recovered density matrix:\\n\")\n", + "print(rho_compressed)\n", + "\n", + "start_lasso = time.time()\n", + "rho_lasso = tomographize(qc, program, qubits, pauli_num=10, t_type=\"lasso\")\n", + "end_lasso = time.time() - start_lasso\n", + "print(\"Compressed tomography took %gs\" % end_lasso)\n", + "print(\"Recovered density matrix:\\n\")\n", + "print(rho_lasso)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Compare results to true output obtained using wavefunction simulator" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pyquil.api import WavefunctionSimulator\n", + "wf_sim = WavefunctionSimulator()\n", + "wf = wf_sim.wavefunction(program)\n", + "psi = wf.amplitudes\n", + "\n", + "rho_true = np.outer(psi, psi.T.conj())\n", + "print(np.around(rho_true, decimals=3))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualize using Hinton plots" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt\n", + "from forest.benchmarking.plotting import hinton\n", + "fig, (ax1, ax2, ax3) = plt.subplots(1, 3)\n", + "hinton(rho_true, ax=ax1)\n", + "hinton(rho_linear, ax=ax2)\n", + "hinton(rho_compressed, ax=ax3)\n", + "ax1.set_title('Analytical Linear')\n", + "ax2.set_title('Estimated Linear')\n", + "ax3.set_title('Estimated Compressed')\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Calculate matrix norm between true and estimated rho" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(\"Linear norm:\")\n", + "print(np.linalg.norm(rho_linear - rho_true))\n", + "print(\"Compressed norm:\")\n", + "print(np.linalg.norm(rho_compressed - rho_true))\n", + "print(\"Lasso norm:\")\n", + "print(np.linalg.norm(rho_lasso - rho_true))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot graph of results for various measurement values" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "max_pauli_num = 4 ** len(qubits)\n", + "num_trials = 5\n", + "\n", + "linear_norms = []\n", + "compressed_norms = []\n", + "\n", + "print(\"Analyzing performance of linear vs. compressed on program:\")\n", + "print(program)\n", + "\n", + "for i in range(1, max_pauli_num):\n", + " print(\"Running iteration %d/%d\" % (i, max_pauli_num))\n", + " linear_norm_mean = 0.0\n", + " compressed_norm_mean = 0.0\n", + " for j in range(num_trials):\n", + " rho_linear = tomographize(qc, program, qubits, pauli_num=i, t_type=\"linear_inv\")\n", + " rho_compressed = tomographize(qc, program, qubits, pauli_num=i, t_type=\"compressed_sensing\")\n", + " linear_norm_mean += np.linalg.norm(rho_linear - rho_true)\n", + " compressed_norm_mean += np.linalg.norm(rho_compressed - rho_true)\n", + " \n", + " linear_norm_mean /= num_trials\n", + " compressed_norm_mean /= num_trials\n", + " \n", + " linear_norms.append(linear_norm_mean)\n", + " compressed_norms.append(compressed_norm_mean)\n", + "\n", + "plt.plot(linear_norms, label='linear')\n", + "plt.plot(compressed_norms, label='compressed')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "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.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/forest/benchmarking/tomography.py b/forest/benchmarking/tomography.py index dd1506a6..0d997637 100644 --- a/forest/benchmarking/tomography.py +++ b/forest/benchmarking/tomography.py @@ -3,17 +3,21 @@ from operator import mul from typing import Callable, Tuple, List, Sequence import warnings +import random import numpy as np from scipy.linalg import logm, pinv, eigh +import cvxpy as cp import forest.benchmarking.distance_measures as dm from forest.benchmarking.superoperator_tools import vec, unvec, proj_choi_to_physical from forest.benchmarking.utils import n_qubit_pauli_basis +from forest.benchmarking.compilation import basic_compile from pyquil import Program from pyquil.operator_estimation import ExperimentSetting, \ TomographyExperiment as PyQuilTomographyExperiment, ExperimentResult, SIC0, SIC1, SIC2, SIC3, \ plusX, minusX, plusY, minusY, plusZ, minusZ, TensorProductState, zeros_state +from pyquil.operator_estimation import group_experiments, measure_observables from pyquil.paulis import sI, sX, sY, sZ, PauliSum, PauliTerm, is_identity from pyquil.unitary_tools import lifted_pauli as pauli2matrix, lifted_state_operator as state2matrix @@ -163,6 +167,98 @@ def linear_inv_state_estimate(results: List[ExperimentResult], rho = pinv(measurement_matrix) @ expectations return unvec(rho) +def compressed_sensing_state_estimate(results: List[ExperimentResult], + qubits: List[int]) -> np.ndarray: + """ + Estimate a quantum state using compressed sensing + + [FLAMMIA] Quantum Tomography via Compressed Sensing: Error Bounds, Sample Complexity, and Efficient Estimators + Steven T. Flammia et. al. + (2012). + https://arxiv.org/pdf/1205.2300.pdf + + :param results: A tomographically complete list of results. + :param qubits: All qubits that were tomographized. This specifies the order in + which qubits will be kron'ed together. + :return: A point estimate of the quantum state rho. + """ + pauli_num = len(results) + qubit_num = len(qubits) + d = 2 ** qubit_num + pauli_list = [] + expectation_list = [] + y = np.zeros((m,1)) + + for i in range(pauli_num): + #Convert the Pauli term into a tensor + r = results[i] + p_tensor = pauli2matrix(r.setting.out_operator, qubits) + e = r.expectation + #A[i] = e * scale_factor + pauli_list.append(p_tensor) + expectation_list.append(e) + + s = cp.Variable((d,d),complex = True) + obj = cp.Minimize(cp.norm(s, 'nuc')) + constraints = [cp.trace(s) == 1] + + for i in range(len(pauli_list)): + trace_bool = cp.trace(cp.matmul(pauli_list[i], s)) - expectation_list[i] == 0 + constraints.append(trace_bool) + + # Form and solve problem. + prob = cp.Problem(obj, constraints) + prob.solve() + rho = s.value + return rho + +def lasso_state_estimate(results: List[ExperimentResult], + qubits: List[int]) -> np.ndarray: + """ + Estimate a quantum state using compressed sensing + + [FLAMMIA] Quantum Tomography via Compressed Sensing: Error Bounds, Sample Complexity, and Efficient Estimators + Steven T. Flammia et. al. + (2012). + https://arxiv.org/pdf/1205.2300.pdf + + :param results: A tomographically complete list of results. + :param qubits: All qubits that were tomographized. This specifies the order in + which qubits will be kron'ed together. + :return: A point estimate of the quantum state rho. + """ + pauli_num = len(results) + qubit_num = len(qubits) + d = 2 ** qubit_num + pauli_list = [] + expectation_list = [] + y = np.zeros((m,1)) + + mu = 4 * pauli_num / np.sqrt(1000 * pauli_num) + + for i in range(m): + #Convert the Pauli term into a tensor + r = results[i] + p_tensor = pauli2matrix(r.setting.out_operator, qubits) + e = r.expectation + y[i] = e + pauli_list.append(p_tensor) + expectation_list.append(e) + + x = cp.Variable((d,d),complex = True) + A = cp.vstack([cp.trace(cp.matmul(pauli_list[i], x)) * np.sqrt(d / m) for i in range(m)]) + + #Minimize trace norm + obj = cp.Minimize(0.5 * cp.norm((A - y), 2) + mu * cp.norm(x, 'nuc')) + constraints = [cp.trace(x) == 1] + + # Form and solve problem. + prob = cp.Problem(obj, constraints) + prob.solve() + rho = x.value + + return rho + def iterative_mle_state_estimate(results: List[ExperimentResult], qubits: List[int], epsilon=.1, entropy_penalty=0.0, beta=0.0, tol=1e-9, maxiter=10_000) \ @@ -611,3 +707,50 @@ def _grad_cost(A, n, estimate, eps=1e-6): p = np.clip(p, a_min=eps, a_max=None) eta = n / p return unvec(-A.conj().T @ eta) + +def tomographize(qc, program: Program, qubits: List[int], num_shots=1000, t_type='lasso', pauli_num=None): + """ + Runs tomography on the state generated by program and estimates the state. Can be used as a debugger + + :param qc: the quantum computer to run the debugger on + :param program: which program to run the tomography on + :param quibits: whihc qubits to run the tomography debugger on + :param num_shots: the number of times to run each tomography experiment to get the expected value + :param t_type: which tomography type to use. Possible values: "linear_inv", "mle", "compressed_sensing", "lasso" + :param pauli_num: the number of pauli matrices to use in the tomography + :return: the density matrix as an ndarray + """ + + # if no pauli_num is specified use the maximum + if pauli_num==None: + pauli_num=len(Qubits) + + #Generate experiments + qubit_experiments = generate_state_tomography_experiment(program=program, qubits=qubits) + + exp_list = [] + #Experiment holds all 2^n possible pauli matrices for the given number of qubits + #Now take pauli_num random pauli matrices as per the paper's advice + if (pauli_num > len(qubit_experiments)): + print("Cannot sample more Pauli matrices thatn d^2!") + return None + exp_list = random.sample(list(qubit_experiments), pauli_num) + input_exp = PyQuilTomographyExperiment(settings=exp_list, program=program) + + #Group experiments if possible to minimize QPU runs + input_exp = group_experiments(input_exp) + + #NOTE: Change qvm depending on whether we are simulating qvm + qc.compiler.quil_to_native_quil = basic_compile + + results = list(measure_observables(qc=qc, tomo_experiment=input_exp, n_shots=num_shots)) + + if t_type == 'compressed_sensing': + return compressed_sensing_state_estimate(results=results, qubits=qubits) + elif t_type == 'mle': + return iterative_mle_state_estimate(results=results, qubits=qubits) + elif t_type == "linear_inv": + return linear_inv_state_estimate(results=results, qubits=qubits) + elif t_type == "lasso": + return lasso_state_estimate(results=results, qubits=qubits) +