Removing ancilla qubits in Qiskit circuit

Question

I want to remove the ancilla qubits from my quantum circuit in Qiskit. My MWE below contains four qubits: psi0, psi1, anc0, anc1. I want to keep psi0 and psi1 qubits and remove ancillas anc0, anc1. Is there a way to do this?

from qiskit import QuantumRegister, QuantumCircuit, AncillaRegister, transpile
from qiskit_aer import StatevectorSimulator
import numpy as np

reg1 = QuantumRegister(2, 'psi')
reg2 = AncillaRegister(2, 'anc')
circ = QuantumCircuit(reg1, reg2)
circ.h(reg1[0])
circ.h(reg1[1])
circ.prepare_state(1/np.sqrt(2)*np.array([1,0,0,1]),reg2)
circ.cry(np.pi/2.2, control_qubit=reg2[1], target_qubit=reg1[1])
circ.cry(np.pi/1.5, control_qubit=reg2[0], target_qubit=reg1[0])
circ.prepare_state(1/np.sqrt(2)*np.array([1,0,0,1]),reg2).inverse()

simulator = StatevectorSimulator()
job = transpile(circ, simulator)  # execute circuit on qasm simulator
result = simulator.run(job,shots=128).result()  # results from the job
out_state = result.get_statevector()
print(out_state)

circ.draw('mpl')

I'm aware of a function that removes idle wires, but it loses the register's info. Is there a different function that specifically removes ancilla qubits (and keeps the register's info)?

Edit: Provided clearer circuit example. To clarify: Yes it will be ok to keep the ancilla qubits until the end of the circuit. I just need psi qubits in the next part of the calculation. I want to do a statevector simulation of psi, not anc, at the end. Ideally I don't want to initialize a separate circuit using the output density matrix -- I prefer to use one self-contained circuit.

Edit 2: The comments below suggested simulating the full statevector of all qubits, then taking the trace to remove ancillas. But the partial_trace function can only be applied to density matrices (see this post).

• What is an equivalent method to trace over statevectors?

• Is it possible to automatically determine all the ancillas and trace them out?

• If there are many ancillas, will it be computationally expensive to simulate the full statevector rather than removing the ancillas before simulation?

Answer 1

As I mentioned in the comments, if your psi qubits are getting entangled with the anc qubits, there is no way to recover a Statevector just for psi because the state will not be pure.

Instead, what you can extract is the density matrix of psi by tracing out the anc register. That being said, if the anc qubits are properly untentangled from psi you can always convert that DensityMatrix object back into a Statevector object.

Here are two examples, one where you have an unentangled state between psi and anc (So you can recover the Statevector for psi), and one where they are entangled.

Example 1

from qiskit import QuantumRegister, QuantumCircuit, AncillaRegister
from qiskit_aer import AerSimulator

reg1 = QuantumRegister(1, 'psi')
reg2 = AncillaRegister(1, 'anc')
circ = QuantumCircuit(reg1, reg2)
circ.h(reg1)
circ.x(reg2)
circ.cz(reg1,reg2)
circ.h(reg1)
circ.x(reg2)
circ.save_density_matrix(reg1)

ρ = AerSimulator().run(circ).result().data()['density_matrix']
ψ = ρ.to_statevector()
print(ψ)

This will print state Statevector([0, 1],dims=(2,)) because the state is separable.

Example 2

Change the cz gate above for a cx gate and then this line of code ψ = ρ.to_statevector() will give an error because the density matrix extracted from the simulation does not correspond to one of a pure state.