Tutorial

QCircuits is a library for simulating quantum circuits. Its primary classes are State, representing the (quantum) state of the computer, a unit vector in a complex vector space, and Operator, representing quantum gates, i.e. unitary operators on those vector spaces. This tutorial describes how to use QCircuits to prepare operators and states, apply operators to states, measure states, etc., in order to implement quantum algorithms.

First, install QCircuits with pip:

pip install qcircuits

In the following, it is assumed that QCircuits has been imported with:

>>> import qcircuits as qc

States

States in QCircuits represent the state of a \(d\)-qubit quantum computer. This section introduces some of the ways in which states can be prepared in QCircuits.

States in an all-zero or all-one computational basis state may be prepared with the zeros() and ones() functions. E.g., to prepare the state \(|\phi⟩ = |000⟩\):

>>> phi = qc.zeros(3)  # the state |000⟩

And to prepare the state \(|\phi⟩ = |11⟩\):

>>> phi = qc.ones(2)  # the state |11⟩

The bitstring() function allows one to prepare a state in an arbitrary computational basis state. E.g., to prepare the state \(|\phi⟩ = |01001⟩\):

>>> phi = qc.bitstring(0, 1, 0, 0, 1)  # the state |01001⟩

A single qubit may be prepared with the qubit() function. A qubit \(|\phi⟩ = \alpha |0⟩ + \beta |1⟩\) may be prepared by providing \(\alpha\) and \(\beta\) such that \(\lvert\alpha\rvert^2 + \lvert\beta\rvert^2 = 1\), e.g.,

>>> phi = qc.qubit(alpha=0.5, beta=np.sqrt(3)/2)  # the state 0.5|0⟩ + 0.866|1⟩

Alternatively, a qubit can be prepared in the state \(e^{i \omega} \big( \cos(\theta/2) |0⟩ + e^{i \phi} \sin(\theta/2) |1⟩ \big)\), where \(\omega\) is the global phase of the state, using the same function, e.g.,

>>> phi = qc.qubit(theta=math.pi, phi=0, global_phase=math.pi/2)

The four Bell states can be prepared using the bell_state() function. This takes two binary arguments x and y, and produces the Bell state \(|\beta_{xy}⟩ = \big( |0, y⟩ + (-1)^x |1, 1-y⟩ \big)/\sqrt{2}\). E.g., the Bell state \(|\beta_{00}⟩ = \frac{|00⟩ + |11⟩}{\sqrt{2}}\) can be prepared:

>>> beta = qc.bell_state(0, 0)  # the state (|0⟩ + |1⟩)/1.414

The positive_superposition() function may be used to prepare a d-qubit state in the positive equal superposition of the computational states. E.g., to construct the 2-qubit state \(|\phi⟩ = \big(|00⟩ + |01⟩ + |10⟩ + |11⟩ \big) / 2\):

>>> phi = qc.positive_superposition(d=2)

States can also be prepared by applying operators to states or taking the tensor product of states, each of which is described in later sections.

How are States Represented?

Internally, QCircuits encodes a d-qubit state with an array of shape (2, 2, …, 2), with d axes in total, representing a tensor with d contravariant indices. E.g., a 3-qubit state is represented by an array of shape (2, 2, 2), and indexing into this array with indices i, j, k gets the probability amplitude for the computational basis vector \(|ijk⟩\). The shape and the rank (number of axes) can be accessed with the State.shape and State.rank properties.

>>> phi = qc.zeros(3)
>>> print(phi)
3-qubit state. Tensor:
[[[1.+0.j 0.+0.j]
  [0.+0.j 0.+0.j]]

 [[0.+0.j 0.+0.j]
  [0.+0.j 0.+0.j]]]

>>> print(phi.shape)
(2, 2, 2)

>>> print(phi.rank)
3

>>> print(phi[0, 0, 0])  # the probability amplitude of our state for |000⟩
(1+0j)

>>> print(phi[0, 0, 1])  # the probability amplitude of our state for |001⟩
0j

A d-qubit state can be constructed by providing this array. E.g., a 3-qubit state can be constructed by providing a (2, 2, 2) shape array:

>>> phi = qc.State([[[1., 0.],  # the state |000⟩
...                  [0., 0.]],
...                 [[0., 0.],
...                  [0., 0.]]])

An alternative and common representation of a d-qubit state is as a column vector of length \(2^d\). This column-vector representation can be obtained with the State.to_column_vector() method:

>>> phi = qc.bitstring(0, 1, 0)
>>> phi.to_column_vector()
array([0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j])

States can be constructed from the column vector representation using the State.from_column_vector() static method:

>>> phi = qc.State.from_column_vector(   # the state |010⟩
...     [0., 0., 1., 0., 0., 0., 0., 0.]
... )

Note that while QCircuits allows the user to extract the tensor data from a state, i.e., observe a state in full, rather than just taking a measurement, in true quantum computation this is not possible.

Operators

Operators in QCircuits represent quantum ‘gates’ for \(d\)-qubit quantum computers, i.e., unitary linear operators on a \(2^d\) dimensional complex vector space.

This section describes some of the built-in operators that may be used.

The PauliX(), PauliY(), and PauliZ() functions return instances of the common X, Y, and Z gates. Here they are shown acting on some computational basis vectors as expected:

>>> X = qc.PauliX()  # The X gate, or NOT gate
>>> phi = qc.zeros(1)
>>> result = X(phi)  # apply the X gate to the state |0⟩
>>> print(result)    # the result is the state |1⟩
1-qubit state. Tensor:
[0.+0.j 1.+0.j]
>>> Y = qc.PauliY()
>>> phi = qc.zeros(1)
>>> result = Y(phi) # apply the Y gate to the state |0⟩
>>> print(result)   # the result is the state i|1⟩
1-qubit state. Tensor:
[0.+0.j 0.+1.j]
>>> Z = qc.PauliZ()
>>> result = Z(qc.zeros(1))  # apply the Z gate to the state |0⟩
>>> print(result)  # the result is the state |0⟩
1-qubit state. Tensor:
[1.+0.j 0.+0.j]

>>> result = Z(qc.ones(1))  # apply the Z gate to the state |1⟩
>>> print(result)  # the result is the state -|1⟩
1-qubit state. Tensor:
[-0.-0.j -1.-0.j]

We have seen in the above examples that operators are applied to states by function application, i.e., U(v), where U is an operator and v a state. Operator application will be described in more detail later in the tutorial.

An instance of the Hadamard gate can be obtained with the Hadamard() function. Here we see an example of applying the Hadamard operator to the state \(|0⟩\):

>>> H = qc.Hadamard()
>>> phi = qc.zeros(1)  # the state |0⟩
>>> result = H(phi)    # the state (|0⟩ + |1⟩) / sqrt(2)
>>> print(result)
1-qubit state. Tensor:
[0.70710678+0.j 0.70710678+0.j]

The above functions, which return instances of the the X, Y, Z, and Hadamard gates, take an integer argument d. The returned d-qubit operator applies the gate in question to each qubit independently. E.g., for the X gate, the returned operator is \(X^{\otimes d}\).

>>> X = qc.PauliX(d=3)  # The 3-qubit operator applying X to each qubit
>>> phi = qc.bitstring(0, 1, 0)  # the state |010⟩
>>> result = X(phi)  # the state |101⟩
>>> print(result)
1-qubit state. Tensor:
[[[0.+0.j 0.+0.j]
  [0.+0.j 0.+0.j]]

 [[0.+0.j 1.+0.j]
  [0.+0.j 0.+0.j]]]

By default, dimensionality of the operator and the state it is applied to must match. I.e., a d-qubit operator must be applied to a d-qubit state. Ways of applying n-qubit operators to d-qubit states where n is less than d will be discussed in the sections on tensor products and operator application below.

The CNOT() function returns an instance of the CNOT operator, i.e., the 2-qubit operator that applies the X operator to the second (target) qubit if the first (control) qubit is in state \(|1⟩\).

>>> CNOT = qc.CNOT()
>>> print(CNOT(qc.bitstring(0, 0)))   # |00⟩ -> |00⟩
2-qubit state. Tensor:
[[1.+0.j 0.+0.j]
 [0.+0.j 0.+0.j]]

>>> print(CNOT(qc.bitstring(0, 1)))   # |01⟩ -> |01⟩
2-qubit state. Tensor:
[[0.+0.j 1.+0.j]
 [0.+0.j 0.+0.j]]

>>> print(CNOT(qc.bitstring(1, 0)))   # |10⟩ -> |11⟩
2-qubit state. Tensor:
[[0.+0.j 0.+0.j]
 [0.+0.j 1.+0.j]]

>>> print(CNOT(qc.bitstring(1, 1)))   # |11⟩ -> |10⟩
2-qubit state. Tensor:
[[0.+0.j 0.+0.j]
 [1.+0.j 0.+0.j]]

Often, we want to swap the role of the qubits, flipping the first qubit if the second is set, or more generally, for a d-qubit state we may want to apply the 2-qubit CNOT operator on any two qubits. How this is done is described in the section on operator application below.

The ControlledU() function takes a d-qubit operator as an argument, and returns the (d+1)-qubit controlled-U operator: if the first qubit is set, the operator U is applied to the following d qubits.

>>> phi0 = qc.bitstring(0, 0, 0)  # the 3-qubit state |000⟩
>>> phi1 = qc.bitstring(1, 0, 0)  # the 3-qubit state |100⟩
>>> H = qc.Hadamard(d=2)  # a 2-qubit operator, applying the Hadamard operator to each qubit
>>> c_H = qc.ControlledU(H)  # a 3-qubit operator, applying the Hadamard operator
...                          # to qubits 2 and 3 if qubit 1 is set
>>> print(c_H(phi0))  # the state is left unchanged
3-qubit state. Tensor:
[[[1.+0.j 0.+0.j]
  [0.+0.j 0.+0.j]]

 [[0.+0.j 0.+0.j]
  [0.+0.j 0.+0.j]]]

>>> print(c_H(phi1))  # the 2-qubit H operator is applied
1-qubit state. Tensor:
[[[0. +0.j 0. +0.j]
  [0. +0.j 0. +0.j]]

 [[0.5+0.j 0.5+0.j]
  [0.5+0.j 0.5+0.j]]]

The U_f() function takes two arguments: a function f and an integer d. The function f must be a boolean function of d-1 boolean arguments. This returns a d-qubit operator whose action is to flip the last qubit if the result of applying the boolean function to the first d-1 qubits is one. An example of its use can be found in the Deutsch-Jozsa algorithm in the examples page.

For a full list of available operators, see Operator.

How are Operators Represented?

Internally, QCircuits encodes an operator for a d-qubit system with an array of shape (2, 2, …, 2), with 2d axes in total, representing a tensor with d contravariant indices and d covariant indices. E.g., an operator for a 2-qubit system is represented by an array of shape (2, 2, 2, 2). The shape and the rank (number of axes) can be accessed with the Operator.shape and Operator.rank properties.

>>> H = qc.Hadamard()
>>> print(H)

Operator for 1-qubit state space. Tensor:
[[ 0.70710678+0.j  0.70710678+0.j]
 [ 0.70710678+0.j -0.70710678+0.j]]

>>> print(H.shape)
(2, 2)

>>> print(H.rank)
2

We use the convention that the covariant and contravariant indices alternate. The result is that indexing into the array representing operator U with indices i, j, k, … in the odd-numbered places gets the state the computational basis vector \(|ijk\ldots⟩\) is taken to by the operator.

>>> print(U[:, 0])  # Indexing in with 0 here gets the tensor representation of
...                 # the operator applied to the state |0⟩
[ 0.53114041-0.31105474j -0.69143236-0.37822758j]

>>> phi = qc.zeros(1)
>>> print(U(phi))
[ 0.53114041-0.31105474j -0.69143236-0.37822758j]

>>> print(V[:, 1, :, 0])  # Indexing in with 1, 0 here gets the tensor representation
...                       # of the operator applied to the state |10⟩
array([[-0.66947579+0.19664594j, -0.37841556-0.24010317j],
       [ 0.30464249-0.40638463j,  0.16243857-0.16716121j]])

>>> phi = qc.bitstring(1, 0)
>>> print(V(phi))
array([[-0.66947579+0.19664594j, -0.37841556-0.24010317j],
       [ 0.30464249-0.40638463j,  0.16243857-0.16716121j]])

An operator can be constructed by providing an array of the appropriate shape. E.g., the two qubit Hadamard gate \(H\otimes H\) can be constructed by providing the (2, 2, 2, 2)-shape array:

>>> H = qc.Operator([[[[ 0.5,  0.5],
...                    [ 0.5, -0.5]],
...                   [[ 0.5,  0.5],
...                    [ 0.5, -0.5]]],
...
...                  [[[ 0.5,  0.5],
...                    [ 0.5, -0.5]],
...                   [[-0.5, -0.5],
...                    [-0.5,  0.5]]]])

An alterantive and common representation of d-qubit operators is as a \(2^d \times 2^d\) matrix. This matrix representation can be accessed with the Operator.to_matrix() method. E.g., for the two-qubit Hadamard gate:

>>> print(H.to_matrix())
[[ 0.5+0.j  0.5+0.j  0.5+0.j  0.5+0.j]
 [ 0.5+0.j -0.5+0.j  0.5+0.j -0.5+0.j]
 [ 0.5+0.j  0.5+0.j -0.5+0.j -0.5+0.j]
 [ 0.5+0.j -0.5+0.j -0.5+0.j  0.5+0.j]]

Operators can be constructed from this matrix representation using the Operator.from_matrix() static method:

>>> H = qc.Operator(
...     [[ 0.5,  0.5,  0.5,  0.5],
...      [ 0.5, -0.5,  0.5, -0.5],
...      [ 0.5,  0.5, -0.5, -0.5],
...      [ 0.5, -0.5, -0.5,  0.5]]
... )

Applying Operators to States

Operators are applied to states as follows.

>>> H = qc.Hadamard()  # the single-qubit Hadamard gate
>>> x = qc.zeros(1)    # a single qubit in the |0⟩ state
>>> y = H(x)           # apply the Hadamard gate to the qubit to obtain an equal superposition state

When applying an operator to a state, a qubit_indices argument can be supplied in order to permute the order of the state for the operator application. For example, the CNOT() gate flips the second bit if the first bit is set:

>>> CNOT = qc.CNOT()        # the CNOT gate
>>> x = qc.bitstring(1, 1)  # the |11⟩ state
>>> y = CNOT(x)             # the resulting |10⟩ state

Passing in the qubit_indices argument allows us to specify the roles of the two qubits explicitly. For example, the following swaps the roles of the qubits during operator application, so that the first qubit flips if the second is set.

>>> y = CNOT(x, qubit_indices=[1, 0])  # the resulting |01⟩ state

The qubit_indices argument can also be used to apply an \(m\)-qubit operator to an \(n\)-qubit state where \(n>m\). For example, the following code applies a two qubit Hadamard operator to the first and third qubits of a six qubit state.

>>> H = qc.Hadamard(d=2)   # the two-qubit Hadamard gate
>>> x = qc.zeros(6)        # the six-qubit state |000000⟩ state
>>> y = H(x, qubit_indices=[0, 2])  # resulting state |+0+000⟩

This is useful, as it is much more efficient than expanding the \(m\)-qubit operator to a \(n\)-qubit operator by taking the tensor product with the identity and then applying the result.

Tensor Products

If a quantum system \(A\) is in state \(|\psi⟩\), and system \(B\) is in state \(|\phi⟩\), then the combined system \(A\otimes B\) is in state \(|\psi⟩ \otimes |\phi⟩\), where \(\otimes\) is the tensor product. If operator \(U\) is applied to system \(A\) and operator V applied to system \(B\), then this can be described by a single operator \(U\otimes V\) applied to system \(A\otimes B\).

This tensor product operation can be used via the infix multiplication operator * in QCircuits to produce operators and states for larger systems from operators and states for smaller systems:

>>> psi = qc.bitstring(0)   # the state |0⟩
>>> phi = qc.bitstring(11)  # the state |11⟩
>>> state = psi * phi       # the state |011⟩
>>> print(state)
3-qubit state. Tensor:
[[[0.+0.j 0.+0.j]
  [0.+0.j 1.+0.j]]

 [[0.+0.j 0.+0.j]
  [0.+0.j 0.+0.j]]]
>>> X = qc.PauliX()  # a single qubit X operator
>>> Z = qc.PauliZ(d=2)  # a 2-qubit operator applying Z independently to two qubits
>>> U = X * Z  # a 3-qubit operator, applying X to the first qubit and Z to the second and third

One use of taking the tensor product of operators is applying smaller operators to larger states. E.g., a 1-qubit operator may be applied to one of the qubits of a 2-qubit state by taking the tensor product of the operator with the identity operator.

>>> H = qc.Hadamard()  # the 1-qubit Hadamard operator
>>> I = qc.Identity()  # the 1-qubit identity operator
>>> phi = qc.bitstring(0, 0)  # the state |00⟩
>>> print((H * I)(phi))  # apply the Hadamard operator to the first qubit
...                      # resulting in state (|0⟩ + |1⟩)/sqrt(2) |0⟩
2-qubit state. Tensor:
[[0.70710678+0.j 0.        +0.j]
 [0.70710678+0.j 0.        +0.j]]

>>> print((I * H)(phi))  # apply the Hadamard operator to the second qubit
...                      # resulting in state |0⟩ (|0⟩ + |1⟩)/sqrt(2)
2-qubit state. Tensor:
[[0.70710678+0.j 0.70710678+0.j]
 [0.        +0.j 0.        +0.j]]

Since a d-qubit operator is specified with \(2^{2d}\) complex values, while a d-qubit state is specified with \(2^d\) complex values, this method is not ideal when applying very small operators to very large states. As an example, working with 30-qubit states is plausible on personal hardware, requiring 16 GB of memory. Expanding an operator to a 30-qubit operator to act on this state is not plausible, as this would require 16 exabytes (1 million TB) of memory. The previous section describes an alternative way to apply smaller operators to larger states, by specifying qubit indices.

One can use power notation to take the tensor product of an operator or state with itself \(n\) times, as in the following code that creates an \(n\)-qubit operator from a single-qubit operator:

>>> Op2 = Op**n  # creating an n-qubit operator from a single-qubit operator

Composing Operators

Operators may be applied to other operators to produce new operators, as linear operators are associative, i.e., \(A(B|\phi⟩) = (AB)|\phi⟩\). For example, suppose we start with state \(|00⟩\), and wish to apply the Hadamard gate to the first qubit, then apply the CNOT gate, resulting in a Bell state. We can either apply these operators in sequence, or we can first construct a single 2-qubit operator by composing the operators, and apply the resulting operator to the state.

>>> H = qc.Hadamard()
>>> I = qc.Identity()
>>> CNOT = qc.CNOT()
>>> phi = qc.bitstring(0, 0)

>>> result1 = CNOT((H * I)(phi))  # these result in the same state
>>> U = CNOT(H * I)  # this method produces an operator U that performs both operators
>>> result2 = U(phi)  # the result is a Bell state

Suppose we wish to apply the inverse operator, taking a Bell state to the state \(|00⟩\). First we apply the CNOT operator to the state, then we apply the 1-qubit Hadamard operator to the first qubit. Again, this can be done with operator composition:

>>> H = qc.Hadamard()
>>> I = qc.Hadamard()
>>> CNOT = qc.CNOT()
>>> phi = qc.bell_state(0, 0)
>>> result = (H * I)(CNOT)(phi)  # first H*I is applied to CNOT, composing the
...                              # operators. The result is applied to the state

In this case, though, where we are applying a smaller operator to a larger operator, we can use the same interface as when applying a smaller operator to a larger state, by specifying the qubits on which the operator acts:

>>> H = qc.Hadamard()
>>> CNOT = qc.CNOT()
>>> phi = qc.bell_state(0, 0)
>>> U = H(CNOT, qubit_indices=[0])
>>> result = U(phi)

Measurement

Measurement in the computational basis may be performed on one or more qubits of a state with the State.measure() method, which returns the result of measurement. Post-measurement, the state collapses to the computational basis state corresponding to the result of the measurement of the measured qubits. The State.measure() method has two arguments: a list of indices specifying the qubits to be measured, and a flag specifying whether the measured qubits are to be removed from the state. If no indices are supplied, every qubit is measured.

E.g., the single-qubit state \(|0⟩\) will yield the measurement 0 with certainty:

>>> qc.zeros(1).measure()
(0,)

The 3-qubit state \(|111⟩\) will yield measurement 1, 1 with certainty when the first two qubits are measured. Post-measurement we are left with either a 3- or 1-qubit state, depending on if we choose to remove the measured qubits from the state.

>>> phi = qc.ones(3)
>>> phi.measure(qubit_indices=[0, 1], remove=False)
(1,1)
>>> phi.shape
(2,2,2)

>>> phi = qc.ones(3)
>>> phi.measure(qubit_indices=[0, 1], remove=True)
(1,1)
>>> phi.shape  # if the measured qubits are removed, we are left with state |1⟩
(2,)

Measuring the state \((|0⟩ + |1⟩)/\sqrt{2}\) yields 0 or 1 with equal probability. Note that in the following code we must prepare a fresh state each time we measure to see different measurement outcomes, because the measured state is left either in the state \(|0⟩\) or in the state \(|1⟩\).

>>> H = qc.Hadamard()
>>> H(qc.zeros(1)).measure()
(0,)

>>> H(qc.zeros(1)).measure()
(1,)

>>> H(qc.zeros(1)).measure()
(0,)

Density Operators

Density operators are a useful way of representing mixed states, i.e., statistical ensembles of quantum states that can arise when there is uncertainty about the actual quantum states that are being manipulated. If a system is known to be in one of a collection of states \(|\phi_i⟩\), each with probability \(p_i\), then the density operator representation of this state is \(\rho = \sum_i p_i |\phi_i⟩⟨\phi_i|\).

A DensityOperator object can be created with the DensityOperator.from_ensemble() static method, by supplying a list of State objects and a list of probabilities, matching in length, which must sum to one:

>>> state1 = qc.bitstring(0, 1, 0)  # the state |010⟩
>>> state2 = qc.bitstring(1, 0, 0)  # the state |100⟩
>>> mixed_state = qc.DensityOperator.from_ensemble(
...     [state1, state2],
...     ps=[0.3, 0.7]
... )  # the statistical ensemble known to be in states |010⟩ and |100⟩
...    # each with probability 1/2.

The collection of states must have matching rank. I.e., one cannot form a mixture of a two qubit and a three qubit state. If the probability vector is omitted a uniform mixture is assumed.

A DensityOperator can be created from a single pure state using the State.density_operator() method:

>>> state = qc.bitstring(0, 0, 0)  # the state |000⟩
>>> rho = state.density_operator()

Operators can be applied to mixed states in the same way that they are applied to pure states:

>>> H = qc.Hadamard(d=3)  # the three-qubit Hadamard operator
>>> state1 = qc.bitstring(0, 1, 0)  # the state |010⟩
>>> state2 = qc.bitstring(1, 0, 0)  # the state |100⟩
>>> rho = qc.DensityOperator.from_ensemble(
...     [state1, state2],
...     ps=[0.3, 0.7]
... )  # the statistical ensemble known to be in states |010⟩ and |100⟩
...    # each with probability 1/2.
>>> result = H(rho)  # apply the operator to the state

Under the hood, this takes the density operator \(\rho\) to \(H \rho H^{\dagger}\), which is the resulting density operator if the operator \(H\) is applied to the state whose uncertainty is being represented.

If the operator is for a lower-dimensional space than the density operator represents, or if you want to permute the roles of the qubits in the operator application, a qubit_indices argument can be supplied as when applying operators to pure states:

>>> H = qc.Hadamard(d=2)  # the two-qubit Hadamard operator
>>> state1 = qc.bitstring(0, 1, 0)  # the state |010⟩
>>> state2 = qc.bitstring(1, 0, 0)  # the state |100⟩
>>> rho = qc.DensityOperator.from_ensemble(
...     [state1, state2],
...     ps=[0.3, 0.7]
... )  # the statistical ensemble known to be in states |010⟩ and |100⟩
...    # each with probability 1/2.
>>> result = H(rho, qubit_indices=[0, 2])  # apply the two-qubit operator to
...                                        # qubits 0 and 2

Measurement in the computational basis may be performed on one or more qubits of a mixed state with the DensityOperator.measure() method, which returns the result of measurement.

>>> state1 = qc.zeros(1)  # the state |0⟩
>>> state2 = qc.ones(1)   # the state |1⟩
>>> mixed_state = qc.DensityOperator.from_ensemble(
...     [state1, state2],
...     ps=[0.5, 0.5]
... )  # the statistical ensemble known to be in states |0⟩ and |1⟩
...    # each with probability 1/2.
>>> mixed_state.measure()  # measure the qubit
(0,)

The DensityOperator.measure() method has two arguments: a list of indices specifying the qubits to be measured, and a flag specifying whether the measured qubits are to be removed from the state. If no indices are supplied, every qubit is measured. E.g., in the following example, the system is known to be in state \(|00⟩\) or state \(|11⟩\) with equal probability. We measure the first qubit, opting to then remove the collapsed qubit from the state, and we measure \(0\). Post-measurement, the density operator represents the state \(|0⟩\) (for the second qubit) with certainty.

>>> state1 = qc.bitstring(0, 0)  # the state |00⟩
>>> state2 = qc.bitstring(1, 1)  # the state |11⟩
>>> mixed_state = qc.DensityOperator.from_ensemble(
...     [state1, state2],
...     ps=[0.5, 0.5]
... )  # the statistical ensemble known to be in states |00⟩ and |11⟩
...    # each with probability 1/2.
>>> mixed_state.measure(qubit_indices=[0], remove=True)  # measure the state, and remove the
...                                                      # qubit from the density operator
(0,)
>>> print(mixed_state)
... Density operator for 1-qubit state space. Tensor:
    [[1.+0.j 0.+0.j]
     [0.+0.j 0.+0.j]]

Generally, the post-measurement state is as follows. Each of the ensemble of states the density operator represents collapses to the computational basis state corresponding to the result of the measurement of the measured qubits. The mixture probabilities, representing our state of belief of the current state of the system, are updated using Bayes rule. Equivalently, the post-measurement density operator is \(P_m \rho P_m^{\dagger} / p(m)\), where \(m\) is the outcome of measurement, \(p(m)\) is the probability of that outcome, and \(P_m\) is the projector onto the computational basis states corresponding to the measurement outcome, and \(\rho\) is the pre-measurement state.

One can take the tensor product of density operators for subsystems, in the same way as for states and operators, to obtain the density operator for the larger composite system:

>>> state1 = qc.DensityOperator.from_ensemble(list_of_states1)
>>> state2 = qc.DensityOperator.from_ensemble(list_of_states2)
>>> print(state1.rank)
2
>>> print(state2.rank)
4
>>> composite_system = state1 * state2
>>> print(composite_system.rank)
6

Reduced Density Operators and Purification

The reduced density operator of a state can be computed using the State.reduced_density_operator() or DensityOperator.reduced_density_operator() methods, for pure and mixed states respectively. This will give the reduced density operator for the subsystem of the given qubit indices, by tracing out the remaining qubits. This is equivalent to the density operator obtained by measuring the traced out qubits and forgetting the result of the measurement.

>>> state = qc.bell_state(0, 0)
>>> rho = state.reduced_density_operator(qubit_indices=[0])
>>> print(rho)
Density operator for 1-qubit state space. Tensor:
[[0.5+0.j 0. +0.j]
 [0. +0.j 0.5+0.j]]

A mixed state can be purified with the DensityOperator.purify() method. This will produce a 2d qubit pure state from a d qubit mixed state, with the measurement probabilities of the first d qubits in agreement with the measurement probabilities of the mixed state.

>>> state = rho.purify()
>>> print(state)
State([[0.70710678+0.j 0.        +0.j]
       [0.        +0.j 0.70710678+0.j]])

Warning: The No-Cloning Theorem

The no cloning theorem says that, in general, given a quantum system in a given state, one cannot clone the state such that another system is in the same state without modifying the state of the original system.

There are several ways to violate the no-cloning theorem in the QCircuits library, and users should be aware that none of the following constitute valid steps in a quantum computation.

Firstly, one could simply create a deepcopy of a state.

Secondly, one could extract the underlying tensor of a state and then create another state from it.

Thirdly, one could apply an operator to a state, such as:

>>> y = Op(x)

and then use both \(x\) and \(y\) in subsequent computation. It is clear that this violates the no-cloning theorem, since the operator may be the identity.

Fourthly, one can take the tensor product of a state with itself, or a tensor power, as in:

>>> x * x
>>> x**5

I have decided not to disallow any of these, since each has genuine use-cases. For example if we know how to prepare a single qubit in a given state, then we know how to prepare a multi-qubit system with each qubit in that state, and the tensor power is a short-hand way of achieving this.

Entanglement / Schmidt Number

A state \(|\phi⟩\) of a composite quantum system \(A\otimes B\) has a Schmidt decomposition: \(|\phi⟩ = \sum_i \lambda_i |i_A⟩|i_B⟩\), where the Schmidt coefficients \(\lambda_i\) are non-negative and the \(|i_A⟩\) and \(|i_B⟩\) are orthonormal bases for systems \(A\) and \(B\). The number of non-zero Schmidt coefficients is a measure of the entanglement of the two systems, and is called the Schmidt number.

In QCircuits, the Schmidt number of a multi-qubit state may be computed with respect to a partitioning of the qubits into two subsystems with the State.schmidt_number() method, by supplying a list of qubit indices for one of the subsystems.

E.g., the Bell State \((|00⟩ + |11⟩)/\sqrt{2}\), which is already written as a Schmidt decomposition, has Schmidt number 2:

>>> print(qc.bell_state(0, 0).schmidt_number(indices=[0]))
2

The state \((|00⟩ + |01⟩ + |10⟩ + |11⟩)/2\) can be written as a product \((|0⟩ + |1⟩)/\sqrt{2} \otimes (|0⟩ + |1⟩)/\sqrt{2}\), so it has a Schmidt number of 1:

>>> print(qc.positive_superposition(d=2).schmidt_number(indices=[0]))
1

Any state that is the result of a tensor product of states with respect to two subsystems with have a Schmidt number 1 with respect to those two subsystems:

>>> state_A = qc.bitstring(0, 1, 0)
>>> state_B = qc.bitstring(0, 0, 0, 1)
>>> state_AB = state_A * state_B
>>> print(state_AB.schmidt_number(indices=[0, 1, 2]))
1

Examples

For examples of the use of QCircuits to implement quantum algorithms, see the examples page.