Operator tools

examples/tomography_state.ipynb

@@ -265,10 +265,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from forest.benchmarking.tomography import project_density_matrix\n",
+    "from forest.benchmarking.operator_tools.project_state_matrix import project_state_matrix_to_physical\n",
     "rho_unphys = np.random.uniform(-1, 1, (2, 2)) \\\n",
     "    * np.exp(1.j * np.random.uniform(-np.pi, np.pi, (2, 2)))\n",
-    "rho_phys = project_density_matrix(rho_unphys)\n",
+    "rho_phys = project_state_matrix_to_physical(rho_unphys)\n",
     "\n",
     "fig, (ax1, ax2) = plt.subplots(1, 2)\n",
     "hinton(rho_unphys, ax=ax1)\n",
@@ -288,7 +288,7 @@
     "# https://doi.org/10.1103/PhysRevLett.108.070502\n",
     "\n",
     "eigs = np.diag(np.array(list(reversed([3.0/5, 1.0/2, 7.0/20, 1.0/10, -11.0/20]))))\n",
-    "phys = project_density_matrix(eigs)\n",
+    "phys = project_state_matrix_to_physical(eigs)\n",
     "np.allclose(phys, np.diag([0, 0, 1.0/5, 7.0/20, 9.0/20]))"
    ]
   },

forest/benchmarking/operator_tools/__init__.py

@@ -0,0 +1,8 @@
+from .apply_superoperator import *
+from .channel_approximation import *
+from .compose_superoperators import *
+from .project_superoperators import *
+from .superoperator_transformations import *
+from .random_operators import *
+from .validate_operator import *
+from .validate_superoperator import *

forest/benchmarking/operator_tools/apply_superoperator.py

@@ -0,0 +1,82 @@
+"""A module containing tools for applying superoperators to states.
+
+We have arbitrarily decided to use a column stacking convention.
+
+For more information about the conventions used, look at the file in
+/docs/Superoperator representations.md
+
+Further references include:
+
+[GRAPTN] Tensor networks and graphical calculus for open quantum systems
+         Wood et al.
+         Quant. Inf. Comp. 15, 0579-0811 (2015)
+         (no DOI)
+         https://arxiv.org/abs/1111.6950
+
+[MATQO] On the Matrix Representation of Quantum Operations
+        Nambu et al.,
+        arXiv: 0504091 (2005)
+        https://arxiv.org/abs/quant-ph/0504091
+
+[DUAL] On duality between quantum maps and quantum states
+       Zyczkowski et al.,
+       Open Syst. Inf. Dyn. 11, 3 (2004)
+       https://dx.doi.org/10.1023/B:OPSY.0000024753.05661.c2
+       https://arxiv.org/abs/quant-ph/0401119
+
+"""
+from typing import Sequence
+import numpy as np
+from forest.benchmarking.utils import partial_trace
+
+
+def apply_kraus_ops_2_state(kraus_ops: Sequence[np.ndarray], state: np.ndarray) -> np.ndarray:
+    r"""
+    Apply a quantum channel, specified by Kraus operators, to state.
+
+    The Kraus operators need not be square.
+
+    :param kraus_ops: A list or tuple of N Kraus operators, each operator is M by dim ndarray
+    :param state: A dim by dim ndarray which is the density matrix for the state
+    :return: M by M ndarray which is the density matrix for the state after the action of kraus_ops
+    """
+    if isinstance(kraus_ops, np.ndarray):  # handle input of single kraus op
+        if len(kraus_ops[0].shape) < 2:  # first elem is not a matrix
+            kraus_ops = [kraus_ops]
+
+    dim, _ = state.shape
+    rows, cols = kraus_ops[0].shape
+
+    if dim != cols:
+        raise ValueError("Dimensions of state and Kraus operator are incompatible")
+
+    new_state = np.zeros((rows, rows))
+    for M in kraus_ops:
+        new_state += M @ state @ np.transpose(M.conj())
+
+    return new_state
+
+
+def apply_choi_matrix_2_state(choi: np.ndarray, state: np.ndarray) -> np.ndarray:
+    r"""
+    Apply a quantum channel, specified by a Choi matrix (using the column stacking convention),
+    to a state.
+
+    The Choi matrix is a dim**2 by dim**2 matrix and the state rho is a dim by dim matrix. The
+    output state is
+
+    rho_{out} = Tr_{A_{in}}[(rho^T \otimes Id) Choi_matrix ],
+
+    where T denotes transposition and Tr_{A_{in}} is the partial trace over input Hilbert space H_{
+    A_{in}}; the Choi matrix representing a process mapping rho in H_{A_{in}} to rho_{out}
+    in H_{B_{out}} is regarded as an operator on the space H_{A_{in}} \otimes H_{B_{out}}.
+
+
+    :param choi: a dim**2 by dim**2 matrix
+    :param state: A dim by dim ndarray which is the density matrix for the state
+    :return: a dim by dim matrix.
+    """
+    dim = int(np.sqrt(np.asarray(choi).shape[0]))
+    dims = [dim, dim]
+    tot_matrix = np.kron(state.transpose(), np.identity(dim)) @ choi
+    return partial_trace(tot_matrix, [1], dims)

forest/benchmarking/operator_tools/channel_approximation.py

@@ -0,0 +1,52 @@
+"""A module containing tools for approximating channels.
+
+We have arbitrarily decided to use a column stacking convention.
+
+For more information about the conventions used, look at the file in
+/docs/Superoperator representations.md
+
+Further references include:
+
+[GRAPTN] Tensor networks and graphical calculus for open quantum systems
+         Wood et al.
+         Quant. Inf. Comp. 15, 0579-0811 (2015)
+         (no DOI)
+         https://arxiv.org/abs/1111.6950
+
+[MATQO] On the Matrix Representation of Quantum Operations
+        Nambu et al.,
+        arXiv: 0504091 (2005)
+        https://arxiv.org/abs/quant-ph/0504091
+
+[DUAL] On duality between quantum maps and quantum states
+       Zyczkowski et al.,
+       Open Syst. Inf. Dyn. 11, 3 (2004)
+       https://dx.doi.org/10.1023/B:OPSY.0000024753.05661.c2
+       https://arxiv.org/abs/quant-ph/0401119
+
+"""
+import numpy as np
+
+
+def pauli_twirl_chi_matrix(chi_matrix: np.ndarray) -> np.ndarray:
+    r"""
+    Implements a Pauli twirl of a chi matrix (aka a process matrix).
+
+    See the folloiwng reference for more details
+
+    [SPICC] Scalable protocol for identification of correctable codes
+            Silva et al.,
+            PRA 78, 012347 
2008
+            http://dx.doi.org/10.1103/PhysRevA.78.012347
+            https://arxiv.org/abs/0710.1900
+
+    Note: Pauli twirling a quantum channel can give rise to a channel that is less noisy; use with
+    care.
+
+    :param chi_matrix:  a dim**2 by dim**2 chi or process matrix
+    :return: dim**2 by dim**2 chi or process matrix
+    """
+    return np.diag(chi_matrix.diagonal())
+
+
+# TODO: Honest approximations for Channels that act on one or MORE qubits.

forest/benchmarking/operator_tools/compose_superoperators.py

@@ -0,0 +1,42 @@
+"""A module containing tools for composing superoperators.
+"""
+from typing import Sequence
+import numpy as np
+
+
+def tensor_channel_kraus(k2: Sequence[np.ndarray],
+                         k1: Sequence[np.ndarray]) -> Sequence[np.ndarray]:
+    r"""
+    Given the Kraus representaion for two channels, $\mathcal E$ and $\mathcal F$, acting on
+    different systems this function returns the Kraus operators representing the composition of
+    these independent channels.
+
+    Suppose $\mathcal E$ and $\mathcal F$ both have one Kraus operator K_1 = X and K_2 = H,
+    that is they are unitary. Then, with respect to the tensor product structure
+
+            $H_2 \otimes H_1$
+
+    of the individual systems this function returns $K_{\rm tot} = H \otimes X$.
+
+    :param k1: The list of Kraus operators on the first system.
+    :param k2: The list of Kraus operators on the second system.
+    :return: A list of tensored Kraus operators.
+    """
+    # TODO: make this function work for an arbitrary number of Kraus operators
+    return [np.kron(k2l, k1j) for k1j in k1 for k2l in k2]
+
+
+def compose_channel_kraus(k2: Sequence[np.ndarray],
+                          k1: Sequence[np.ndarray]) -> Sequence[np.ndarray]:
+    r"""
+    Given two channels, K_1 and K_2, acting on the same system in the Kraus representation this
+    function return the Kraus operators representing the composition of the channels.
+
+    It is assumed that K_1 is applied first then K_2 is applied.
+
+    :param k2: The list of Kraus operators that are applied second.
+    :param k1: The list of Kraus operators that are applied first.
+    :return: A combinatorially generated list of composed Kraus operators.
+    """
+    # TODO: make this function work for an arbitrary number of Kraus operators
+    return [np.dot(k2l, k1j) for k1j in k1 for k2l in k2]

forest/benchmarking/operator_tools/project_state_matrix.py

@@ -0,0 +1,52 @@
+import numpy as np
+import functools
+from scipy.linalg import eigh
+
+
+def project_state_matrix_to_physical(rho: np.ndarray) -> np.ndarray:
+    """
+    Project a possibly unphysical estimated density matrix to the closest (with respect to the
+    2-norm) positive semi-definite matrix with trace 1, that is a valid quantum state.
+
+    This is the so called "wizard" method. It is described in the following reference:
+
+    [MLEWIZ] Efficient Method for Computing the Maximum-Likelihood Quantum State from
+             Measurements with Additive Gaussian Noise
+             Smolin et al.,
+             Phys. Rev. Lett. 108, 070502 (2012)
+             https://doi.org/10.1103/PhysRevLett.108.070502
+             https://arxiv.org/abs/1106.5458
+
+    :param rho: the density (state) matrix with shape (N, N)
+    :return rho_projected: The closest positive semi-definite trace 1 matrix to rho.
+    """
+    # Rescale to trace 1 if the matrix is not already
+    rho_impure = rho / np.trace(rho)
+
+    dimension = rho_impure.shape[0]  # the dimension of the Hilbert space
+    [eigvals, eigvecs] = eigh(rho_impure)
+
+    # If matrix is already trace one PSD, we are done
+    if np.min(eigvals) >= 0:
+        return rho_impure
+
+    # Otherwise, continue finding closest trace one, PSD matrix
+    eigvals = list(eigvals)
+    eigvals.reverse()
+    eigvals_new = [0.0] * len(eigvals)
+
+    i = dimension
+    accumulator = 0.0  # Accumulator
+    while eigvals[i - 1] + accumulator / float(i) < 0:
+        accumulator += eigvals[i - 1]
+        i -= 1
+    for j in range(i):
+        eigvals_new[j] = eigvals[j] + accumulator / float(i)
+    eigvals_new.reverse()
+
+    # Reconstruct the matrix
+    rho_projected = functools.reduce(np.dot, (eigvecs,
+                                              np.diag(eigvals_new),
+                                              np.conj(eigvecs.T)))
+
+    return rho_projected