Density Matrix and Mixed State Evolution

Overview

TensorCircuit provides two methods of simulating noisy, mixed state quantum evolution. Full density matrix simulation of \(n\) qubits is provided by using tc.DMCircuit(n), and then adding quantum operations – both unitary gates as well as general quantum operations specified by Kraus operators – to the circuit. Relative to pure state simulation of \(n\) qubits via tc.Circuit, full density matrix simulation is twice as memory-intensive, and thus the maximum system size simulatable will be half of what can be simulated in the pure state case. A less memory intensive option is to use the standard tc.Circuit(n) object and stochastically simulate open system evolution via Monte Carlo trajectory methods.

Setup

import numpy as np
import tensorcircuit as tc

K = tc.set_backend("tensorflow")

Density Matrix Simulation with tc.DMCircuit

We illustrate this method below, by considering a simple circuit on a single qubit, which takes as input the mixed state corresponding to a probabilistic mixture of the \(\vert{0}\rangle\) state and the maximally mixed state \(\rho(\alpha) = \alpha\vert 0\rangle \langle 0\vert + (1-\alpha)I/2.\)

This state is then passed through a circuit that applies an \(X\) gate, followed by a quantum operation corresponding to an amplitude damping channel \(\mathcal{E}_\gamma\) with parameter \(\gamma\). This has Kraus operators \(K_0 = \begin{pmatrix} 1 & 0 \\ 0 & \sqrt{1-\gamma} \end{pmatrix}, \quad K_1 = \begin{pmatrix} 0 & \sqrt{\gamma} \\ 0 & 0 \end{pmatrix}\) This circuit thus causes the evolution \(\rho(\alpha) \xrightarrow[]{X} X\rho(\alpha)X\xrightarrow[]{\mathcal{E}_\gamma}\sum_{i=0}^1 K_i X\rho(\alpha)X K_i^\dagger\)

To simulate this in TensorCircuit, we first create a tc.DMCircuit (density matrix circuit) object and set the input state using the dminputs optional argument (note that if a pure state input is provided to tc.DMCircuit, this should be done via the inputs optional argument).

\(\rho(\alpha)\) has matrix form \(\rho(\alpha) = \begin{pmatrix} \frac{1+\alpha}{2} & \\ & \frac{1-\alpha}{2} \end{pmatrix},\) and thus (taking \(\alpha=0.6\)) we initialize the density matrix circuit as follows.

To implement a general quantum operation such as the amplitude damping channel, we use general_kraus, supplied with the corresponding list of Kraus operators.

def rho(alpha):
    return np.array([[(1 + alpha) / 2, 0], [0, (1 - alpha) / 2]])


input_state = rho(0.6)
dmc = tc.DMCircuit(1, dminputs=input_state)

dmc.x(0)


def amp_damp_kraus(gamma):
    K0 = np.array([[1, 0], [0, np.sqrt(1 - gamma)]])
    K1 = np.array([[0, np.sqrt(gamma)], [0, 0]])
    return K0, K1


K0, K1 = amp_damp_kraus(0.3)
dmc.general_kraus([K0, K1], 0)  # apply channel with Kraus operators [K0,K1] to qubit 0

# get the output density matrix
dmc.state()

<tf.Tensor: shape=(2, 2), dtype=complex64, numpy=
array([[0.44+0.j, 0.  +0.j],
       [0.  +0.j, 0.56+0.j]], dtype=complex64)>

# evaluate the expectation as a circuit object
print(dmc.expectation_ps(z=[0]), dmc.measure(0))

tf.Tensor((-0.11999999+0j), shape=(), dtype=complex64) (<tf.Tensor: shape=(1,), dtype=float32, numpy=array([1.], dtype=float32)>, -1.0)

In the above example, we input the Kraus operators for the amplitude damping channel manually, in order to illustrate the general approach to implementing quantum channels. In fact, TensorCircuit includes built-in methods for returning the Kraus operators for a number of common channels, including the amplitude damping, depolarizing, phase damping, and reset channels.

# a set of built-in quantum channels

for k in dir(tc.channels):
    if k.endswith("channel"):
        print(k)

amplitudedampingchannel
depolarizingchannel
phasedampingchannel
resetchannel

dmc = tc.DMCircuit(2)
dmc.h(0)
gamma = 0.2
K0, K1 = tc.channels.phasedampingchannel(gamma)
dmc.general_kraus([K0, K1], 0)
dmc.state()

<tf.Tensor: shape=(4, 4), dtype=complex64, numpy=
array([[0.49999997+0.j, 0.        +0.j, 0.4472136 +0.j, 0.        +0.j],
       [0.        +0.j, 0.        +0.j, 0.        +0.j, 0.        +0.j],
       [0.4472136 +0.j, 0.        +0.j, 0.49999994+0.j, 0.        +0.j],
       [0.        +0.j, 0.        +0.j, 0.        +0.j, 0.        +0.j]],
      dtype=complex64)>

# or we can directly use the following API for shorthand

dmc = tc.DMCircuit(2)
dmc.h(0)
gamma = 0.2
dmc.phasedamping(0, gamma=0.2)
dmc.state()

<tf.Tensor: shape=(4, 4), dtype=complex64, numpy=
array([[0.49999997+0.j, 0.        +0.j, 0.4472136 +0.j, 0.        +0.j],
       [0.        +0.j, 0.        +0.j, 0.        +0.j, 0.        +0.j],
       [0.4472136 +0.j, 0.        +0.j, 0.49999994+0.j, 0.        +0.j],
       [0.        +0.j, 0.        +0.j, 0.        +0.j, 0.        +0.j]],
      dtype=complex64)>

AD and JIT Compatibility

tc.DMCircuit, like tc.Circuit is also compatible with ML paradigms such as AD, jit, and vmap. See the example below.

n = 3
nbatch = 2


def loss(params, noisep):
    c = tc.DMCircuit(n)
    for i in range(n):
        c.rx(i, theta=params[i])
    for i in range(n):
        c.depolarizing(i, px=noisep, py=noisep, pz=noisep)
    return K.real(K.sum([c.expectation_ps(z=[i]) for i in range(n)]))


loss_vvg = K.jit(
    K.vectorized_value_and_grad(loss, argnums=(0, 1), vectorized_argnums=(0))
)

vs, (gparams, gnoisep) = loss_vvg(0.1 * K.ones([nbatch, n]), 0.1 * K.ones([]))

vs.shape, gparams.shape, gnoisep.shape

(TensorShape([2]), TensorShape([2, 3]), TensorShape([]))

Note how the noise parameter can also be differentiated and jitted!

Monte Carlo Noise Simulation with tc.Circuit

For pure state inputs, Monte Carlo methods can be used to sample noisy quantum evolution using tc.Circuit instead of tc.DMCircuit where the mixed state is effectively simulated with an ensemble of pure states.

As for density matrix simulation, quantum channels \(\mathcal{E}\) can be added to a circuit object by providing a list of their associated Kraus operators \(\{K_i\}\). The API is the same as for the full density matrix simulation.

input_state = np.array([1, 1] / np.sqrt(2))
c = tc.Circuit(1, inputs=input_state)
c.general_kraus(tc.channels.phasedampingchannel(0.5), 0)
c.state()

<tf.Tensor: shape=(2,), dtype=complex64, numpy=array([0.+0.j, 1.+0.j], dtype=complex64)>

In this framework though, the output of a channel acting on \(\vert{\psi}\rangle\) , i.e. \(\mathcal{E} ( \vert{\psi}\rangle\langle{\psi}\vert) = \sum_i K_i \vert{\psi}\rangle\langle{\psi}\vert K_i^ \dagger\) is viewed as an ensemble of states \(\frac{K_i\vert{\psi}\rangle}{\sqrt{\langle{\psi}\vert K_i^\dagger K_i \vert{\psi}\rangle}}\) that each occur with probability \(p_i = \langle{\psi}\vert K_i^\dagger K_i \vert{\psi}\rangle\). Thus, the code above stochastically produces the output of a single qubit initialized in state \(\vert{\psi}\rangle=\frac{\vert{0}\rangle+\vert{1}\rangle}{\sqrt{2}}\) being passed through a phase damping channel with parameter \(\gamma=0.5\).

The Monte Carlo simulation of channels where the Kraus operators are all unitary matrices (up to a constant factor) can be handled with additional efficiency by using unitary_kraus instead of general_kraus.

px, py, pz = 0.1, 0.2, 0.3
c.unitary_kraus(tc.channels.depolarizingchannel(px, py, pz), 0)

<tf.Tensor: shape=(), dtype=int32, numpy=3>

Note the int tensor returned above indicates in this trajectory, which operator is applied on the circuit.

Externalizing the Randomness

The general_kraus and unitary\_kraus examples above both handle randomness generation from inside the respective methods. That is, when the list \([K_0, K_1, \ldots, K_{m-1}]\) of Kraus operators is supplied to general_kraus or unitary_kraus, the method partitions the interval \([0,1]\) into \(m\) contiguous intervals \([0,1] = I_0 \cup I_1 \cup \ldots I_{m-1}\) where the length of \(I_i\) is equal to the relative probability of obtaining outcome \(i\). Then a uniformly random variable \(x\) in \([0,1]\) is generated from within the method, and outcome \(i\) selected based on which interval \(x\) lies in.

In TensorCircuit, we have full backend agnostic infrastructure for random number generation and management. However, the interplay between jit, random number, and backend switch is often subtle if we rely on the random number generation inside these methods. See advance.html#randoms-jit-backend-agnostic-and-their-interplay for details.

In some situations, it may be preferable to first generate the random variable from outside the method, and then pass the value generated into general_kraus or unitary_kraus. This can be done via the optional status argument:

px, py, pz = 0.1, 0.2, 0.3
x = 0.5
print(c.unitary_kraus(tc.channels.depolarizingchannel(px, py, pz), 0, status=x))
x = 0.8
print(c.unitary_kraus(tc.channels.depolarizingchannel(px, py, pz), 0, status=x))

tf.Tensor(1, shape=(), dtype=int32)
tf.Tensor(3, shape=(), dtype=int32)

This is useful, for instance, when one wishes to use vmap to batch compute multiple runs of a Monte Carlo simulation. This is illustrated in the example below, where vmap is used to compute 10 runs of the simulation in parallel.

def f(x):
    c = tc.Circuit(1)
    c.h(0)
    c.unitary_kraus(tc.channels.depolarizingchannel(0.1, 0.2, 0.3), 0, status=x)
    return c.expectation_ps(x=[0])


f_vmap = K.vmap(f, vectorized_argnums=0)
X = K.implicit_randn(10)
f_vmap(X)

<tf.Tensor: shape=(10,), dtype=complex64, numpy=
array([ 0.99999994+0.j,  0.99999994+0.j,  0.99999994+0.j, -0.99999994+0.j,
        0.99999994+0.j,  0.99999994+0.j,  0.99999994+0.j,  0.99999994+0.j,
       -0.99999994+0.j,  0.99999994+0.j], dtype=complex64)>