How to find an equivalent circuit without ancilla qubits?

Question

$\newcommand{\ket}[1]{|#1\rangle}$

I have the following quantum circuit:

(The inner qubits are both initialized to $|i\rangle$. $U$ is a arbitrary quantum gate.)

But I am only interested in the outcome of the qubits $q_0$ and $q_3$. Is there a way to reduce this circuit to only 2 qubits?

My first idea was to replace all 2-qubit gates by single qubit operations and once "isolated" remove the inner 2 qubits. The paper Constructing a virtual two-qubit gate by sampling single-qubit operations describes a strategy to decompose a two-qubit gate to a sequence of single-qubit operations (with sampling overhead), which leads me to believe that this circuit can be reduced to only 2 qubits.

Does my apprach make sense? Is there a general way to get rid of ancilla qubits?

---

Edit

Since my question lacked clarity, I now provide my calculation to be more specific in what I want.

For my own convenience I reordered the qubits according to the circuit below:

The $U$ gate is an arbitrary unitary $U=\left[\begin{matrix}u_{0} & u_{1}\\u_{2} & u_{3}\end{matrix}\right]$. I represent the combined state of $q_3$ and $q_0$ as $a\ket{00} + b\ket{01} + c\ket{10} + d\ket{11} =\left[\begin{matrix}a\\b\\c\\d\end{matrix}\right]$, which results in the overall initial state $\ket{\psi}=\ket{ii}\otimes\left[\begin{matrix}a\\b\\c\\d\end{matrix}\right]$.

I use sympy to calculate the circuit $C_{2}{\left(Z_{1}\right)} H_{2} C_{2}{\left(U_{0}\right)} \text{CNOT}_{3,2} \text{CNOT}_{1,3}\ket{\psi}$

from sympy import *
from sympy.physics.quantum.tensorproduct import TensorProduct
from sympy.physics.quantum.gate import Z, H
from qiskit import QuantumCircuit
from qiskit.quantum_info import Operator

# define variables
a, b, c, d = var('a b c d', complex=True)
u0, u1, u2, u3 = var('u0 u1 u2 u3', complex=True)
u = Matrix([[u0, u1], [u2, u3]])
cx = Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
cz = Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]])

# define state & its density matrix
state = TensorProduct(Matrix([1, I]) / sqrt(2), Matrix([1, I]) / sqrt(2), Matrix([a, b, c, d]))
state_d = state * state.H
# CNOT(q_3, i_0)
op_qc = QuantumCircuit(4)
op_qc.cx(1, 3)
op = Matrix(Operator(op_qc))
t = simplify(op * state_d * op.H)
# CNOT(i_0, i_1)
op = TensorProduct(cx, eye(4))
t = simplify(op * t * op.H)
# CU(i_1, q_0)
op = TensorProduct(eye(2), eye(2), eye(2), eye(2)) / 2 \
    + TensorProduct(eye(2), Z().get_target_matrix(), eye(2), eye(2)) / 2 \
    + TensorProduct(eye(2), eye(2), eye(2), u) / 2 \
    - TensorProduct(eye(2), Z().get_target_matrix(), eye(2), u) / 2
t = simplify(op * t * op.H)
# H(i_1)
op = TensorProduct(eye(2), H().get_target_matrix(), eye(4))
t = simplify(op * t * op.H)
# CZ(i_1, q_3)
op = TensorProduct(eye(2), cz, eye(2))
t = simplify(op * t * op.H)
# trace out qubits i_0 and i_1
r = TensorProduct(Matrix([1, 0, 0, 0]).H, eye(4)) * t * TensorProduct(Matrix([1, 0, 0, 0]), eye(4)) + \
    TensorProduct(Matrix([0, 1, 0, 0]).H, eye(4)) * t * TensorProduct(Matrix([0, 1, 0, 0]), eye(4)) + \
    TensorProduct(Matrix([0, 0, 1, 0]).H, eye(4)) * t * TensorProduct(Matrix([0, 0, 1, 0]), eye(4)) + \
    TensorProduct(Matrix([0, 0, 0, 1]).H, eye(4)) * t * TensorProduct(Matrix([0, 0, 0, 1]), eye(4))
r = simplify(r)

The resulting density matrix for $q_3$ and $q_0$ is:

$\left[\begin{matrix}0.5 a u_{0} \overline{a} \overline{u_{0}} + 0.5 a u_{0} \overline{b} \overline{u_{1}} + 0.5 a \overline{a} + 0.5 b u_{1} \overline{a} \overline{u_{0}} + 0.5 b u_{1} \overline{b} \overline{u_{1}} & 0.5 a u_{0} \overline{a} \overline{u_{2}} + 0.5 a u_{0} \overline{b} \overline{u_{3}} + 0.5 a \overline{b} + 0.5 b u_{1} \overline{a} \overline{u_{2}} + 0.5 b u_{1} \overline{b} \overline{u_{3}} & 0.5 a u_{0} \overline{c} - 0.5 a \overline{c} \overline{u_{0}} - 0.5 a \overline{d} \overline{u_{1}} + 0.5 b u_{1} \overline{c} & 0.5 a u_{0} \overline{d} - 0.5 a \overline{c} \overline{u_{2}} - 0.5 a \overline{d} \overline{u_{3}} + 0.5 b u_{1} \overline{d}\\0.5 a u_{2} \overline{a} \overline{u_{0}} + 0.5 a u_{2} \overline{b} \overline{u_{1}} + 0.5 b u_{3} \overline{a} \overline{u_{0}} + 0.5 b u_{3} \overline{b} \overline{u_{1}} + 0.5 b \overline{a} & 0.5 a u_{2} \overline{a} \overline{u_{2}} + 0.5 a u_{2} \overline{b} \overline{u_{3}} + 0.5 b u_{3} \overline{a} \overline{u_{2}} + 0.5 b u_{3} \overline{b} \overline{u_{3}} + 0.5 b \overline{b} & 0.5 a u_{2} \overline{c} + 0.5 b u_{3} \overline{c} - 0.5 b \overline{c} \overline{u_{0}} - 0.5 b \overline{d} \overline{u_{1}} & 0.5 a u_{2} \overline{d} + 0.5 b u_{3} \overline{d} - 0.5 b \overline{c} \overline{u_{2}} - 0.5 b \overline{d} \overline{u_{3}}\\- 0.5 c u_{0} \overline{a} + 0.5 c \overline{a} \overline{u_{0}} + 0.5 c \overline{b} \overline{u_{1}} - 0.5 d u_{1} \overline{a} & - 0.5 c u_{0} \overline{b} + 0.5 c \overline{a} \overline{u_{2}} + 0.5 c \overline{b} \overline{u_{3}} - 0.5 d u_{1} \overline{b} & 0.5 c u_{0} \overline{c} \overline{u_{0}} + 0.5 c u_{0} \overline{d} \overline{u_{1}} + 0.5 c \overline{c} + 0.5 d u_{1} \overline{c} \overline{u_{0}} + 0.5 d u_{1} \overline{d} \overline{u_{1}} & 0.5 c u_{0} \overline{c} \overline{u_{2}} + 0.5 c u_{0} \overline{d} \overline{u_{3}} + 0.5 c \overline{d} + 0.5 d u_{1} \overline{c} \overline{u_{2}} + 0.5 d u_{1} \overline{d} \overline{u_{3}}\\- 0.5 c u_{2} \overline{a} - 0.5 d u_{3} \overline{a} + 0.5 d \overline{a} \overline{u_{0}} + 0.5 d \overline{b} \overline{u_{1}} & - 0.5 c u_{2} \overline{b} - 0.5 d u_{3} \overline{b} + 0.5 d \overline{a} \overline{u_{2}} + 0.5 d \overline{b} \overline{u_{3}} & 0.5 c u_{2} \overline{c} \overline{u_{0}} + 0.5 c u_{2} \overline{d} \overline{u_{1}} + 0.5 d u_{3} \overline{c} \overline{u_{0}} + 0.5 d u_{3} \overline{d} \overline{u_{1}} + 0.5 d \overline{c} & 0.5 c u_{2} \overline{c} \overline{u_{2}} + 0.5 c u_{2} \overline{d} \overline{u_{3}} + 0.5 d u_{3} \overline{c} \overline{u_{2}} + 0.5 d u_{3} \overline{d} \overline{u_{3}} + 0.5 d \overline{d}\end{matrix}\right]$

---

And here is my calculation for the circuit suggested in the comments:

from sympy import *
from sympy.physics.quantum.tensorproduct import TensorProduct
from sympy.physics.quantum.gate import Z

a, b, c, d = var('a b c d', complex=True)
u0, u1, u2, u3 = var('u0 u1 u2 u3', complex=True)
cu = Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, u0, u1], [0, 0, u2, u3]])

# define state & its density matrix
state = Matrix([[a], [b], [c], [d]])
state_d = state * state.H
# CU(q_3, q_0)
t = simplify(cu * state_d * cu.H)
# Z(q_3)
op = TensorProduct(Z().get_target_matrix(), eye(2))
r = simplify(op * t * op)

Result:

$\left[\begin{matrix}a \overline{a} & a \overline{b} & - a \left(\overline{c} \overline{u_{0}} + \overline{d} \overline{u_{1}}\right) & - a \left(\overline{c} \overline{u_{2}} + \overline{d} \overline{u_{3}}\right)\\b \overline{a} & b \overline{b} & - b \left(\overline{c} \overline{u_{0}} + \overline{d} \overline{u_{1}}\right) & - b \left(\overline{c} \overline{u_{2}} + \overline{d} \overline{u_{3}}\right)\\- \left(c u_{0} + d u_{1}\right) \overline{a} & - \left(c u_{0} + d u_{1}\right) \overline{b} & \left(c u_{0} + d u_{1}\right) \left(\overline{c} \overline{u_{0}} + \overline{d} \overline{u_{1}}\right) & \left(c u_{0} + d u_{1}\right) \left(\overline{c} \overline{u_{2}} + \overline{d} \overline{u_{3}}\right)\\- \left(c u_{2} + d u_{3}\right) \overline{a} & - \left(c u_{2} + d u_{3}\right) \overline{b} & \left(c u_{2} + d u_{3}\right) \left(\overline{c} \overline{u_{0}} + \overline{d} \overline{u_{1}}\right) & \left(c u_{2} + d u_{3}\right) \left(\overline{c} \overline{u_{2}} + \overline{d} \overline{u_{3}}\right)\end{matrix}\right]$

Answer 1

As drawn, it's impossible to remove qubits 2 and 3 from your circuit.

The second qubit has to interact with the first qubit somehow, because it ends up flipped conditioned on the first qubit.

The fourth qubit must interact with the second and third qubits somehow because $U$ is applied exactly when an odd number of qubits 1,2,3 were ON at the start.

The only way it would be possible to simplify the circuit would be if the second and third qubits had known values at the start or end. Alternatively, knowing $U$ was self-inverse would be helpful for breaking up the controls.

Here's an equivalent form of the circuit that makes these relationships clearer (the "Zpar" controls are parity controls; they are satisfied as a group when an odd number are ON):

Answer 2

The density matrix for $q_0$ and $q_3$ as written in the question can be computed from the following circuit:

This circuit only requires 1 ancilla qubit, which needs to be traced out.