How to decompose a multi-target controlled gate?

Question

I'm trying to replicate with qiskit the results of this paper in which basically they implement a quantum version of the Principal Component Analysis applying Quantum Phase Estimation algorithm in order to find eigenvalues and eigenvectors from a covariance matrix.

Everything works fine for the 2x2 covariance matrix case but when they extend the circuit to the 4x4 case it's not very clear how they modify the phase kickback part of the quantum circuit. In the picture that shows the general circuit for QPE, some multi-target controlled-U gates are shown:

In the 2x2 case this can be made using simple cU3 gates repeated several times as shown in this answer, where the parameters $\theta, \phi, \lambda$ that are needed for the cU3 gate are calculated using the qiskit function euler_angles_1q directly from the $e^{2\pi i\rho}$ 2x2 matrix where $\rho$ is the normalized covariance matrix:

$\rho=\frac{\sigma}{tr(\sigma)}$

($\sigma$ is the 2x2 covariance matrix for which we want to find the eigenvalues and the eigenvectors)

But for the 4x4 case they provide a circuit that does the phase kickback with several U3 gates

My question is:

1. Since I am a beginner I don't understand how a controlled multi-target gate is decomposed in this chain of gates. Could you please explain me the idea behind this decomposition or give me some resources to start from in order to understand it?
2. Is it possible to decompose a multi-target cU3 gate using only 2 single target cU3 gates?

Answer 1

An example of constructing (with help of Qiskit) a controlled version of some simple 4x4 unitary matrix:

$$ U = \begin{pmatrix} \mathrm{e}^{i g_1} & 0 & 0 & 0 \\ 0 & \mathrm{e}^{i g_2} & 0 & 0 \\ 0 & 0 & \mathrm{e}^{i g_3} & 0 \\ 0 & 0 & 0 & \mathrm{e}^{i g_4} \\ \end{pmatrix} $$

where $g_s$ are some given constants. It can be shown that this corresponds to the following circuit:

from qiskit import *

g_1 = 1.2
g_2 = 1.5
g_3 = 0.7
g_4 = -0.9

q = QuantumRegister(2, "q")
circuit = QuantumCircuit(q)

circuit.cu1(g_4, q[0], q[1])
circuit.x(q[1])
circuit.cu1(g_3, q[0], q[1])
circuit.x(q[0])
circuit.cu1(g_1, q[0], q[1])
circuit.x(q[1])
circuit.cu1(g_2, q[0], q[1])
circuit.x(q[0])

to obtain the controlled version of the circuit we will use Qiskit's get_controlled_circuit() method:

from qiskit.aqua.utils.controlled_circuit import get_controlled_circuit

control_qubit = QuantumRegister(1, "a")
controlled_circuit = QuantumCircuit(q, control_qubit)

controlled_circuit = get_controlled_circuit(circuit, control_qubit[0])
print(controlled_circuit.qasm())

How I understand the steps of the get_controlled_circuit() method: 1) it takes your given circuit 2) changes your gates to Qiskit's basis gates (u1, u2, u3, cx). 3) for each given basis gate it modifies and replaces it with the controlled version of it (it has a "dictionary of methods" that implements the controlled circuits for all basis gates).

So, to sum up. Firstly one needs to decompose the given $U$ operator into gates, then for each gate construct the controlled version of it. In other words, if one has a circuit gate_1 gate_2 gate_3, the controlled version of it will look like this controlled_gate_1 controlled_gate_2 controlled_gate_3. If the control qubit is in the $| 1 \rangle$ state all gates will be executed, and if it is in the $| 0 \rangle$ state nothing will happen (as is expected by definition of the control operations). And if you want to implement a controlled version of $U^n$ operator, you can just apply $n$ controlled $U$ circuits.

Final Notes: To obtain a circuit from a given unitary matrix is not a trivial problem and I don't know a general solution for that. Usually, in QPE algorithm implementations, one is given a Hamiltonian $H$ operator expressed by a sum of Pauli product terms and the problem is to simulate $e^{iHt}$ unitary operator to estimate the phase $\varphi$, from $e^{iHt} | \psi \rangle = e^{iEt} | \psi \rangle = e^{i 2 \pi \varphi} | \psi \rangle$, where $| \psi \rangle$ is some eigenstate of $H$. The $H$ is given in this form:

$$H = P_1 + P_2 + P_3 + P_3$$

where $P_s$ are the Pauli terms (like $X \otimes Y \otimes Y$ or $Z \otimes Y$). Simulation of the $e^{i P t}$ for a given $P$ Pauli term can be done as described in J.D. Whitfield et al [1] paper. If we have all circuits for Pauli terms we will be able to simulate the whole Hamiltonian with the use of Trotterization techniques (more info can be found in [1] paper).