QuTiP Example: Stinespring Dilation of Quantum Channels

Christopher Granade
University of Sydney $\newcommand{\ket}[1]{\left|#1\right\rangle}$ $\newcommand{\bra}[1]{\left\langle#1\right|}$ $\newcommand{\cnot}{{\scriptstyle \rm CNOT}}$ $\newcommand{\TT}{\operatorname{T}}$ $\newcommand{\Tr}{\operatorname{Tr}}$ $\newcommand{\Uni}{\operatorname{U}}$ $\newcommand{\Chan}{\operatorname{C}}$ $\newcommand{\Hil}{\mathcal{H}}$

Introduction

In this notebook, we will demonstrate the to_stinespring function, which converts a Qobj representing the quantum channel $\Phi$ to a pair of partial isometries that describe the action of $\Phi$.

In introducing the Stinespring dilation, it is helpful to first adopt some notation. Let $\Hil_X$, $\Hil_Y$ and $\Hil_Z$ be finite-dimensional Hilbert spaces. Let $\Uni(\Hil_X, \Hil_Y)$ be the set of partial isometries from $\Hil_X$ to $\Hil_Y$, and $\TT(\Hil_X, \Hil_Y)$ be the set of maps (that is, superoperators) from operators on $\Hil_X$ to those acting on $\Hil_Y$. Finally, let $\Chan(\Hil_X, \Hil_Y) \subset \TT(\Hil_X, \Hil_Y)$ be the set of channels; that is, those maps which are completely positive and trace-preserving.

For a quantum map $\Phi \in \TT(\Hil_X, \Hil_Y)$, then, the Stinespring dilation of $\Phi$ is a pair of partial isometries $A, B \in \Uni(\Hil_X, \Hil_Y \otimes \Hil_Z)$ for some space $\Hil_Z$ such that \begin{equation} \Phi(X) = \Tr_Z (A X B^\dagger). \label{eq:stinespring-action} \end{equation} If $\Phi$ is completely positive, then $A = B$. We do not insist on complete positivity, however, as it is common to consider the Stinespring dilation for the difference between two quantum channels, $\Phi - \Phi'$. By construction, such a map is not completely positive, nor is it trace-preserving.

Informally, the ancillary space $\Hil_Z$ serves the same role as the Kraus index. If $\Phi$ is a channel with Kraus operators $\{K_i : i \in \{0, \dots, k - 1\}\}$, then $A = B = \sum_i K_i \otimes \ket{i}$, such that summation over $i$ is performed by the partial trace in \eqref{eq:stinespring-action}.

The Stinespring representation is useful for evaluating norms and for building system/environment models of a quantum channel.

Further Resources

Preamble

Features

We enable a few features such that this notebook runs in both Python 2 and 3.


In [1]:
from __future__ import division, print_function

Imports


In [2]:
import numpy as np
import qutip as qt

from qutip.ipynbtools import version_table

Plotting Support


In [3]:
%matplotlib inline

Settings


In [4]:
qt.settings.colorblind_safe = True

Dilations of CPTP Maps

We'll start by considering the case in which $\Phi(X) = U X U^\dagger$ for a random unitary operator $U \in \Uni(\Hil_X, \Hil_X)$.


In [5]:
phi = qt.to_super(qt.rand_unitary_haar(4))

In [6]:
qt.visualization.hinton(phi)


Out[6]:
(<matplotlib.figure.Figure at 0x10e019c18>,
 <matplotlib.axes._subplots.AxesSubplot at 0x1167fd630>)

By construction, $\Phi$ has one Kraus operator--- the random unitary itself. Because of numerical precision in finding eigenvalues of the Choi matrix, a decomposition will result in many operators with very small norms.


In [7]:
for K in qt.to_kraus(phi):
    print(K.norm())


1.44965724751e-08
4.0
2.17301265733e-08
1.52222408748e-08
1.5017315825e-08
1.76642230882e-08
1.37895174864e-08
1.40385925527e-08
1.07348372039e-08
9.50757741411e-09
6.71512203218e-09
1.14139195174e-08
1.01210695111e-08
6.00690608844e-09
3.38230626595e-09
4.01926307771e-10

Because $\Phi$ only has one Kraus operator, the Stinespring dilation will use a 1-dimensional ancillary space.


In [8]:
A, B = qt.to_stinespring(phi)
print(A.dims)


[[4, 1], [4]]

If we ask for a random channel of Choi rank 2 (that is, with 2 Kraus operators), then the corresponding space will be two-dimensional instead.


In [9]:
A, B = qt.to_stinespring(qt.rand_super_bcsz(4, rank=2))
print(A.dims)


[[4, 2], [4]]

Moreover, since we have demanded that $A$ and $B$ are the Stinespring pair for a channel, we can easily verify that $A = B$.


In [10]:
A - B


Out[10]:
Quantum object: dims = [[4, 2], [4]], shape = (8, 4), type = other\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\\end{array}\right)\end{equation*}

Differences of Random Channels

Next, let's consider two random qubit channels from the BCSZ distribution.


In [11]:
phi1 = qt.rand_super_bcsz(2)
phi2 = qt.rand_super_bcsz(2)

Looking at the Pauli-basis Hinton diagram for the difference between these channels gives some insight into how this map behaves.


In [12]:
qt.visualization.hinton(phi1 - phi2)


Out[12]:
(<matplotlib.figure.Figure at 0x10e03b320>,
 <matplotlib.axes._subplots.AxesSubplot at 0x119ffa828>)

In particular, note that $\Delta\Phi := \Phi_1 - \Phi_2$ annilates the traceful parts of its input, leaving an operator with negative eigenvalues.


In [13]:
d_phi = qt.to_super(phi1 - phi2)
rho_in = qt.rand_dm(2)
rho = qt.vector_to_operator(d_phi * qt.operator_to_vector(rho_in))
print("Trace: {}\nEigenvalues:{}".format(rho.tr(), rho.eigenenergies()))
rho


Trace: -3.885780586188048e-16
Eigenvalues:[-0.39781299  0.39781299]
Out[13]:
Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}-0.334 & (-0.211+0.051j)\\(-0.211-0.051j) & 0.334\\\end{array}\right)\end{equation*}

Because $\Delta \Phi$ is not completely positive, we cannot obtain a Kraus decomposition where the left and right operators are the same:


In [14]:
Ks = qt.to_kraus(d_phi)

In [15]:
rho - sum(K * rho_in * K.dag() for K in Ks)


Out[15]:
Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}-0.993 & (-0.141+0.001j)\\(-0.141-0.001j) & -0.327\\\end{array}\right)\end{equation*}

On the other hand, if we allow the left and right Kraus operators to be different, or consider a Stinespring dilation, we get the right answer:


In [16]:
A, B = qt.to_stinespring(d_phi)
rho - (A * rho_in * B.dag()).ptrace((0,))


Out[16]:
Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0\\0.0 & 0.0\\\end{array}\right)\end{equation*}

Constructing System-Environment Models

For a channel $\Phi \in \Chan(\Hil_X, \Hil_X)$, once we have found the Stinespring dilation $A$, we then represent $\Phi$ as a unitary on the full space, $U \in \Uni(\Hil_X \otimes \Hil_Z)$, interpreting $\Hil_Z$ as a preparation of an environment. Concretely, we want to find $U$ such that for a preparation $\ket{\psi}$ of the environment, \begin{equation} \Phi(X) = \Tr_Z(U \rho \otimes \ket{\psi}\bra{\psi} U^\dagger). \end{equation} By convention, we'll choose a basis for $\Hil_Z$ such that $\ket{\psi} = \ket{0}$. Then, if we let $V = A \otimes \bra{0}$, \begin{equation} \Tr_Z(V \rho \otimes \ket{0}\bra{0} V^\dagger) = \Tr_Z(A \rho \otimes \bra{0}\ket{0}\bra{0}\ket{0} A) = \Tr_Z(A \rho A^\dagger) \end{equation} as desired.


In [17]:
A, B = qt.to_stinespring(qt.rand_super_bcsz(3))
V = qt.tensor(A, qt.basis(A.dims[0][-1], 0).dag())
# This adds an annoying left index, so let's drop it now.
del V.dims[0][-1]
V


Out[17]:
Quantum object: dims = [[3, 9], [3, 9]], shape = (27, 27), type = oper, isherm = False\begin{equation*}\left(\begin{array}{*{11}c}-0.135 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\-0.208 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.171 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\-0.266 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.233 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\\vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots\\(-0.057-0.048j) & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\(-0.098-0.087j) & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\(0.135+0.086j) & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\(0.035+0.025j) & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\(-0.060+0.028j) & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\\end{array}\right)\end{equation*}

In [18]:
rho = qt.rand_dm_ginibre(3)
(A * rho * A.dag()).ptrace((0,)) - (V * qt.tensor(rho, qt.ket2dm(qt.basis(A.dims[0][-1], 0))) * V.dag()).ptrace((0,))


Out[18]:
Quantum object: dims = [[3], [3]], shape = (3, 3), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0\\\end{array}\right)\end{equation*}

It thus remains to extend $V$ to act non-trivially on the full space. This is done by noting that the singular values of $V$ are each one or zero, convienently partitioning $\Hil_X \otimes \Hil_Z$ into the null space of $V$ and its complement.


In [19]:
Vu, Vs, Vv = np.linalg.svd(V.data.todense())
U = qt.Qobj(Vu, dims=V.dims) * qt.Qobj(Vv, dims=V.dims)
U * U.dag()


Out[19]:
Quantum object: dims = [[3, 9], [3, 9]], shape = (27, 27), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}1.000 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 1.000 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 1.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 1.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 1.000 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\\vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 1.000 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 1.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 1.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 1.000 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 1.000\\\end{array}\right)\end{equation*}

We finish by verifying that preparing an environment state, evolving under the system/environment unitary $U$, then partial tracing over the environment produces the same map as the original channel.


In [20]:
(A * rho * A.dag()).ptrace((0,)) - (U * qt.tensor(rho, qt.ket2dm(qt.basis(A.dims[0][-1], 0))) * U.dag()).ptrace((0,))


Out[20]:
Quantum object: dims = [[3], [3]], shape = (3, 3), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0\\\end{array}\right)\end{equation*}

Epilouge


In [21]:
version_table()


Out[21]:
SoftwareVersion
QuTiP4.3.0.dev0+6e5b1d43
Numpy1.13.1
SciPy0.19.1
matplotlib2.0.2
Cython0.25.2
Number of CPUs2
BLAS InfoINTEL MKL
IPython6.1.0
Python3.6.2 |Anaconda custom (x86_64)| (default, Jul 20 2017, 13:14:59) [GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)]
OSposix [darwin]
Thu Jul 20 22:58:26 2017 MDT