Decomposition of arbitrary 2 qubit operator

Question

As you know, universal quantum computing is the ability to construct a circuit from a finite set of operations that can approximate to arbitrary accuracy any unitary operation.

There also exist some results proving that exact decompositions of particular unitary operations can be found. For instance, a method was provided here (see section V) to construct a quantum circuit for a general two-qubit gate, based on the method given by Kraus and Cirac to decompose any desired two-qubit gate.

I want to understand this method such that I learn how to find the quantum circuit implementing any $4\times4$ unitary matrix. To help with that, I devised the following challenge: You are given an unitary matrix, $U$, over two qubits and asked to follow either the previously mentioned methods or another method you know about to come up with the quantum circuit implementing this unitarity. I know what one circuit representation is because I built it myself to obtain this matrix. Our goal should be to find the decomposition with the fewer number of gates, however, other decomposition methods will also be useful to learn. Here's the matrix

$$ U\left(x\right)=\left(\begin{array}{cccc} 0 & \left(-\frac{1}{2}+\frac{i}{2}\right)e^{ix/2} & 0 & \left(\frac{1}{2}+\frac{i}{2}\right)e^{ix/2}\\ 0 & \frac{ie^{-ix/2}}{\sqrt{2}} & 0 & -\frac{e^{-ix/2}}{\sqrt{2}}\\ \left(-\frac{1}{2}-\frac{i}{2}\right)e^{-ix/2} & 0 & \left(\frac{1}{2}-\frac{i}{2}\right)e^{-ix/2} & 0\\ -\frac{e^{ix/2}}{\sqrt{2}} & 0 & \frac{ie^{ix/2}}{\sqrt{2}} & 0 \end{array}\right) $$

If you are aware of any numerical methods that can achieve this task I also encourage you to post the solution you find and the steps to find it. E.g., as answered at another question on this site, qubiter has a module which apparently can decompose any n-qubit unitary into cnots and single qubit rotations using the cosine-sine decomposition, and this is probably not the only numerical method looking into this problem.

I believe an exhaustive exploration of these methods will be a helpful reference to a big audience, but let us start simple: can you give a circuit decomposition for $U\left(x\right)$?

Universal gate set:

You may consider the following set of basic gates: hadamard ($H$), phase ($S$), $\pi/8$ rotation ($T$), $cX$, $S^{\dagger}$, $T^{\dagger}$, and rotations

$$ R_{x}\left(\theta\right)=\left(\begin{array}{cc} \cos\left(\frac{\theta}{2}\right) & -i\sin\left(\frac{\theta}{2}\right)\\ -i\sin\left(\frac{\theta}{2}\right) & \cos\left(\frac{\theta}{2}\right) \end{array}\right),\quad R_{y}\left(\theta\right)=\left(\begin{array}{cc} \cos\left(\frac{\theta}{2}\right) & -\sin\left(\frac{\theta}{2}\right)\\ \sin\left(\frac{\theta}{2}\right) & \cos\left(\frac{\theta}{2}\right) \end{array}\right),\quad R_{z}\left(\theta\right)=\left(\begin{array}{cc} e^{-\frac{i\theta}{2}} & 0\\ 0 & e^{\frac{i\theta}{2}} \end{array}\right) $$

With an eye on a real implementation, you may also consider IBM QX universal gate set made of a $cX$ together with one-qubit rotational and phase gates

$$ \begin{split}V_{1}(\lambda) & =\begin{pmatrix}1 & 0\\ 0 & e^{i\lambda} \end{pmatrix}\\ V_{2}(\phi,\lambda) & =\frac{1}{\sqrt{2}}\begin{pmatrix}1 & -e^{i\lambda}\\ e^{i\phi} & e^{i(\lambda+\phi)} \end{pmatrix}\\ V_{3}(\theta,\phi,\lambda) & =\begin{pmatrix}\cos\left(\frac{\theta}{2}\right) & -e^{i\lambda}\sin\left(\frac{\theta}{2}\right)\\ e^{i\phi}\sin\left(\frac{\theta}{2}\right) & e^{i(\lambda+\phi)}\cos\left(\frac{\theta}{2}\right) \end{pmatrix} \end{split} $$

Actually, $V_{1}$ and $V_{2}$ are just special cases of $V_{3}$, hence IBM's universal set can be reduced to just two gates.

Example:

Suppose you wanted to find the circuit implementation for an unitarity $U\left(x,y\right)$ given by the following matrix

$$ \left(\begin{array}{cccc} e^{-\frac{1}{2}i\pi y}\cos\frac{x}{2} & 0 & 0 & -e^{-\frac{1}{2}i\pi y}\sin\frac{x}{2}\\ 0 & -ie^{\frac{i\pi y}{2}}\sin\frac{x}{2} & -ie^{\frac{i\pi y}{2}}\cos\frac{x}{2} & 0\\ \sqrt[4]{-1}e^{-\frac{1}{2}i\pi y}\sin\frac{x}{2} & 0 & 0 & \sqrt[4]{-1}e^{-\frac{1}{2}i\pi y}\cos\frac{x}{2}\\ 0 & -(-1)^{3/4}e^{\frac{i\pi y}{2}}\cos\frac{x}{2} & (-1)^{3/4}e^{\frac{i\pi y}{2}}\sin\frac{x}{2} & 0 \end{array}\right) $$

You can check $U\left(x,y\right)$ can be decomposed into the following product of operations

$$ U\left(x,y\right)=\left(T\otimes S^{\dagger}\right)\cdot CNOT_{01}\cdot\left(R_{y}\left(x\right)\otimes R_{z}\left(y\pi\right)\right)\cdot CNOT_{10} $$

and the circuit representation is

Answer 1

The circuit that I came up with is: up to a global phase factor of $-e^{i x/2 + i \pi/4}$. The phase gate "x" that I used is essentially $R_z(x/2)=V_1(x)$ up to a global phase factor, so this is directly applicable to the first gates set. There's some possibility to move around some of the single-qubit gates, but I think this option minimises the depth of the circuit, which may be a more relevant number than the absolute number of gates, as it's more related to the time a circuit would take to implement (and hence how susceptible to noise it might be). Of course, I don't pretend that all gates would take the same length of time, and that's a completely different set of rules that one would be optimising over.

This same circuit can be rewritten using the other gate set

To come up with this decomposition, I didn't use any particular method beyond the standard "I stared at it for a long time" technique. I then wrote down a circuit which was always going to be sub-optimal (including a swap gate, and a two-qubit diagonal gate with arbitrary phases on the diagonal), but then rearranged the circuit to bring some elements together so that they cancelled.