# Matrix representation of quantum circuits - notations and gotchas
> A case study by building a tensor network to match qiskit conventions

- toc: true 
- badges: true
- comments: true
- categories: [qiskit, tensor networks, quantum concepts]
- image: images/upsidedown.jpg

# Intro
Usually, for experimenting with quantum circuits I use `qiskit`. As any higher level environment it is very convenient for common tasks, but may turn out too inflexible for unusual use cases. A somewhat opposite approach is to use much lower level tools to gain in flexibility at the expense of convenience. Currently I want to use Google's `tensornetwork` [package](https://tensornetwork.readthedocs.io/en/latest/) for simulations and training of quantum circuits, but this requires building many things that are for free in `qiskit` from scratch. It is also necessary to become explicit about conventions for matrix representation of quantum circuits. As long as you stay within a single framework this may not be an issue. However for debugging purposes as well as for comparison between different frameworks this may become unavoidable. Thus, I always anticipated, that a day will come when I need to face my fears and order all terms in a tensor product by hands. Now it seems I'm past the difficult part and I'm better writing this down in case I would need to do something similar in the future.

## Defining the problem
OK, so what is the problem? Consider the following simple circuit built with `qiskit`:


```python
#collapse
import numpy as np
from qiskit import QuantumCircuit
from qiskit.quantum_info import Operator, Statevector

qc = QuantumCircuit(2)
qc.x(0)
qc.y(1)
qc.cx(0,1)
qc.draw(output='mpl')
```

It is not hard or ambiguous to interpret what this circuit does by inspecting the diagram. Say the input state is $q_0=|0\rangle$, $q_1=|1\rangle$. After $X$ acts on $q_0$ it becomes $q_0\to X |0\rangle=|1\rangle$. Similarly, $q_1$ after $Y$ becomes $q_1\to Y|1\rangle=-i |0\rangle$. Since now $q_0$ is "on" the CNOT gate switches the state of $q_1$ further to $q_0 \to -i|1\rangle$. So the end result is that $q_0=|0\rangle, q_1=|1\rangle$ is transformed to $q_0=|1\rangle, q_1=-i|1\rangle$. Or perhaps a picture says it better



Similarly, we can work out what the circuit does for other computational basis states which by linearity fully fixes the action of the circuit. Although quite explicit, this is a clumsy description. This is why the matrix notation is usually used. And indeed, we can obtain the matrix corresponding to our quantum circuit quite easily in `qiskit`:


```python
U_qs = Operator(qc).data
U_qs
```




    array([[0.+0.j, 0.+0.j, 0.+0.j, 0.-1.j],
           [0.+1.j, 0.+0.j, 0.+0.j, 0.+0.j],
           [0.+0.j, 0.+1.j, 0.+0.j, 0.+0.j],
           [0.+0.j, 0.+0.j, 0.-1.j, 0.+0.j]])



It is important to realize that a number of conventions must be chosen before such explicit matrix representation can be written down. In particular, I will emphasize two points I tripped over while studying this: ordering of the qubit states in the tensor product or "vertical ordering" and ordering of operators or "horizontal ordering". 



In the rest of the post I will clarify what are the conventions used in `qiskit` and how to reproduce the circuit with the `tensornetwork` library.

# States: vertical ordering

## Single qubit states
First we need to give matrix representations to two basis states of a single qubit. Here I think it is quite uncontroversial to choose
\begin{align}
|0\rangle = \begin{pmatrix}1\\0\end{pmatrix},\qquad |1\rangle = \begin{pmatrix}0\\1\end{pmatrix} \label{kets}
\end{align}

These are the "ket" vectors. Their "bra" counterparts are
\begin{align}
\langle 0| = \begin{pmatrix}1 & 0\end{pmatrix}, \qquad \langle 1| = \begin{pmatrix}0 & 1\end{pmatrix} \label{bras}
\end{align}

With these, the following operators can be computed
\begin{align}
|0\rangle\langle 0| = \begin{pmatrix}1 & 0 \\ 0 & 0\end{pmatrix},\qquad |0\rangle\langle 1| = \begin{pmatrix}0 & 1 \\ 0 & 0\end{pmatrix} \nonumber\\ |1\rangle\langle 0| = \begin{pmatrix}0 & 0 \\ 1 & 0\end{pmatrix},\qquad |1\rangle\langle 1| = \begin{pmatrix}0 & 0 \\ 0 & 1\end{pmatrix} \label{ketbras}
\end{align}
## Multiple qubit states
When there is more than a single qubit things become a bit more interesting and potentially confusing. For example, the combined Hilbert space of two qubits $\mathcal{H}_2$ is a tensor product of single-qubit Hilbert spaces $\mathcal{H}_2 = \mathcal{H}_1 \otimes \mathcal{H}_1$ but we need to decide which qubit goes first and which goes second. In `qiskit` a convention is adopted that additional qubits join from the *left*, i.e. when we have two qubits as here


```python
#collapse
qc01 = QuantumCircuit(2)
qc01.draw(output='mpl')
```


The state of the system is $|q_1\rangle\otimes |q_0\rangle$ (this is of course only true literally for [non-entangled states](https://idnm.github.io/blog/quantum%20concepts/qiskit/2021/07/12/Entanglement.html) but we can define everything only on the computational basis states ). OK, but how do we translate this into the matrix representation? The states in the tensor product of vector spaces can be represented by the [Kronecker product](https://en.wikipedia.org/wiki/Kronecker_product) which is not symmetric with respect to permutation arguments. Best way to explain how Kronecker product works is, as usual, through examples:

\begin{align}
\begin{pmatrix} 1 \\ 0 \end{pmatrix} \otimes \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} a \\ b \\ 0 \\ 0 \end{pmatrix},\qquad \begin{pmatrix} 0 \\ 1 \end{pmatrix} \otimes \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} 0\\ 0\\ a \\ b \end{pmatrix}
\end{align}
Result for generic left vector can be obtained by linearity
\begin{align}
\begin{pmatrix} x \\ y \end{pmatrix} \otimes \begin{pmatrix} a \\ b \end{pmatrix} = x \begin{pmatrix} 1 \\ 0 \end{pmatrix} \otimes \begin{pmatrix} a \\ b\end{pmatrix} +y\begin{pmatrix} 0 \\ 1 \end{pmatrix} \otimes \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} x a\\ x b\\ y a \\ y b \end{pmatrix}  = \begin{pmatrix} x \begin{pmatrix} a\\ b\end{pmatrix} \\ y  \begin{pmatrix} a \\  b\end{pmatrix} \end{pmatrix}
\end{align}

The last notation here is a bit informal but it shows what happens. One just substitutes the right vector into all elements of the left vector, multiplied by the corresponding components of the left vector. The Kronecker product is defined in the same way for matrices of arbitrary size, not just for two vectors.

So, now we can compute matrix representations of states in the computation basis of two-qubit system

\begin{align}
|00\rangle = \begin{pmatrix}1\\0 \end{pmatrix} \otimes \begin{pmatrix}1\\0 \end{pmatrix} = \begin{pmatrix}1\\0\\0\\0\end{pmatrix},\quad |01\rangle = \begin{pmatrix}1\\0 \end{pmatrix} \otimes \begin{pmatrix}0\\1 \end{pmatrix} = \begin{pmatrix}0\\1\\0\\0\end{pmatrix} \label{01}\\
|10\rangle = \begin{pmatrix}0\\1\end{pmatrix} \otimes \begin{pmatrix}1\\0 \end{pmatrix} = \begin{pmatrix}0\\0\\1\\0\end{pmatrix},\quad |11\rangle = \begin{pmatrix}0\\1\end{pmatrix} \otimes \begin{pmatrix}0\\1 \end{pmatrix} = \begin{pmatrix}0\\0\\0\\1\end{pmatrix}
\end{align}

There is a useful relation between the index of the non-zero element $n$ in the four-dimensional representation and the computational basis bitstring $q_1q_0$, namely $n=2q_1+q_0$. I.e. the bitstring $q_1q_0$ is the binary representation of the index $n$. This extends to arbitrary number of qubits, for example since $101$ is $5$ in binary representation it follows
\begin{align}
|101\rangle = \begin{pmatrix}0\\0\\0\\0\\0\\1\\0\\0 \end{pmatrix} \label{101}
\end{align}
(try to obtain this from the two tensor products!)

Don't believe me? OK, let's check! In `qiskit` there is a convenient function to construct a vector representation from a bit string which we will take advantage of. First start with a two-qubit example:


```python
s01 = Statevector.from_label('01')
s01.data
```




    array([0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j])



 Comparing to \eqref{01} we find agreement. Similarly, 


```python
s101 = Statevector.from_label('101')
s101.data
```




    array([0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j])



Again, this is in agreement with \eqref{101}. 

However, I am not sure that this relation is sufficient to justify the ordering of the tensor products. To me it is much more natural to read the circuit from top to bottom and construct the Hilbert spaces accordingly, say $\mathcal{H}_0\otimes \mathcal{H}_1 \otimes \mathcal{H}_2 \dots$ instead of $\cdots \mathcal{H}_2\otimes \mathcal{H}_1\otimes \mathcal{H}_0$. Later I will change the ordering of the tensor product to my liking, but for now we stick with the `qiskit` one. Now, with conventions for states in place we can proceed to operators.

# Operators: horizontal ordering
One can say that convention for states representation and ordering of tensor products is a "vertical" convention. There is also a "horizontal" convention which might be potentially confusing. Consider the following circuit


```python
#collapse
qc123 = QuantumCircuit(1)
qc123.rx(1, 0)
qc123.ry(2, 0)
qc123.rz(3, 0)
qc123.draw(output='mpl')
```

Here, the operator $R_x$ is appplied first, the operator $R_y$ second and $R_z$ last. So in mathematical notation the circuit corresponds to $R_z R_y R_x$ and *not* to $R_x R_y R_z$. I think that the circuit notation is actually better. We think and write from left to right, this is also a direction that time flows on paper. When another thing happens, we write it to the right and it would be convenient to apply the corresponding operator also to the right. I heard real mathematicians complain about that issue, but I guess we are stuck with it for now.

# Paper-and-pencil computation
With the set up in place we can compute the circuit of interest by hands. For convenience I plot it here once again:


```python
#collapse
qc.draw(output='mpl')
```

OK, so what is the unitary matrix corresponding to this circuit? It is 
\begin{align}
U = CNOT_{01} \cdot (Y\otimes X)
\end{align}
Here
\begin{multline}
CNOT_{01} = \mathbb{1}\otimes |0\rangle\langle 0|+X\otimes |1\rangle\langle 1|=\\\begin{pmatrix}1&0\\0&1\end{pmatrix}\otimes \begin{pmatrix}1&0\\0&0\end{pmatrix}+\begin{pmatrix}0&1\\1&0\end{pmatrix}\otimes \begin{pmatrix}0&0\\0&1\end{pmatrix}=\begin{pmatrix}1&0&0&0\\0&0&0&1\\0&0&1&0\\0&1&0&0\end{pmatrix}
\end{multline}
and
\begin{align}
Y\otimes X = \begin{pmatrix} 0& -i\\i&0\end{pmatrix} \otimes \begin{pmatrix} 0& 1\\1&0\end{pmatrix}=\begin{pmatrix}0&0&0&-i\\0&0&-i&0\\0&i&0&0\\i&0&0&0\end{pmatrix}
\end{align}
Multiplying them together gives
\begin{align}
U = \begin{pmatrix}0 & 0 & 0 & -i \\ i&0&0&0 \\ 0 & i & 0 & 0 \\ 0 & 0 & -i & 0\end{pmatrix}
\end{align}
Alright, so this is indeed the matrix that `qiskit` computes:


```python
U_qs
```




    array([[0.+0.j, 0.+0.j, 0.+0.j, 0.-1.j],
           [0.+1.j, 0.+0.j, 0.+0.j, 0.+0.j],
           [0.+0.j, 0.+1.j, 0.+0.j, 0.+0.j],
           [0.+0.j, 0.+0.j, 0.-1.j, 0.+0.j]])



We can now also check that that the states evolve as we expected. For example recall that we computed that our quantum circuit maps $q_0 =|0\rangle, q_1 =|1\rangle$ to $q_0 =|1\rangle, q_1 =|1\rangle$ with an overall phase $-i$. Agreement with  `qiskit` can be checked as follows:


```python
qs_state = Statevector.from_label('10').evolve(qc).data 
our_state = -1j*Statevector.from_label('11').data

np.allclose(qs_state, our_state)
```




    True



# Implementation with `tensornetworks`
I will not give a proper introduction to tensor networks but just make some digressions I think should be helpful as we go along.

First thing we will need are the matrices defining $X, Y$ and $CNOT$ gates. Let us introduce them.


```python
X = np.array([[0, 1], [1, 0]])
Y = np.array([[0, -1j], [1j, 0]])
CNOT = np.array([[1, 0, 0, 0],
                 [0, 0, 0, 1],
                 [0, 0, 1, 0],
                 [0, 1, 0, 0]]).reshape(2,2,2,2)
```

## An aside about reshaping
Note that as usually written, $CNOT$ is a $4\times4$ matrix. Since as a quantum gate it acts on two qubits, so it should rather be a four-legged tensor. This is the purpose of the reshaping operation. At first the reshaping might be a bit tricky, so let me illustrate it with an example. Introduce two $4\times4$ matrices and define their product:


```python
A = np.random.rand(4,4)
B = np.random.rand(4,4)

AB = A @ B
```

Now define the corresponding four-legged tensors.


```python
import tensornetwork as tn

a = tn.Node(A.reshape(2,2,2,2))
b = tn.Node(B.reshape(2,2,2,2))
```

By contracting the legs (or "edges" in terminology of `tensornetworks`) appropriately, we can reproduce the matrix multiplication. First the code:


```python
a[2] ^ b[0]
a[3] ^ b[1]

ab = tn.contractors.greedy([a, b], output_edge_order=[a[0], a[1], b[2], b[3]]).tensor
```

We can check that the contraction performed in this way exactly reproduces the matrix multiplication of original $4\times4$ matrices:


```python
np.allclose(AB, ab.reshape(4,4))
```




    True



This can be interpreted graphically as follows. First, the reshaping procedure can be thought of as splitting each of two four-dimensional legs of the original matrix into two two-dimensional ones



The labels on the legs have nothing to do with qubit states, these are just indices of edges as assigned by `tn.Node` operation on our matrices. The matrix multiplication of the original matrices in terms of four-legged tensors then can be drawn as follows



The index arrangements in the last part explain why we connected the edges in our code the way we did. This is something to watch out for. For example, connecting edges of two identity tensors in the wrong way may produce a $SWAP$ gate.

## Tensor product ordering
The matrix representation of a tensor diagram like this



also comes with a convention for the ordering of tensor products. In `tensornetwork` as well as in my opinion it is natural to order top-down, i.e. the above diagram is $U\otimes \mathbb{1}$ instead of $\mathbb{1}\otimes U$ as is adopted in `qiskit`. 


## Circuit from tensor network
Alright, not we are in a position to reproduce the circuit unitary from the tensor network with nodes `x`, `y` and `cnot`:


```python

# Make tensors from matrices
x, y, cnot = list(map(tn.Node, [X, Y, CNOT]))

# Connect edges properly
cnot[2] ^ y[0]
cnot[3] ^ x[0]

# Perform the contraction ~ matrix multiplication
U_tn = tn.contractors.greedy([cnot, x, y], output_edge_order=[cnot[0], cnot[1], y[1], x[1]]).tensor
```

This way of contracting the edges corresponds to the following diagram:



Note that this is basically the original circuit with both the vertical and the horizontal directions reversed. The horizontal reversal is due to mathematical vs circuit notation (circuit is better!) and the vertical reversal is due to the mismatch between `qiskit` and `tensornetwork` ordering of tensor product (`tensornetwork`'s is better!). We can check that the unitary we obtain from this tensor network agrees with `qiskit`'s


```python
np.allclose(U_tn.reshape(4,4), U_qs)
```




    True



## A better way
I find all this misalignment very inconvenient and hard to debug. Ideally I want to look at the quantum circuit and construct the corresponding tensor network just as I read a text: from left to right and from top to bottom. Here I propose a solution which seems much more satisfactory to me. We will deal with horizontal reversal by first defining edges and then applying gates to them. This way we can read the circuit from left to right and simply add new gates, just as in `qiskit`. I will not try to revert the vertical direction directly, because I find it hard to think upside down. Instead, for comparison with `qiskit` I will use a built-in `reverse_bits` method.

So let's start by defining a function that applies a given gate to the collection of qubits (this is a slight modification of an [example](https://tensornetwork.readthedocs.io/en/latest/quantum_circuit.html) from `tensornetwork` docs) :


```python
def apply_gate(qubits, gate_tensor, positions):
    
    gate = tn.Node(gate_tensor)
    assert len(gate.edges) == 2*len(positions), 'Gate size does not match positions provided.'
    

    for i, p in enumerate(positions):
        # Connect RIGHT legs of the gate to the active qubits
        gate[i+len(positions)] ^ qubits[p] 
        # Reassing active qubits to the corresponding LEFT legs of the gate 
        qubits[p] = gate[i]
```

Importantly, here, in contrast to the official docs, we append the gate from the *left*, so that a sequence of application of some $G_1$ followed by $G_2$ is equivalent to the application of $G_2\cdot G_1$. Now there is one more subtlety. Previously we used matrix representation of $CNOT$ assuming that the uppermost qubit comes last in the tensor product. Now that we decided to turn this convention upside down our matrix representation of $CNOT$ must be $CNOT =|0\rangle\langle 0|\otimes \mathbb{1}+|1\rangle\langle 1|\otimes X$ or explicitly


```python
CNOT = np.array([[1, 0, 0, 0],
                 [0, 1, 0, 0],
                 [0, 0, 0, 1],
                 [0, 0, 1, 0]]).reshape(2,2,2,2)
```

With that we are ready to reconstruct our original circuit in a convenient way:


```python

# The context manager  `NodeCollection` is a bit of a magic trick
# which keeps track of all tensors in the network automatically.

all_nodes = []
with tn.NodeCollection(all_nodes): 
    # I do not know how to create 'abstract' edges in `tensornetworks`.
    # Instead, I create an identity tensor and use its edges to apply new gates to.
    id0 = tn.Node(np.identity(4).reshape(2,2,2,2))
    qubits0 = id0.edges[2:4]
    qubits = id0.edges[0:2]

    apply_gate(qubits, X, [0])
    apply_gate(qubits, Y, [1])
    apply_gate(qubits, CNOT, [0,1])
```

Now let us check!


```python
U_tn = tn.contractors.greedy(all_nodes, output_edge_order=qubits+qubits0).tensor.reshape(4,4)
U_reversed_qs = Operator(qc.reverse_bits()).data

np.allclose(U_tn, U_reversed_qs)
```




    True



Wohoo, it worked! If that looked simple to you I'm happy. It took me several hours of debugging to finally match the two matrices. Just to make sure, let me conclude with a more complicated example.


```python
qc3 = QuantumCircuit(3)
qc3.x(0)
qc3.cx(0, 1)
qc3.y(1)
qc3.x(2)
qc3.cx(2, 1)
qc3.y(2)

qc3.draw(output='mpl')
```

As you can see, constructing the tensor network analog now works more or less identically:


```python
all_nodes = []
with tn.NodeCollection(all_nodes):
    id0 = tn.Node(np.identity(8).reshape(2,2,2,2,2,2))
    qubits0 = id0.edges[3:6]
    qubits = id0.edges[0:3]
    
    # The essential part
    apply_gate(qubits, X, [0])
    apply_gate(qubits, CNOT, [0, 1])
    apply_gate(qubits, Y, [1])
    apply_gate(qubits, X, [2])
    apply_gate(qubits, CNOT, [2, 1])
    apply_gate(qubits, Y, [2])
```

And now we compare:


```python
U3_tn = tn.contractors.greedy(all_nodes, output_edge_order=qubits+qubits0).tensor.reshape(8,8)
U3_qs_reversed = Operator(qc3.reverse_bits()).data

np.allclose(U3_tn, U3_qs_reversed)
```




    True



Alright, this resounding **True** is the best way to conclude that comes to mind. I own many thanks to [Ilia Luchnikov](https://github.com/LuchnikovI) for the help with `tensornetwork` library. Any questions are welcome in the comments!

5e5c55b04e2aef9e44ba5005c80efa9b2f51ea3c
ipynb
Jupyter Notebook
_notebooks/2021-08-18-Matrix representation of quantum circuits.ipynb
idnm/blog
a9e976ea45fe077b7b13a5fa3680fab1affc2c48
Apache-2.0
_notebooks/2021-08-18-Matrix representation of quantum circuits.ipynb
idnm/blog
a9e976ea45fe077b7b13a5fa3680fab1affc2c48
Apache-2.0
_notebooks/2021-08-18-Matrix representation of quantum circuits.ipynb
idnm/blog
a9e976ea45fe077b7b13a5fa3680fab1affc2c48
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES