Skip to content

Upcoming version: density matrix simulation

Jump to bottom
Pablo Le Hénaff edited this page Apr 28, 2023 · 18 revisions

This page is a work in progress about a not-yet-released version of QX-simulator.

The new simulator has a lot of new features, and the API has changed. The overall direction is the handling of mixed states in a efficient manner, so that having to iterate a simulation is useless and in fact this feature has been removed.

The new output of the simulator consists in different things:

  • the (sparse) state vector, in case the simulation ends on a pure state; this is returned as a Python dictionary;
  • the density matrix, in case the simulation ends on a mixed state; this is returned as a numpy array, and is disabled when using more than 64 qubits because of memory constraints... In the future a sparse version could be available;
  • the probabilities of obtaining each possible measurement register, as a Python dictionary mapping bitstrings to floating-point probabilities.

The input of the simulator is still a cQasm file or string, however, you can also specify your own quantum operations through (parametrized) Kraus operators. More details and examples below.

To raise some ambiguity and have an easier API for specifying quantum (and classical) operations, classical bits are no longer implicit when measuring a qubit. That is, in the new version, the usual measure operation of a qubit is written, for instance, measure q[0], b[0] to store the output in b[0].

One reason for this, is that with the full generality offered by the new API, you can specify measurements that don't affect the measurement register. This would conflict with the previous measure operation which always sets the bit result.

In the future or just for compatibility, we could add a notation or even add syntactic sugar somehow to support the old way. This will need to be discussed with the users.

Example output of a simple mixed state

>>> import qxelarator
>>> qxelarator.execute_string("""
version 1.0
qubits 1

h q[0]
measure q[0], b[0]
""")

Density matrix:
[[0.5+0.j 0. +0.j]
 [0. +0.j 0.5+0.j]]

Measurement register probabilities:
{'1': 0.4999999999999999, '0': 0.4999999999999999}

Example output of a simple pure state

>>> qxelarator.execute_string("""
version 1.0
qubits 1

h q[0]
""")

State:
{'0': (0.7071067811865475+0j), '1': (0.7071067811865475+0j)}

Measurement register probabilities:
{'0': 0.9999999999999998}

Output

In case the simulation fails for whatever reason, a SimulationError object is returned with a string content, like before.

>>> output = qxelarator.execute_string("""bl;""")
>>> output
Quantum simulation error: Error: <unknown>:1:1: syntax error
>>> output.message
'Error: <unknown>:1:1: syntax error'

Error messages have been improved a lot:

>>> qxelarator.execute_string("""version 1.0;qubits 1;not q[0]""").message
'Requested operation 'not' with signature (Qubit), but the registered operation with that name has signature (ClassicalBit)'

If the simulation succeeds, there are three fields: state, densityMatrix and results. For a pure state, SimulationResult.densityMatrix is None. However, the inference of purity of the state is not perfect: you might get a density matrix even though the state is pure.

>>> output = qxelarator.execute_string("""version 1.0;qubits 1;x q[0]""")
>>> output.state
{'1': (1+0j)}
>>> output.densityMatrix
>>> output.results
{'0': 1.0}

For a mixed state:

>>> qxelarator.execute_string("""version 1.0;qubits 2;rx q[0], 1.23;measure q[0], b[0]""").densityMatrix
array([[0.66711886+0.j, 0.        +0.j, 0.        +0.j, 0.        +0.j],
       [0.        +0.j, 0.33288114+0.j, 0.        +0.j, 0.        +0.j],
       [0.        +0.j, 0.        +0.j, 0.        +0.j, 0.        +0.j],
       [0.        +0.j, 0.        +0.j, 0.        +0.j, 0.        +0.j]])

Writing quantum operations

With the new version of QX-simulator, every single instruction in the circuit can be specified from Python (or C++). This includes quantum channels with parameters, which need to be passed as lambdas. This is all done with the operations parameter to qxelarator.execute_string. A valid operations parameter is a Python dict whose keys are strings (the names of the operations) and whose values are tuples consisting of a signature (a string) and a Python callable returning a list of Kraus operators (numpy arrays).

For example, the amplitude damping channel as described by Nielsen & Chuang can be defined as such:

import numpy as np
import math

def amplitude_damping_channel(gamma):
  k0 = np.array([[1, 0                   ],
                 [0, math.sqrt(1 - gamma)]])
  k1 = np.array([[0, math.sqrt(gamma)],
                 [0, 0               ]])
  return [k0, k1]

amplitude_damping_channel_signature = "qd" # first argument is a qubit, second is a double/floating point number

def x():
  return [np.array([[0, 1], [1, 0]])] # Needs to be an array even though there is a single Kraus operator.

operations = {"ampl_damp": (amplitude_damping_channel_signature, amplitude_damping_channel), "x": ("q", x)}

simulationOutput = qxelarator.execute_string("version 1.0; qubits 1; x q[0]; ampl_damp q[0], 0.8", operations=operations)
print(simulationOutput.densityMatrix)

which should print:

Density matrix:
[[0.8+0.j 0. +0.j]
 [0. +0.j 0.2+0.j]]

Measurement register probabilities:
{'0': 0.9999999999999999}

Note that unfortunately, at the moment, passing operations means having to redefine ALL operations without being able to use the default ones (which is why the Pauli X matrix had to be defined above). This is annoying and can be fixed soon.

Signature of quantum operations

As already seen above, quantum operations need to be specified with a "signature" which defines the types of arguments they can be called with in the input cQasm file or string, from left to right. The types correspond to characters in the strings and are as follows:

  • "q" is a qubit, e.g. q[2],
  • "b" is a classical bit (from the measurement register), e.g. b[4],
  • "d" is a floating-point value, e.g. 3.14159
  • "i" is a signed integer, e.g. 3 or -8.

So, defining an operation "op" with signature "iqd" means that op 3, q[2], 0.23 is a valid call, however op q[2], 2.3 is not. The order, types and number of the operands must match the signature exactly. Otherwise, a SimulationError object is returned.

Static vs dynamic operands

Static operands are integers "i" and floating-point numbers "d". Their value is directly present in the input circuit and they are passed in the proper ordering to the Kraus-operators generation function specified in the operations parameter. Other operand types (qubit and classical bits) are dynamic operands. Their state is only known during execution of the circuit.

So, for example, an operation with signature "qqd" has two dynamic qubit operands and a single float static operand and therefore its generation function need to be callable with a single float parameter. An example of such an operation could be a controlled rotation:

import typing
def crx(theta: float):
  u = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, math.cos(theta/2), -1j*math.sin(theta / 2)], [0, 0, -1j*math.sin(theta/2), math.cos(theta/2)]])
  return [u]

operations = {"x": ("q", x), "crx": ("qqd", crx)}

print(qxelarator.execute_string("version 1.0;qubits 2; x q[0]; crx q[0], q[1], 2.5421", operations=operations))

The result is:

State:
{'01': (0.29527785388641337+0j), '11': -0.9554114239448018j}

Measurement register probabilities:
{'00': 0.9999999999999999}

Classical bits as operands and custom measurement operators

In the new simulator, classical bits are no longer implicit in the input cQasm. This is to remove ambiguity and allow specifying operations completely, including their classical effects. Classical bits as dynamic operands will behave like qubits, that is, they are operated on by matrices. They need to be thought of as ancilla qubits storing the results of operations. When retrieving the final density matrix, they are traced out. This approach is required to account for intermediate measurements inside the circuit, in the density matrix simulation framework, without having to use an increasing number of ancilla qubits which would be required by the deferred measurement technique.

For example, if a custom operation takes as operands a qubit and a classical bit, its Kraus matrices need to be of size 4 * 4. This allows to specify precisely which output 0 or 1 corresponds to which Kraus operation, or you could even define operations which always set a classical bit to 1.

In this regard, the classical not gate is defined using the same matrix as the quantum Pauli-X gate.

To see all the power of this approach, let's define measurement operators for the bell basis, as described on Wikipedia: $$\left| \Phi^+ \right&gt; = \frac{1}{\sqrt{2}} \left( \left| 00 \right&gt; + \left| 11 \right&gt; \right)$$ $$\left| \Phi^- \right&gt; = \frac{1}{\sqrt{2}} \left( \left| 00 \right&gt; - \left| 11 \right&gt; \right)$$ $$\left| \Psi^+ \right&gt; = \frac{1}{\sqrt{2}} \left( \left| 01 \right&gt; + \left| 10 \right&gt; \right)$$ $$\left| \Psi^- \right&gt; = \frac{1}{\sqrt{2}} \left( \left| 01 \right&gt; - \left| 10 \right&gt; \right)$$

Assume we want to write the result of the measurement in two classical bits $b_1b_0$ such that measuring $\left| \Phi^+ \right&gt;$ (respectively $\left| \Phi^- \right&gt;$, $\left| \Psi^+ \right&gt;$ and $\left| \Psi^- \right&gt;$) gives $00$ (respectively $01$, $10$ and $11$). This can be written with the QX-simulator API as follows:

phiPlus = 1/math.sqrt(2) * np.array([1, 0, 0, 1])
phiMinus = 1/math.sqrt(2) * np.array([1, 0, 0, -1])
psiPlus = 1/math.sqrt(2) * np.array([0, 1, 1, 0])
psiMinus = 1/math.sqrt(2) * np.array([0, 1, -1, 0])

def oneOne(i, j):
    res = np.zeros((4, 4))
    res[i, j] = 1
    return res

k0 = [np.kron(np.outer(phiPlus, phiPlus.conj().T), oneOne(0, j)) for j in range(4)]
k1 = [np.kron(np.outer(phiMinus, phiMinus.conj().T), oneOne(1, j)) for j in range(4)]
k2 = [np.kron(np.outer(psiPlus, psiPlus.conj().T), oneOne(2, j)) for j in range(4)]
k3 = [np.kron(np.outer(psiMinus, psiMinus.conj().T), oneOne(3, j)) for j in range(4)]

measure_bell = ("qqbb", lambda: [i for l in [k0, k1, k2, k3] for i in l])

Note that this quantum operation operates on 2 qubits and 2 classical bits, and makes use of 16 Kraus operators.

Measuring quantum states with this new operation works as expected:

h q[0] 
cnot q[0],q[1]

measure_bell q[0], q[1], b[1], b[0] # 00

x q[0]
h q[0] 
cnot q[0],q[1]

measure_bell q[0], q[1], b[1], b[0] # 01

Kraus operators specification

As described in textbooks, a quantum operation $\mathcal{E}$ is defined by its associated Kraus operators $E_1, \dots, E_n$ such that for a density operator $\rho$: $$\mathcal{E}(\rho) = \sum_i E_i \rho E_i^\dagger$$ and is only valid if those satisfy the following inequality: $$\sum_i E_i^\dagger E_i \le I\ \text{,}$$ which means that the maximum of the eigenvalues of the left member (spectral radius) is not greater than one. This is called the trace-preserving, or trace-non-increasing property.

This is checked by the simulator using a power-iteration eigenvalues algorithm; if we try to specify an operation that doesn't satisfy this requirement, a SimulationError is output:

>>> qxelarator.execute_string("version 1.0;qubits 1;myop q[0]", operations = {"myop": ("q", lambda: [np.array([[1, 1], [0, 1]])])})
Quantum simulation error: Kraus operators for operation myop are not non-trace-increasing

Note that in some settings it is the case that quantum operations are always considered trace-preserving, that is, the inequality in the above expression is in fact an equality. In the future, we could add an option to enable trace-decreasing quantum operations, and disable it by default.

When trace-decreasing operations are used, the state vectors and density matrices obtained by the simulation don't have the usual properties of squared norms of amplitudes summing to 1 or trace equal to 1, respectively. Currently in the debug version of QX-simulator (therefore not the released Python package) there are assertions which check for those properties, which means that those will crash the program if such trace-decreasing operations are used. The desired behavior is to be discussed with users.

Performance

The new simulator is designed to offer zero-cost density matrix simulation when all operations are unitaries and computational-basis measurements. That is, when a circuit consists of only unitaries and measure_z, the runtime of a simulation should be in the same order of magnitude compared to the previous state-vector simulator that was non-deterministically branching as mixed-states were created.

This means that the previous state-vector simulator is obsolete and completely superseeded by this new release.

Clone this wiki locally