Skip to content

Commit

Permalink
Operator tools (#140)
Browse files Browse the repository at this point in the history
This is a major module re-org plus adding new module to check plain old operators are unitary etc.

* move the random_operators file to the module operator_tools

* move the superoperator tools

* fix the legend not showing and add padding to the title

* superoperator_transformations

* notebook update

* type checking and docstring updates

* add superoperator_transformations

* tests

* superoperator_validation

* apply superop

* proj superops

* unused imports

* add and remove files

* validate operators

* improvements

* test superop transforms

* test apply superop

* test channel approx

* first pass at operator_tools refactor

* Spruce up permute_tesnor_factor and expand functionality.

* Move project_state_matrix_to_physical out of tomography.

* Fix import

* Update names and imports for project_state_matrix

* Remove nonsense code.

* compose kraus + tests

* improve

* big update

* format improvments

* update

* to fix the failed test add # NBVAL_IGNORE_OUTPUT

* it should have been # NBVAL_RAISES_EXCEPTION

* Typo fixing.

* matplotlib >= 3.1.0
  • Loading branch information
joshcombes authored Jun 14, 2019
1 parent 63805e9 commit a374475
Show file tree
Hide file tree
Showing 26 changed files with 2,433 additions and 406 deletions.
242 changes: 197 additions & 45 deletions examples/random_operators.ipynb

Large diffs are not rendered by default.

660 changes: 573 additions & 87 deletions examples/superoperator_tools.ipynb

Large diffs are not rendered by default.

6 changes: 3 additions & 3 deletions examples/tomography_state.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand All @@ -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]))"
]
},
Expand Down
8 changes: 8 additions & 0 deletions forest/benchmarking/operator_tools/__init__.py
Original file line number Diff line number Diff line change
@@ -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 *
82 changes: 82 additions & 0 deletions forest/benchmarking/operator_tools/apply_superoperator.py
Original file line number Diff line number Diff line change
@@ -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)
52 changes: 52 additions & 0 deletions forest/benchmarking/operator_tools/channel_approximation.py
Original file line number Diff line number Diff line change
@@ -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.
42 changes: 42 additions & 0 deletions forest/benchmarking/operator_tools/compose_superoperators.py
Original file line number Diff line number Diff line change
@@ -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]
52 changes: 52 additions & 0 deletions forest/benchmarking/operator_tools/project_state_matrix.py
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit a374475

Please sign in to comment.