Repository: quantumlib/OpenFermion
Pull request #88 (2017-10-27 22:17:30)

Adds basic single qubit error channels for simulation convenience
Adds the ability to apply simple 1-qubit error channels numerically to given density matrices with tests.  Useful for testing error resilience in VQE and QSE methods.

jarrodmcc: Comments should be taken care of

Patch:
diff --git a/src/openfermion/utils/__init__.py b/src/openfermion/utils/__init__.py
index d4476a6..9d5bb20 100644
--- a/src/openfermion/utils/__init__.py
+++ b/src/openfermion/utils/__init__.py
@@ -10,6 +10,9 @@
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
 
+from ._channel_state import (amplitude_damping_channel, dephasing_channel,
+                             depolarizing_channel)
+
 from ._grid import Grid
 
 from ._operator_utils import (commutator, count_qubits,
diff --git a/src/openfermion/utils/_channel_state.py b/src/openfermion/utils/_channel_state.py
new file mode 100644
index 0000000..2cdd291
--- /dev/null
+++ b/src/openfermion/utils/_channel_state.py
@@ -0,0 +1,184 @@
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+"""Module to manipulate basic models of quantum channels"""
+
+from functools import reduce
+from itertools import chain
+from numpy import array, dot, eye, kron, log2, sqrt
+
+
+def _verify_channel_inputs(density_matrix, probability, target_qubit):
+    """Verifies input parameters for channels
+
+    Args:
+        density_matrix (numpy.ndarray): Density matrix of the system
+        probability (float): Probability error is applied p \in [0, 1]
+        target_qubit (int): target for the channel error.
+
+    Returns:
+        new_density_matrix(numpy.ndarray): Density matrix with the channel
+            applied.
+    """
+    n_qubits = int(log2(density_matrix.shape[0]))
+
+    if (len(density_matrix.shape) != 2 or
+            density_matrix.shape[0] != density_matrix.shape[1]):
+        raise ValueError("Error in input of density matrix to channel.")
+    if (probability < 0) or (probability > 1):
+        raise ValueError("Channel probability must be between 0 and 1.")
+    if (target_qubit < 0) or (target_qubit >= n_qubits):
+        raise ValueError("Target qubits must be within number of qubits.")
+
+
+def _lift_operator(operator, n_qubits, target_qubit):
+    """Lift a single qubit operator into the n_qubit space by kron product
+
+    Args:
+        operator (ndarray): Single qubit operator to lift into full space
+        n_qubits (int): Number of total qubits in the space
+        target_qubit (int): Qubit to act on
+
+    Return:
+        new_operator(Sparse Operator): Operator representing the embedding in
+            the full space.
+    """
+    new_operator = (
+        reduce(kron,
+               chain(
+                   (eye(2) for i in range(0, target_qubit)),
+                   [operator],
+                   (eye(2) for i in range(target_qubit + 1, n_qubits)))))
+    return new_operator
+
+
+def amplitude_damping_channel(density_matrix, probability, target_qubit,
+                              transpose=False):
+    """Apply an amplitude damping channel
+
+    Applies an amplitude damping channel with a given probability to the target
+    qubit in the density_matrix.
+
+    Args:
+        density_matrix (numpy.ndarray): Density matrix of the system
+        probability (float): Probability error is applied p \in [0, 1]
+        target_qubit (int): target for the channel error.
+        transpose (bool): Transpose channel operators, useful for acting on
+            Hamiltonians in variational channel state models
+
+    Returns:
+        new_density_matrix(numpy.ndarray): Density matrix with the channel
+            applied.
+    """
+    _verify_channel_inputs(density_matrix, probability, target_qubit)
+    n_qubits = int(log2(density_matrix.shape[0]))
+
+    E0 = _lift_operator(array([[1.0, 0.0],
+                               [0.0, sqrt(1.0 - probability)]], dtype=complex),
+                        n_qubits, target_qubit)
+
+    E1 = _lift_operator(array([[0.0, sqrt(probability)],
+                               [0.0, 0.0]], dtype=complex),
+                        n_qubits, target_qubit)
+
+    if transpose:
+        E0 = E0.T
+        E1 = E1.T
+
+    new_density_matrix = (dot(E0, dot(density_matrix, E0.T)) +
+                          dot(E1, dot(density_matrix, E1.T)))
+
+    return new_density_matrix
+
+
+def dephasing_channel(density_matrix, probability, target_qubit,
+                      transpose=False):
+    """Apply a dephasing channel
+
+    Applies an amplitude damping channel with a given probability to the target
+    qubit in the density_matrix.
+
+    Args:
+        density_matrix (numpy.ndarray): Density matrix of the system
+        probability (float): Probability error is applied p \in [0, 1]
+        target_qubit (int): target for the channel error.
+        transpose (bool): Transpose channel operators, useful for acting on
+            Hamiltonians in variational channel state models
+
+    Returns:
+        new_density_matrix (numpy.ndarray): Density matrix with the channel
+            applied.
+    """
+    _verify_channel_inputs(density_matrix, probability, target_qubit)
+    n_qubits = int(log2(density_matrix.shape[0]))
+
+    E0 = _lift_operator(sqrt(1.0 - probability/2.) * eye(2),
+                        n_qubits, target_qubit)
+    E1 = _lift_operator(sqrt(probability/2.) *
+                        array([[1.0, 0.0], [1.0, -1.0]]),
+                        n_qubits, target_qubit)
+
+    if transpose:
+        E0 = E0.T
+        E1 = E1.T
+
+    new_density_matrix = (dot(E0, dot(density_matrix, E0.T)) +
+                          dot(E1, dot(density_matrix, E1.T)))
+
+    return new_density_matrix
+
+
+def depolarizing_channel(density_matrix, probability, target_qubit,
+                         transpose=False):
+    """Apply a depolarizing channel
+
+    Applies an amplitude damping channel with a given probability to the target
+    qubit in the density_matrix.
+
+    Args:
+        density_matrix (numpy.ndarray): Density matrix of the system
+        probability (float): Probability error is applied p \in [0, 1]
+        target_qubit (int): target for the channel error.
+        transpose (bool): Transpose channel operators, useful for acting on
+            Hamiltonians in variational channel state models
+
+    Returns:
+        new_density_matrix (numpy.ndarray): Density matrix with the channel
+            applied.
+    """
+    _verify_channel_inputs(density_matrix, probability, target_qubit)
+    n_qubits = int(log2(density_matrix.shape[0]))
+
+    E0 = _lift_operator(sqrt(1.0 - probability) * eye(2),
+                        n_qubits, target_qubit)
+    E1 = _lift_operator(sqrt(probability/3.) * array([[0.0, 1.0],
+                                                      [1.0, 0.0]]),
+                        n_qubits, target_qubit)
+    E2 = _lift_operator(sqrt(probability/3.) * array([[0.0, -1.0j],
+                                                      [1.0j, 0.0]]),
+                        n_qubits, target_qubit)
+    E3 = _lift_operator(sqrt(probability/3.) * array([[1.0, 0.0],
+                                                      [0.0, -1.0]]),
+                        n_qubits, target_qubit)
+
+    if transpose:
+        E0 = E0.T
+        E1 = E1.T
+        E2 = E2.T
+        E3 = E3.T
+
+    new_density_matrix = (dot(E0, dot(density_matrix, E0.T)) +
+                          dot(E1, dot(density_matrix, E1.T)) +
+                          dot(E2, dot(density_matrix, E2.T)) +
+                          dot(E3, dot(density_matrix, E3.T)))
+
+    return new_density_matrix


Test patch:
diff --git a/src/openfermion/utils/_channel_state_test.py b/src/openfermion/utils/_channel_state_test.py
new file mode 100644
index 0000000..aca937f
--- /dev/null
+++ b/src/openfermion/utils/_channel_state_test.py
@@ -0,0 +1,148 @@
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+from __future__ import absolute_import
+
+from numpy import dot, zeros
+from scipy.linalg import norm
+
+import unittest
+
+from ._channel_state import *
+
+
+class ChannelTest(unittest.TestCase):
+
+    def setUp(self):
+        """Initialize a few density matrices"""
+        zero_state = array([[1], [0]], dtype=complex)
+        one_state = array([[0], [1]], dtype=complex)
+        one_one_state = kron(one_state, one_state)
+        zero_zero_state = kron(zero_state, zero_state)
+        cat_state = 1./sqrt(2) * (zero_zero_state + one_one_state)
+
+        self.density_matrix = dot(one_one_state, one_one_state.T)
+        self.cat_matrix = dot(cat_state, cat_state.T)
+
+    def test_amplitude_damping(self):
+        """Test amplitude damping on a simple qubit state"""
+
+        # With probability 0
+        test_density_matrix = (
+            amplitude_damping_channel(self.density_matrix, 0, 1))
+        self.assertAlmostEquals(norm(self.density_matrix -
+                                     test_density_matrix), 0.0)
+
+        test_density_matrix = (
+            amplitude_damping_channel(self.density_matrix, 0, 1,
+                                      transpose=True))
+        self.assertAlmostEquals(norm(self.density_matrix -
+                                     test_density_matrix), 0.0)
+
+        # With probability 1
+        correct_density_matrix = zeros((4, 4), dtype=complex)
+        correct_density_matrix[2, 2] = 1
+
+        test_density_matrix = (
+            amplitude_damping_channel(self.density_matrix, 1, 1))
+
+        self.assertAlmostEquals(norm(correct_density_matrix -
+                                     test_density_matrix), 0.0)
+
+    def test_dephasing(self):
+        """Test dephasing on a simple qubit state"""
+
+        # Check for identity on |11> state
+        test_density_matrix = (
+            dephasing_channel(self.density_matrix, 1, 1))
+        self.assertAlmostEquals(norm(self.density_matrix -
+                                     test_density_matrix), 0.0)
+
+        test_density_matrix = (
+            dephasing_channel(self.density_matrix, 1, 1,
+                              transpose=True))
+
+        correct_matrix = array([[0., 0., 0., 0.],
+                                [0., 0., 0., 0.],
+                                [0., 0., 0.5, -0.5],
+                                [0., 0., -0.5, 1.]])
+        self.assertAlmostEquals(norm(correct_matrix -
+                                     test_density_matrix), 0.0)
+
+        # Check for correct action on cat state
+        # With probability = 0
+        test_density_matrix = (
+            dephasing_channel(self.cat_matrix, 0, 1))
+        self.assertAlmostEquals(norm(self.cat_matrix -
+                                     test_density_matrix), 0.0)
+        # With probability = 1
+
+        correct_matrix = array([[0.50,  0.25,  0.00,  0.00],
+                                [0.25,  0.25, 0.00, -0.25],
+                                [0.00,  0.00,  0.00,  0.00],
+                                [0.00, -0.25, 0.00, 0.50]])
+        test_density_matrix = (
+            dephasing_channel(self.cat_matrix, 1, 1))
+        self.assertAlmostEquals(norm(correct_matrix -
+                                     test_density_matrix), 0.0)
+
+    def test_depolarizing(self):
+        """Test depolarizing on a simple qubit state"""
+
+        # With probability = 0
+        test_density_matrix = (
+            depolarizing_channel(self.cat_matrix, 0, 1))
+        self.assertAlmostEquals(norm(self.cat_matrix -
+                                     test_density_matrix), 0.0)
+
+        test_density_matrix = (
+            depolarizing_channel(self.cat_matrix, 0, 1,
+                                 transpose=True))
+        self.assertAlmostEquals(norm(self.cat_matrix -
+                                     test_density_matrix), 0.0)
+
+        # With probability 1 on both qubits
+        correct_density_matrix = array([[0.0555555, 0.,  0.0,  0.27777778],
+                                        [0.00, 0.00, -0.22222222, 0.000],
+                                        [0.00, -0.22222222, 0., 0.],
+                                        [0.27777778, 0.0, 0.0, 0.05555556]])
+
+        test_density_matrix = (
+            depolarizing_channel(self.cat_matrix, 1, 0))
+        test_density_matrix = (
+            depolarizing_channel(test_density_matrix, 1, 1))
+
+        self.assertAlmostEquals(norm(correct_density_matrix -
+                                     test_density_matrix), 0.0, places=6)
+
+        # Depolarizing channel should be self-adjoint
+        test_density_matrix = (
+            depolarizing_channel(self.cat_matrix, 1, 0,
+                                 transpose=True))
+        test_density_matrix = (
+            depolarizing_channel(test_density_matrix, 1, 1,
+                                 transpose=True))
+
+        self.assertAlmostEquals(norm(correct_density_matrix -
+                                     test_density_matrix), 0.0, places=6)
+
+    def test_verification(self):
+        """Verify basic sanity checking on inputs"""
+        with self.assertRaises(ValueError):
+            _ = amplitude_damping_channel(self.density_matrix, 2, 1)
+
+        with self.assertRaises(ValueError):
+            _ = amplitude_damping_channel(self.density_matrix, 0.5, 3)
+
+        with self.assertRaises(ValueError):
+            bad_density = zeros((3, 4))
+            _ = amplitude_damping_channel(bad_density, 0.5, 3)