Decomposing unitary matrix into Q# quantum gates

#	User	Rating
1	tourist	4009
2	jiangly	3823
3	Benq	3738
4	Radewoosh	3633
5	jqdai0815	3620
6	orzdevinwang	3529
7	ecnerwala	3446
8	Um_nik	3396
9	ksun48	3390
10	gamegame	3386

#	User	Contrib.
1	cry	167
2	Um_nik	163
3	maomao90	162
3	atcoder_official	162
5	adamant	159
6	-is-this-fft-	158
7	awoo	157
8	TheScrasse	154
9	Dominater069	153
9	nor	153

This post is part of the Q# Advent Calendar 2020. Check it out for more great posts about quantum computing and Q#!

Hello Codeforces!

In this post I am going to tell you about an algorithm which implements arbitrary unitary matrix as a sequence of elementary quantum gates in Q# programming language.

Introduction

Algorithms for a quantum computer are usually described as a sequence of quantum gates (such as Pauli X, Hadamard or CNOT) acting on register of qubits, followed by a sequence of measurements. Any quantum gate can be described as a unitary matrix acting on a vector representing the current state of a quantum computer. Applying multiple gates multiplies state by a matrix which is matrix product of matrices for individual gates. Therefore, any sequence of quantum gates can be described by a unitary matrix.

But what about going the other way: given a unitary matrix, represent it as a sequence of quantum gates? Turns out it's possible if we can use gates from universal set. One famous result is that you can implement any unitary operation by using only one-qubit gates and CNOT gate (which means these gates form a universal set).

In this post we will consider different universal set of gates: set of gates which you can describe using Q# language. In other words, I am going to describe an algorithm, which for arbitrary unitary matrix produces Q# source code implementing action of that matrix on a register of qubits. This algorithm was inspired by problem B2 from Microsoft Q# Coding Contest - Winter 2019, to solve which I had to write Q# code implementing a particular 4x4 unitary matrix.

What operations does Q# support?

Firstly, Q# has a variety of 1-qubit gates. We will use X gate and rotation gates (R1, Ry, Rz).

Here are matrices of these gates:

$$$ X = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, R_1(\alpha) = \begin{pmatrix} 1 & 0 \\ 0 & e^{i\alpha} \end{pmatrix}, $$$ $$$ R_y(\alpha) = \begin{pmatrix} \cos(\alpha/2) & \sin(\alpha/2) \\ -\sin(\alpha/2) & \cos(\alpha/2) \end{pmatrix}, $$$ $$$ R_z(\alpha) = \begin{pmatrix} e^{i \alpha/2} & 0 \\ 0 & e^{-i \alpha/2} \end{pmatrix}. $$$

Secondly, you can apply a controlled gate. What's a controlled gate? Let's say you have a gate $$$A$$$ acting on qubits $$$q_1, \dots, q_N$$$. Gate A controlled by qubits ($$$p_1, \dots, p_M$$$) is such a gate that acts as $$$A$$$ if all qubits $$$p_i$$$ are in state $$$|1 \rangle$$$, and does nothing otherwise. In Q# it's done using the "Controlled" keyword. For example, Controlled X([control1, control2], target) is equivalent to CCNOT(control1, control2, target).

Of course, Q# has much more features, but we will use only these.

Now, when we know what gates we are allowed to use, we have problem of matrix decomposition: represent given unitary matrix $$$A$$$ of size $$$2^N \times 2^N$$$ as product of matrices representing Q# X, R1, Ry and Rz gates, possibly controlled. We are going to solve this problem in 4 steps.

Step 1. Two-level decomposition

Two-level unitary matrix is a unitary matrix obtained from an identity matrix by changing a 2×2 principal submatrix. In other words, take identity matrix $$$A$$$, select distinct indices $$$i, j$$$ and replace elements $$$A_{i,i}, A_{i,j}, A_{j,i}, A_{j,j}$$$ such that $$$A$$$ remains unitary — you'll get two-level unitary matrix. Any unitary matrix can be decomposed into a sequence of two-level unitary matrices. In this step our goal is to represent unitary matrix $$$A$$$ as a product of two-level unitary matrices.

We are going to achieve this by transforming the original matrix $$$A$$$ row by row. Our goal is to make all elements outside main diagonal zeroes, and we are allowed to do this by multiplying the matrix by a two-level unitary matrix (from the right). This guarantees that the matrix remains unitary. When we zeroed all elements on the first row (except the first one), all elements in the first column will also become zeroes (this is a property of unitary matrices). Then we can forget about the first row and column and repeat the process for the next column. We should repeat this $$$2^N - 2$$$ times, and it will look like this:

$$$ \begin{pmatrix} * & * & * & * \\ * & * & * & * \\ * & * & * & * \\ * & * & * & * \end{pmatrix} \to \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & * & * & * \\ 0 & * & * & * \\ 0 & * & * & * \end{pmatrix} \to \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & * & * \\ 0 & 0 & * & * \end{pmatrix} $$$

In the end we will see that final matrix $$$U_d$$$ if 2-level unitary. Let's denote $$$U_1, U_2, \dots U_D$$$ — all matrices which we applied during this algorithm. Then the whole process can be written as $$$ A \cdot U_1 \cdot U_2 \dots U_D = U_f$$$, from which follows $$$ A= U_f \cdot U_D^+ \cdot U_{D-1}^+ \dots U_1^+$$$, which is desired decomposition ($$$U^+$$$ denotes adjoint matrix, which for unitary matrix is the same as inverse).

How do we turn a matrix element into zero without affecting elements to the top and to the right, and only by multiplying matrix by 2-level unitary matrix? Let's focus only on two rows and denote desired two-level unitary $$$U$$$, then what we want is to find such U that:

$$$ \begin{pmatrix} \dots & a & b & 0 & \dots & 0 \\ \dots & \dots & \dots & \dots &\dots & \dots \end{pmatrix} ~ U = \begin{pmatrix} \dots & c & 0 & 0 & \dots & 0 \\ \dots & \dots & \dots & \dots &\dots & \dots \end{pmatrix}. $$$

$$$U$$$ can be two-level unitary matrix with non-trivial 2x2 submatrix $$$U'$$$, such that $$$(a ~ b ) U' = ( c ~ 0 )$$$.

Let's use the following ansatz:

$$$ U' = \begin{pmatrix} \cos{\theta}e^{i \lambda} & \sin{\theta}e^{i \mu} \\ -\sin{\theta}e^{i \mu} & \cos{\theta}e^{-i \lambda} \end{pmatrix}. $$$

Then one can show that fitting values of real parameters $$$\theta, \lambda, \mu$$$ are $$$\theta = \arctan\left(\left| \frac{b}{a} \right| \right);~ \lambda = - \arg(a); ~\mu = \pi + \arg(b).$$$

The one special case is when $$$a = 0$$$. In this case we should use the Pauli X matrix as $$$U'$$$.

Step 2. Gray codes

Now all two-level matrices in decomposition act on pairs of states $$$(i, i+1)$$$. For our purposes we want them to act on pairs of states differing only in one bit, i.e. $$$(i, i \oplus 2^k)$$$.

Luckily, for any positive integer $$$n$$$ there exists such permutation of numbers $$$0, 1, \dots, 2^n-1$$$, that any two adjacent numbers in it differ only in one bit. Such permutation is called binary-reflected Gray code, and is given by formula $$$\pi_i = i \oplus \left\lfloor i/2 \right\rfloor$$$, where $$$i=0, 1, \dots, 2^n-1$$$. For example, Gray code for $$$n=3$$$ is $$$(0, 1, 3, 2, 6, 7, 5, 4)$$$.

So, let's re-number basis states such that states which differed in one bit now are adjacent (using Gray code). Then let's apply decomposition from step 1 and go back to original basis numbering.

More formally, define a permutation matrix $$$P$$$ of size $$$2^N \times 2^N$$$, such that $$$P_{ij} = \delta_{i, \pi_j}$$$. Define $$$A' = P A P^+$$$. Apply two-level decomposition algorithm to $$$A'$$$, which results in decomposition $$$A' = \prod_{i}A'_i$$$, where $$$A'_i$$$ acts on states whose indices differ in 1. Then $$$A = P^+ A' P = \prod_{i}(P^+ A'_i P)$$$. Each matrix $$$(P^+ A'_i P)$$$ acts on states differing in 1 bit, which is exactly what we need.

Step 3. Implementing 2-level matrix as fully-controlled 1-qubit gate

Let's call a gate “fully controlled gate acting on qubit $$$i$$$” if this gate acts on a qubit $$$i$$$ and is controlled by all other qubits in the register. This gate will act on a certain basis state only if all bits in the index of this state (except maybe $$$i$$$-th) are set to one.

Let $$$U$$$ be a two-level unitary acting matrix on states $$$(i, i \oplus 2^r$$$). Define $$$J_0$$$ — set of all indices $$$j$$$, such that $$$j$$$-th bit of $$$i$$$ is zero, $$$J_1$$$ — set of all indices $$$j$$$, such that $$$j$$$-th bit of $$$i$$$ is one (both sets don't include $$$r$$$). Then we need to apply this two-level unitary only to pair of such states, whose indices have zeroes on positions $$$J_0$$$, and ones on positions $$$J_1$$$. If $$$J_0 = \emptyset$$$, this is simply a fully-controlled gate acting on qubit $$$r$$$.

But if there is some $$$j \in J_0$$$, then we have to apply $$$X$$$ on $$$j$$$-th qubit, then apply $$$U$$$ and then apply $$$X$$$ on $$$j$$$-th qubit again. This will work because $$$X$$$ acting on $$$j$$$-th qubit swaps state $$$i$$$ with state $$$i \oplus 2^j$$$, so if $$$U$$$ does something only with states $$$i$$$ where $$$i[j]=1$$$, then $$$XUX$$$ does the same thing with states $$$i$$$ where $$$i[j]=0$$$.

So, to implement two-level unitary matrix $$$U$$$ we have to apply $$$X$$$ gate to all qubits from $$$J_0$$$, then apply fully-controlled gate on qubit $$$r$$$ and then again apply $$$X$$$ gate to all qubits from $$$J_0$$$.

Step 4. Implementing 1-qubit gate with rotation gates

The only thing left is to decompose arbitrary 1-qubit gate into gates supported by Q#. This decomposition is given by the formula:

$$$ U = R_1(\phi)R_z(\lambda+\mu) R_y(2 \theta) R_z(\lambda — \mu), $$$

where $$$\phi = \arg{\det U}$$$, $$$\theta = \arccos(|U'_{00}|)$$$, $$$\lambda = \arg(U'_{00})$$$, $$$\mu = \arg(U'_{01}) $$$, $$$U' = R_1(- \phi) U$$$.

Example

Let's consider the following matrix:

$$$A = \frac{1}{\sqrt{3}} \begin{pmatrix} 1 & 1 & 1 & 0 \\ 1 & e^{\frac{2\pi i}{3}} & e^{\frac{4 \pi i}{3}} & 0 \\ 1 & e^{\frac{4\pi i}{3}} & e^{\frac{2 \pi i}{3}} & 0 \\ 0 & 0 & 0 & -i \sqrt{3} \end{pmatrix} $$$

Implementing such matrix in Q# was one way to solve this problem. If you feed it to the algorithm, it will produce the following Q# code:

operation ApplyUnitaryMatrix (qs : Qubit[]) : Unit {
  CNOT(qs[1], qs[0]);
  Controlled Ry([qs[0]], (-1.570796326794897, qs[1]));
  X(qs[1]);
  Controlled Ry([qs[1]], (-1.910633236249018, qs[0]));
  X(qs[1]);
  Controlled Rz([qs[0]], (-4.712388980384691, qs[1]));
  Controlled Ry([qs[0]], (-1.570796326794897, qs[1]));
  Controlled Rz([qs[0]], (-1.570796326794896, qs[1]));
  Controlled Rz([qs[1]], (-1.570796326794897, qs[0]));
  Controlled Ry([qs[1]], (-3.141592653589793, qs[0]));
  Controlled Rz([qs[1]], (1.570796326794897, qs[0]));
}

And this is decomposition of $$$A$$$ into two-level unitaries (result of steps 1 and 2):

$$$A = \begin{pmatrix} 1 & 0 & 0 & 0 \\0& 1 & 0& 0 \\ 0 & 0 & 0 & -i \\ 0 & 0 & -i & 0 \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & -\frac{\sqrt{2}}{2} & 0 & -\frac{\sqrt{2}}{2}i \\ 0 & 0 & 1 & 0 \\ 0 & -\frac{\sqrt{2}}{2}i & 0 & -\frac{\sqrt{2}}{2} \end{pmatrix} \begin{pmatrix} \sqrt{\frac{1}{3}} & \sqrt{\frac{2}{3}} & 0 & 0 \\ -\sqrt{\frac{2}{3}} & \sqrt{\frac{1}{3}} & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & \frac{\sqrt{2}}{2} & 0 & \frac{\sqrt{2}}{2} \\ 0 & 0 & 1 & 0 \\ 0 & -\frac{\sqrt{2}}{2} & 0 & \frac{\sqrt{2}}{2} \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{pmatrix} $$$

Notes

This algorithm in theory implements an arbitrary unitary matrix, but in practice it works only when the number of qubits is small (less than 10). Indeed, for $$$N$$$ qubits, the unitary matrix has $$$4^N$$$ entries, so the algorithm can't be faster than $$$O(4^N)$$$. In fact, it produces $$$O(4^N)$$$ gates. And one can show that one can't do better (because a unitary matrix has $$$O(4^N)$$$ degrees of freedom, and every rotation gate adds just one degree of freedom).
Version of this algorithm can be found in any textbook on Quantum Computing, for example in chapter 4.5 of "Quantum Computing and Quantum Information" by Nielsen & Chuang. However, such textbooks usually present it as the theorem that any unitary gate can be decomposed into a CNOT gate and single-qubit gates, and it's hard to go from proof of that theorem to a practical algorithm. What I needed is a practical algorithm which converts unitary matrix to Q# gates, and at the moment when I needed it I didn't find any implementation for such algorithm online, so I wrote my own implementation.
Definitions of $$$R_y$$$ and $$$R_z$$$ gates in this post (and the paper) differ from those in Q# by argument sign.
When I use notation $$$i[j]$$$ in step 3, this means we view index $$$i$$$ as binary string and take $$$j$$$-th bit of it.

References

Algorithm for two-level decomposition is based on this paper.
I implemented this algorithm as a Python library, which can be found in this Github repository. This library takes unitary matrix and produces Q# code implementing that matrix. Also in the repository you can find Jupyter Notebook with more detailed examples.
I described the algorithm in this paper. This blog post is based on that paper, but omits some details, so please refer to the paper for more rigorous derivation.

I hope this blog post was useful to you. You’re welcome to use the library in the next Q# coding contest!