Measurement estimation (#139)

* Estimating sets of pauli sums to fixed variance

Estimates the expected value of a linear combination of pauli terms to a
fixed variance on the estimator.  The variance of the pauli sum is
estimated by sample variance.  Example notebooks will be added later.

grove/measurements/estimation.py

@@ -0,0 +1,177 @@
+"""
+Utilities for estimating expected values of Pauli terms given pyquil programs
+"""
+import numpy as np
+from pyquil.paulis import (PauliSum, PauliTerm, commuting_sets, sI,
+                           term_with_coeff)
+from pyquil.quil import Program
+from pyquil.gates import RY, RX
+
+
+class CommutationError(Exception):
+    pass
+
+
+def remove_imaginary_terms(pauli_sums):
+    """
+    Remove the imaginary component of each term in a Pauli sum
+
+    :param PauliSum pauli_sums: The Pauli sum to process.
+    :return: a purely hermitian Pauli sum.
+    :rtype: PauliSum
+    """
+    if not isinstance(pauli_sums, PauliSum):
+        raise TypeError("not a pauli sum. please give me one")
+    new_term = sI(0) * 0.0
+    for term in pauli_sums:
+        new_term += term_with_coeff(term, term.coefficient.real)
+
+    return new_term
+
+
+def get_rotation_program(pauli_term):
+    """
+    Generate a rotation program so that the pauli term is diagonal
+
+    :param PauliTerm pauli_term: The Pauli term used to generate diagonalizing
+                                 one-qubit rotations.
+    :return: The rotation program.
+    :rtype: Program
+    """
+    meas_basis_change = Program()
+    for index, gate in pauli_term:
+        if gate == 'X':
+            meas_basis_change.inst(RY(-np.pi / 2, index))
+        elif gate == 'Y':
+            meas_basis_change.inst(RX(np.pi / 2, index))
+        elif gate == 'Z':
+            pass
+        else:
+            raise ValueError()
+
+    return meas_basis_change
+
+
+def get_parity(pauli_terms, bitstring_results):
+    """
+    Calculate the eigenvalues of Pauli operators given results of projective measurements
+
+    The single-qubit projective measurement results (elements of
+    `bitstring_results`) are listed in physical-qubit-label numerical order.
+
+    An example:
+
+    Consider a Pauli term Z1 Z5 Z6 I7 and a collection of single-qubit
+    measurement results corresponding to measurements in the z-basis on qubits
+    {1, 5, 6, 7}. Each element of bitstring_results is an element of
+    :math:`\{0, 1\}^{\otimes 4}`.  If [0, 0, 1, 0] and [1, 0, 1, 1]
+    are the two projective measurement results in `bitstring_results` then
+    this method returns a 1 x 2 numpy array with values [[-1, 1]]
+
+    :param List pauli_terms: A list of Pauli terms operators to use
+    :param bitstring_results: A list of projective measurement results.  Each
+                              element is a list of single-qubit measurements.
+    :return: Array (m x n) of {+1, -1} eigenvalues for the m-operators in
+             `pauli_terms` associated with the n measurement results.
+    :rtype: np.ndarray
+    """
+    max_weight_pauli_index = np.argmax(
+        [len(term.get_qubits()) for term in pauli_terms])
+    active_qubit_indices = sorted(
+        pauli_terms[max_weight_pauli_index]._ops.keys())
+    index_mapper = dict(zip(active_qubit_indices,
+                            range(len(active_qubit_indices))))
+
+    results = np.zeros((len(pauli_terms), len(bitstring_results)))
+
+    # convert to array so we can fancy index into it later.
+    # list() is needed to cast because we don't really want a map object
+    bitstring_results = list(map(np.array, bitstring_results))
+    for row_idx, term in enumerate(pauli_terms):
+        memory_index = np.array(list(map(lambda x: index_mapper[x],
+                                         sorted(term.get_qubits()))))
+
+        results[row_idx, :] = [-2 * (sum(x[memory_index]) % 2) +
+                               1 for x in bitstring_results]
+    return results
+
+
+def estimate_pauli_sum(pauli_terms, basis_transform_dict, program,
+                       variance_bound, quantum_resource,
+                       commutation_check=True):
+    """
+    Estimate the mean of a sum of pauli terms to set variance
+
+    The sample variance is calculated by
+
+    .. math::
+        \begin{align}
+        \mathrm{Var}[\hat{\langle H \rangle}] = \sum_{i, j}h_{i}h_{j}
+        \mathrm{Cov}(\hat{\langle P_{i} \rangle}, \hat{\langle P_{j} \rangle})
+        \end{align}
+
+    :param pauli_terms: list of pauli terms to measure simultaneously or a
+                        PauliSum object
+    :param basis_transform_dict: basis transform dictionary where the key is
+                                 the qubit index and the value is the basis to
+                                 rotate into. Valid basis is [I, X, Y, Z].
+    :param program: program generating a state to sample from
+    :param variance_bound:  Bound on the variance of the estimator for the
+                            PauliSum. Remember this is the SQUARE of the
+                            standard error!
+    :param quantum_resource: quantum abstract machine object
+    :param Bool commutation_check: Optional flag toggling a safety check
+                                   ensuring all terms in `pauli_terms`
+                                   commute with each other
+    :return: estimated expected value, covariance matrix, variance of the
+             estimator, and the number of shots taken
+    """
+    if not isinstance(pauli_terms, (list, PauliSum)):
+        raise TypeError("pauli_terms needs to be a list or a PauliSum")
+
+    if isinstance(pauli_terms, PauliSum):
+        pauli_terms = PauliSum.terms
+
+    # check if each term commutes with everything
+    if commutation_check:
+        if len(commuting_sets(sum(pauli_terms))) != 1:
+            raise CommutationError("Not all terms commute in the expected way")
+
+    pauli_for_rotations = PauliTerm.from_list(
+        [(value, key) for key, value in basis_transform_dict.items()])
+    post_rotations = get_rotation_program(pauli_for_rotations)
+
+    coeff_vec = np.array(
+        list(map(lambda x: x.coefficient, pauli_terms))).reshape((-1, 1))
+
+    # upper bound on samples given by IV of arXiv:1801.03524
+    num_sample_ubound = int(np.ceil(np.sum(np.abs(coeff_vec))**2 / variance_bound))
+    results = None
+    sample_variance = np.infty
+    number_of_samples = 0
+    while (sample_variance > variance_bound and
+           number_of_samples < num_sample_ubound):
+        # note: bit string values are returned according to their listed order
+        # in run_and_measure.  Therefore, sort beforehand to keep alpha numeric
+        tresults = quantum_resource.run_and_measure(
+            program + post_rotations,
+            sorted(list(basis_transform_dict.keys())),
+            trials=min(10000, num_sample_ubound))
+        number_of_samples += len(tresults)
+
+        parity_results = get_parity(pauli_terms, tresults)
+
+        # Note: easy improvement would be to update mean and variance on the fly
+        # instead of storing all these results.
+        if results is None:
+            results = parity_results
+        else:
+            results = np.hstack((results, parity_results))
+
+        # calculate the expected values....
+        covariance_mat = np.cov(results)
+        sample_variance = coeff_vec.T.dot(covariance_mat).dot(coeff_vec) / \
+                          results.shape[1]
+
+    return coeff_vec.T.dot(np.mean(results, axis=1)), covariance_mat, \
+           sample_variance, results.shape[1]

grove/tests/measurements/test_estimation.py

@@ -0,0 +1,148 @@
+"""
+Tests for the estimation module
+"""
+import pytest
+from mock import Mock
+import numpy as np
+from scipy.stats import bernoulli
+from pyquil.paulis import sX, sY, sZ, sI, PauliSum, is_zero
+from pyquil.quil import Program
+from pyquil.gates import RY, RX
+from pyquil.api import QVMConnection
+from grove.measurements.estimation import (remove_imaginary_terms,
+                                           get_rotation_program,
+                                           get_parity,
+                                           estimate_pauli_sum,
+                                           CommutationError)
+
+
+def test_imaginary_removal():
+    """
+    remove terms with imaginary coefficients from a pauli sum
+    """
+    test_term = 0.25 * sX(1) * sZ(2) * sX(3) + 0.25j * sX(1) * sZ(2) * sY(3)
+    test_term += -0.25j * sY(1) * sZ(2) * sX(3) + 0.25 * sY(1) * sZ(2) * sY(3)
+    true_term = 0.25 * sX(1) * sZ(2) * sX(3) + 0.25 * sY(1) * sZ(2) * sY(3)
+    assert remove_imaginary_terms(test_term) == true_term
+
+    test_term = (0.25 + 1j) * sX(0) * sZ(2) + 1j * sZ(2)
+    # is_identity in pyquil apparently thinks zero is identity
+    assert remove_imaginary_terms(test_term) == 0.25 * sX(0) * sZ(2)
+
+    test_term = 0.25 * sX(0) * sZ(2) + 1j * sZ(2)
+    assert remove_imaginary_terms(test_term) == PauliSum([0.25 * sX(0) * sZ(2)])
+
+    with pytest.raises(TypeError):
+        remove_imaginary_terms(5)
+
+    with pytest.raises(TypeError):
+        remove_imaginary_terms(sX(0))
+
+
+def test_rotation_programs():
+    """
+    Testing the generation of post rotations
+    """
+    test_term = sZ(0) * sX(20) * sI(100) * sY(5)
+    # note: string comparison of programs requires gates to be in the same order
+    true_rotation_program = Program().inst(
+        [RX(np.pi / 2)(5), RY(-np.pi / 2)(20)])
+    test_rotation_program = get_rotation_program(test_term)
+    assert true_rotation_program.out() == test_rotation_program.out()
+
+
+def test_get_parity():
+    """
+    Check if our way to compute parity is correct
+    """
+    single_qubit_results = [[0]] * 50 + [[1]] * 50
+    single_qubit_parity_results = list(map(lambda x: -2 * x[0] + 1,
+                                           single_qubit_results))
+
+    # just making sure I constructed my test properly
+    assert np.allclose(np.array([1] * 50 + [-1] * 50),
+                       single_qubit_parity_results)
+
+    test_results = get_parity([sZ(5)], single_qubit_results)
+    assert np.allclose(single_qubit_parity_results, test_results[0, :])
+
+    np.random.seed(87655678)
+    brv1 = bernoulli(p=0.25)
+    brv2 = bernoulli(p=0.4)
+    n = 500
+    two_qubit_measurements = list(zip(brv1.rvs(size=n), brv2.rvs(size=n)))
+    pauli_terms = [sZ(0), sZ(1), sZ(0) * sZ(1)]
+    parity_results = np.zeros((len(pauli_terms), n))
+    parity_results[0, :] = [-2 * x[0] + 1 for x in two_qubit_measurements]
+    parity_results[1, :] = [-2 * x[1] + 1 for x in two_qubit_measurements]
+    parity_results[2, :] = [-2 * (sum(x) % 2) + 1 for x in
+                            two_qubit_measurements]
+
+    test_parity_results = get_parity(pauli_terms, two_qubit_measurements)
+    assert np.allclose(test_parity_results, parity_results)
+
+
+def test_estimate_pauli_sum():
+    """
+    Full test of the estimation procedures
+    """
+    quantum_resource = QVMConnection()
+
+    # type checks
+    with pytest.raises(TypeError):
+        estimate_pauli_sum('5', {0: 'X', 1: 'Z'}, Program(), 1.0E-3,
+                           quantum_resource)
+
+    with pytest.raises(CommutationError):
+        estimate_pauli_sum([sX(0), sY(0)], {0: 'X', 1: 'Z'}, Program(), 1.0E-3,
+                           quantum_resource)
+
+    with pytest.raises(TypeError):
+        estimate_pauli_sum(sX(0), {0: 'X', 1: 'Z'}, Program(), 1.0E-3,
+                           quantum_resource)
+
+    # mock out qvm
+    np.random.seed(87655678)
+    brv1 = bernoulli(p=0.25)
+    brv2 = bernoulli(p=0.4)
+    n = 500
+    two_qubit_measurements = list(zip(brv1.rvs(size=n), brv2.rvs(size=n)))
+    pauli_terms = [sZ(0), sZ(1), sZ(0) * sZ(1)]
+
+    fakeQVM = Mock(spec=QVMConnection())
+    fakeQVM.run_and_measure = Mock(return_value=two_qubit_measurements)
+    mean, cov, estimator_var, shots = estimate_pauli_sum(pauli_terms,
+                                                         {0: 'Z', 1: 'Z'},
+                                                         Program(),
+                                                         1.0E-1, fakeQVM)
+    parity_results = np.zeros((len(pauli_terms), n))
+    parity_results[0, :] = [-2 * x[0] + 1 for x in two_qubit_measurements]
+    parity_results[1, :] = [-2 * x[1] + 1 for x in two_qubit_measurements]
+    parity_results[2, :] = [-2 * (sum(x) % 2) + 1 for x in
+                            two_qubit_measurements]
+
+    assert np.allclose(np.cov(parity_results), cov)
+    assert np.isclose(np.sum(np.mean(parity_results, axis=1)), mean)
+    assert np.isclose(shots, n)
+    variance_to_beat = np.sum(cov) / n
+    assert np.isclose(variance_to_beat, estimator_var)
+
+    # Double the shots by ever so slightly decreasing variance bound
+    double_two_q_measurements = two_qubit_measurements + two_qubit_measurements
+    mean, cov, estimator_var, shots = estimate_pauli_sum(pauli_terms,
+                                                         {0: 'Z', 1: 'Z'},
+                                                         Program(),
+                                                         variance_to_beat - \
+                                                         1.0E-8, fakeQVM)
+
+    parity_results = np.zeros((len(pauli_terms), 2 * n))
+    parity_results[0, :] = [-2 * x[0] + 1 for x in double_two_q_measurements]
+    parity_results[1, :] = [-2 * x[1] + 1 for x in double_two_q_measurements]
+    parity_results[2, :] = [-2 * (sum(x) % 2) + 1 for x in
+                            double_two_q_measurements]
+
+    assert np.allclose(np.cov(parity_results), cov)
+    assert np.isclose(np.sum(np.mean(parity_results, axis=1)), mean)
+    assert np.isclose(shots, 2 * n)
+    assert np.isclose(np.sum(cov) / (2 * n), estimator_var)
+

grove/tests/measurements/test_term_grouping.py

@@ -146,3 +146,4 @@ def test_term_grouping():
                 ('Y', 'Y', 'Y', 'Y'): set(map(lambda x: x.id(), yyyy_terms))}
     for key, value in clumped_terms.items():
         assert set(map(lambda x: x.id(), clumped_terms[key])) == true_set[key]
+