Operators are forward-propagated through depth >1 slices

[BryceFuller]

Environment

• OBP=0.1.0:
• Python=3.11.8:
• OSX Sonoma 14.5:

What is happening and why is it wrong?

What's wrong

When a list of slices are passed to the backpropagate method, they are processed in reverse order (last to first). However, if the slices are depth >1 then the instructions within the circuit are currently not being processed in reverse order.

Example 1

Say we have two slices, each made up of two instructions.
`slices' ~= [($U_1$, $U_2$), ($U_3$, $U_4$)]

When we backpropagate an observable $O$, we expect OBP to approximate
$O' = (U^\dagger_1 U^\dagger_2) (U^\dagger_3 U^\dagger_4)\ O\ (U_4 U_3) (U_2 U_1)$

But because the instructions within each slice are not processed in reverse order what we actually compute is
$O' = (U^\dagger_2 U^\dagger_1 U^\dagger_4 U^\dagger_3)\ O\ (U_3 U_4 U_1 U_2)$

Example 2

In the extreme case, which is reproduced below in code, circuits passed in as a single slice are forward-propagated entirely.

`slices' ~= [($U_1$, $U_2$, $U_3$, $U_4$)]

The current code will compute
$O' = (U^\dagger_4 U^\dagger_3 U^\dagger_2 U^\dagger_1)\ O\ (U_1 U_2 U_3 U_4)$

How can we reproduce the issue?

import numpy as np
from copy import deepcopy
from qiskit import QuantumCircuit
from qiskit.compiler import transpile
from qiskit.quantum_info import SparsePauliOp, Operator
from qiskit_aer import AerSimulator
from qiskit_addon_obp import backpropagate
from qiskit_addon_obp.utils.simplify import OperatorBudget
from qiskit_addon_utils.slicing import slice_by_gate_types

Define a circuit with an explicit unitary matrix

unitary_mat = np.array(
    [
        [1/np.sqrt(2), 1/2, 1/np.sqrt(8), 1/np.sqrt(8)],
        [1/np.sqrt(2), -1/2, -1/np.sqrt(8),-1/np.sqrt(8)],
        [0, 1/np.sqrt(2), -1/2, -1/2],
        [0, 0, 1/np.sqrt(2), -1/np.sqrt(2)]
    ])

circ = QuantumCircuit(2)
circ.unitary(Operator(unitary_mat), [0,1], label='split')
circ = transpile(circ, basis_gates=['rz', 'cx', 'cz', 'h'], optimization_level=3)

Define an observable $O$

O = SparsePauliOp.from_list([('IX',1),('XI',1),('XX',1),('IZ',1),('ZI',1),('ZZ',1)])
# Prepare the matrix representation of this observable
O_mat = O.to_matrix()

Compare $U^\dagger O U$ Computed with OBP vs computed directly with matrices

slices = slice_by_gate_types(circ)
op_budget = OperatorBudget()

# O' = U†•O•U computed via backprop
O_prime, remaining_slices, metadata = backpropagate(O, slices, operator_budget=op_budget)
O_prime_mat = O_prime.to_matrix()

# O' = U†•O•U computed directly via matrices
O_prime_mat_direct = unitary_mat.conj().T @ O_mat @ unitary_mat

print("These two matrices are equal, which is correct")
np.round(O_prime_mat_direct - O_prime_mat, 3)

>>>
These two matrices are equal, which is correct

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

Perform the same calculation but for OBP pass in a single slice with the whole circuit

op_budget = OperatorBudget()

# O' = U†•O•U computed via backprop
O_prime, remaining_slices, metadata = backpropagate(O, [circ], operator_budget=op_budget)
O_prime_mat = O_prime.to_matrix()

# O' = U†•O•U computed directly via matrices
O_prime_mat_direct = unitary_mat.conj().T @ O_mat @ unitary_mat

print("These two matrices are now not equal, which is Wrong")
np.round(O_prime_mat_direct - O_prime_mat, 3)
>>>
array([[-1.664-0.j,  1.371+0.j,  1.854-0.j, -0.707+0.j],
       [ 1.371+0.j, -0.75 +0.j,  0.707+0.j, -0.146+0.j],
       [ 1.854-0.j,  0.707+0.j,  0.457+0.j,  0.25 -0.j],
       [-0.707+0.j, -0.146+0.j,  0.25 -0.j,  1.957-0.j]])

If instead we compare the output of OBP with the forward propagated operator $O* = U O U^\dagger$, we see agreement.

# O* = U•O•U† computed directly via matrices
O_star_mat_direct = unitary_mat @ O_mat @ unitary_mat.conj().T
np.round(O_star_mat_direct - O_prime_mat, 3)
>>>
array([[ 0.-0.j,  0.+0.j, -0.-0.j, -0.+0.j],
       [ 0.+0.j,  0.+0.j,  0.+0.j,  0.+0.j],
       [-0.-0.j,  0.+0.j, -0.+0.j, -0.-0.j],
       [-0.+0.j,  0.+0.j, -0.-0.j, -0.-0.j]])

Traceback

No response

Any suggestions?

I think this can be addressed by changing Line 167 of backpropagation.py from

for op_idx, op_node in enumerate(circuit_to_dag(slice_).topological_op_nodes()):

to

for op_idx, op_node in list(enumerate(circuit_to_dag(slice_).topological_op_nodes()))[::-1]:

I realize that this casts the generator to a sequence in order to reverse it, so there likely is a more elegant solution.