Saving statevector on more than one location in a quantum circuit in Qiskit

Question

So, I'm fairly new to Qiskit, and I've been playing around and following the tutorials from the Qiskit textbook. However, there is one thing I fail to understand/implement: for a quantum circuit with multiple gates (e.g. 2-qubit circle, with a simple Hadamard acting on qubit 0 and, let's say, CNOT acting on qubit 1 (controlled by qubit 0)) is it possible somehow to save statevectors more than once, in order to see intermediate state of the system?

Below you can see the code idea:

q1 = QuantumCircuit(2)
q1.save_statevector() # Save initial state
q1.h(0)
q1.save_statevector() # Save state after Hadamard
q1.cx(0, 1)
q1.save_statevector() # Save state after CNOT (also a final state)
job = execute(q1, backend=Aer.get_backend('aer_simulator'), shots=1024)
statevectors = job.result().get_statevector()

However, if I were to try and run this, an error occurs upon reaching execute command. If anyone can provide any insight on this, I would be very grateful.

Answer 1

You can also obtain the states at any point during circuit construction using Statevector, the class from Qiskit's quantum_info module as follows.

First, import the Statevector class,

from qiskit.quantum_info import Statevector

And for your example, the code below will produce all the intermediate states that you want.

qc = QuantumCircuit(2)
st0 = Statevector.from_instruction(qc) 
qc.h(0)
st1 = Statevector.from_instruction(qc)
qc.cnot(0, 1)
st2 = Statevector.from_instruction(qc)

print(st0)
print(st1)
print(st2)


Statevector([1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
        dims=(2, 2))
Statevector([0.70710678+0.j, 0.70710678+0.j, 0.             +0.j,
         0.        +0.j],
        dims=(2, 2))
Statevector([0.70710678+0.j, 0.        +0.j, 0.        +0.j,
         0.70710678+0.j],
        dims=(2, 2))

Answer 2

This is because you have multiple save_statevector() instruction in one circuit.

If you only have one instruction of save_statevector() then this is fine. That is:

from qiskit import transpile 
q1 = QuantumCircuit(2)
q1.h(0)
q1.cx(0, 1)
q1.save_statevector() # Save state after CNOT (also a final state)
backend=Aer.get_backend('aer_simulator') 
q1 = transpile(q1, backend)
result = backend.run(q1).result()
statevector = result.get_statevector()

which will return your statevector here as:

array([0.70710678+0.j, 0.        +0.j, 0.        +0.j, 0.70710678+0.j])

---

If you have more than one place in the circuit that you want to save the statevector of, by using multiple save_statevector() methods in the circuit, then what you can do is copy the circuit at each of these step to another circuit. At the end, you can execute the simulation for each of these copied circuit independently and extract out your statevector through get_statevector() method.

For example: Using your original problem, we can do something like

from qiskit import transpile 
import copy 
main_qc = QuantumCircuit(2)
q1 = copy.deepcopy(main_qc)
q1.save_statevector() # Save initial state

main_qc.h(0)
q2 =  copy.deepcopy(main_qc)
q2.save_statevector() # Save state after Hadamard

main_qc.cx(0, 1)
q3 =  copy.deepcopy(main_qc)
q3.save_statevector() # Save state after CNOT (also a final state)

backend=Aer.get_backend('aer_simulator') 

statevectors_circuits = [ q1, q2, q3]
transpiled_circuits = [transpile(circuit, backend) for circuit in circuits] 
results =[ backend.run(transpile_circuit).result() for transpile_circuit in transpiled_circuits] 

statevectors = [results[i].get_statevector(circuits[i]) for i in range(len(statevectors_circuits))  ]

Now, if you print out the statevectors, you will see:

[array([1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]),
 array([0.70710678+0.j, 0.70710678+0.j, 0.        +0.j, 0.        +0.j]),
 array([0.70710678+0.j, 0.        +0.j, 0.        +0.j, 0.70710678+0.j])]

These are the 3 states that you were trying to save.

Answer 3

It took me some time, but I found a more direct approach, alongside my initial idea. It turns out that statevector can indeed be saved multiple times for a given circuit, where one must provide a unique label each time state is saved. The same reasoning should also be valid for, e.g. density_matrix and unitary type of constructs (availability of specific construct depends on the selected backend). So, for example, consider the following code which stores initial and final statevector and density_matrix, respectively:

# Select simulation backend:
aer_sim = Aer.get_backend('aer_simulator')
# Set number of simulation iterations:
Nshot = 1e5

# Allocate empty circuit with 3 qubits:
qc = QuantumCircuit(3)
# Save initial state:
qc.save_statevector(label=f'init_psi')
# Save initial density matrix, too:
qc.save_density_matrix(label=f'init_rho')
# Set qubits to the GHZ 3-state:
qc.h(0)
qc.cx(0, 1)
qc.cx(0, 2)
# Save final state:
qc.save_statevector(label=f'final_psi')
# Save final density matrix, too:
qc.save_density_matrix(label=f'final_rho')

# Transpile the circuit:
qc_tp = transpile(qc, aer_sim)
# Run the simulation:
job = aer_sim.run(qc_tp, shots=Nshots)
# Get results:
res = job.result()
# Print initial state:
print(res.data(0)['init_psi'])
# Print final state:
print(res.data(0)['final_psi'])
# Print initial density matrix:
print(res.data(0)['init_rho'])
# Print final density matrix:
print(res.data(0)['final_rho'])

In terms of absolute bookkeeping where one might be interested in the incremental changes of the system (unitary, statevector, and density_matrix), an alternative approach would be to follow the suggestion of @HwajungKang and import Statevector, DensityMatrix, and/or Operator objects from ```qiskit.quantum_info``` which are not bounded by the selected simulator backend.