10 hybrid schrodinger heisenberg

Hybrid Schrodinger-Heisenberg expectation values¶

For a circuit $C = C_1 \cdot C_2$, computing $\langle\Psi_0|C^\dagger O C|\Psi_0\rangle$ splits into two halves:

$$\langle\Psi_0|C_2^\dagger \big(C_1^\dagger O C_1\big) C_2|\Psi_0\rangle = \langle\Psi|\tilde{O}|\Psi\rangle, \qquad \tilde{O} = C_1^\dagger O C_1, \quad |\Psi\rangle = C_2|\Psi_0\rangle$$

In the above equation, we can compute $\tilde{\mathcal{O}}$ using Heisenberg propagation, provided natively by propaq, and $|\Psi\rangle$ via Matrix Product States (MPS) using quimb. Finally, we calculate the overlap between the terms of the propagated observable and the MPS. This method is well suited for circuits with bipartite structure, where a certain picture is well suited for each part of the circuit.

Setup¶

A Trotterized transverse-field Ising circuit on $n=16$ qubits: each "layer" applies RZZ across every adjacent pair, followed by RX on every qubit. The observable is a single $Z$ on the middle qubit.

In [1]:

import time

import numpy as np
import qiskit.qasm2 as qasm2
import quimb.tensor as qtn
from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp, Statevector

from propaq.circuits import PauliCircuit
from propaq.datatypes import PauliTermSum
from propaq.hybrid import hybrid_expectation_value
from propaq.truncation import TruncationPolicy
from propaq.propagators.pauli import PauliPropagator

N_QUBITS = 16


def tfim_layer(qc, n, j_coupling, h_field, rng):
    for q in range(n - 1):
        qc.rzz(rng.uniform(0.5, 1.5) * j_coupling, q, q + 1)
    for q in range(n):
        qc.rx(rng.uniform(0.5, 1.5) * h_field, q)


def tfim_circuit(n, layers, seed, j_coupling=0.8, h_field=0.8):
    rng = np.random.default_rng(seed)
    qc = QuantumCircuit(n)
    for _ in range(layers):
        tfim_layer(qc, n, j_coupling, h_field, rng)
    return qc


observable = SparsePauliOp("I" * (N_QUBITS // 2) + "Z" + "I" * (N_QUBITS // 2 - 1))
pauli_observable = PauliTermSum.from_sparse_pauli_op(observable)
truncator = TruncationPolicy(weight_cutoff=10000, coeff_cutoff=1e-12)
propagator = PauliPropagator(None, truncator)

import time import numpy as np import qiskit.qasm2 as qasm2 import quimb.tensor as qtn from qiskit import QuantumCircuit from qiskit.quantum_info import SparsePauliOp, Statevector from propaq.circuits import PauliCircuit from propaq.datatypes import PauliTermSum from propaq.hybrid import hybrid_expectation_value from propaq.truncation import TruncationPolicy from propaq.propagators.pauli import PauliPropagator N_QUBITS = 16 def tfim_layer(qc, n, j_coupling, h_field, rng): for q in range(n - 1): qc.rzz(rng.uniform(0.5, 1.5) * j_coupling, q, q + 1) for q in range(n): qc.rx(rng.uniform(0.5, 1.5) * h_field, q) def tfim_circuit(n, layers, seed, j_coupling=0.8, h_field=0.8): rng = np.random.default_rng(seed) qc = QuantumCircuit(n) for _ in range(layers): tfim_layer(qc, n, j_coupling, h_field, rng) return qc observable = SparsePauliOp("I" * (N_QUBITS // 2) + "Z" + "I" * (N_QUBITS // 2 - 1)) pauli_observable = PauliTermSum.from_sparse_pauli_op(observable) truncator = TruncationPolicy(weight_cutoff=10000, coeff_cutoff=1e-12) propagator = PauliPropagator(None, truncator)

Pure Heisenberg propagation¶

In [2]:

print(f"{'depth':>6} {'n_terms':>9} {'time_s':>8}")
for depth in [1, 2, 3, 4, 5]:
    qc = tfim_circuit(N_QUBITS, depth, seed=1)
    pc = PauliCircuit.from_qiskit(qc)
    t0 = time.perf_counter()
    theta = propagator.propagate(pauli_observable.copy(), pc)
    t1 = time.perf_counter()
    print(f"{depth:6d} {len(theta):9d} {t1 - t0:8.3f}")

print(f"{'depth':>6} {'n_terms':>9} {'time_s':>8}") for depth in [1, 2, 3, 4, 5]: qc = tfim_circuit(N_QUBITS, depth, seed=1) pc = PauliCircuit.from_qiskit(qc) t0 = time.perf_counter() theta = propagator.propagate(pauli_observable.copy(), pc) t1 = time.perf_counter() print(f"{depth:6d} {len(theta):9d} {t1 - t0:8.3f}")

 depth   n_terms   time_s
     1         5    0.067

     2        42    0.138
     3       429    0.176

     4      4862    0.224

     5     58782    0.380

As expected, term growth grows rapidly with depth. Deeper circuits would necessitate aggressive truncation policies, worsening accuracy.

The hybrid split¶

Fix $C_1$ to a small, shallow tail next to the observable, and let $C_2$ grow arbitrarily. $\tilde{O} = C_1^\dagger O C_1$'s term count is now bounded by $C_1$ alone, not by the total circuit depth.

In [3]:

C1_LAYERS = 2
c1 = tfim_circuit(N_QUBITS, C1_LAYERS, seed=100)
pc1 = PauliCircuit.from_qiskit(c1)
theta = propagator.propagate(pauli_observable.copy(), pc1)
print(f"C1 fixed at {C1_LAYERS} layers -> theta has {len(theta)} terms, no matter how deep C2 gets")

C1_LAYERS = 2 c1 = tfim_circuit(N_QUBITS, C1_LAYERS, seed=100) pc1 = PauliCircuit.from_qiskit(c1) theta = propagator.propagate(pauli_observable.copy(), pc1) print(f"C1 fixed at {C1_LAYERS} layers -> theta has {len(theta)} terms, no matter how deep C2 gets")

C1 fixed at 2 layers -> theta has 42 terms, no matter how deep C2 gets

We can do an exact MPS $|\Psi\rangle = C_2 |\Psi_0\rangle$, and then compute the overlap.

In [4]:

print(f"{'total_depth':>11} {'c2_layers':>9} {'max_chi':>7} {'eval_s':>8} {'value':>10}")
for total_depth in [6, 8, 10, 12, 15]:
    c2_layers = total_depth - C1_LAYERS
    c2 = tfim_circuit(N_QUBITS, c2_layers, seed=1)

    circ = qtn.CircuitMPS.from_openqasm2_str(qasm2.dumps(c2))
    max_chi = max(circ.psi.bond_size(i, i + 1) for i in range(N_QUBITS - 1))

    t0 = time.perf_counter()
    value = hybrid_expectation_value(theta, c2, initial_state=0)
    t1 = time.perf_counter()
    print(f"{total_depth:11d} {c2_layers:9d} {max_chi:7d} {t1 - t0:8.4f} {value:10.6f}")

print(f"{'total_depth':>11} {'c2_layers':>9} {'max_chi':>7} {'eval_s':>8} {'value':>10}") for total_depth in [6, 8, 10, 12, 15]: c2_layers = total_depth - C1_LAYERS c2 = tfim_circuit(N_QUBITS, c2_layers, seed=1) circ = qtn.CircuitMPS.from_openqasm2_str(qasm2.dumps(c2)) max_chi = max(circ.psi.bond_size(i, i + 1) for i in range(N_QUBITS - 1)) t0 = time.perf_counter() value = hybrid_expectation_value(theta, c2, initial_state=0) t1 = time.perf_counter() print(f"{total_depth:11d} {c2_layers:9d} {max_chi:7d} {t1 - t0:8.4f} {value:10.6f}")

total_depth c2_layers max_chi   eval_s      value

          6         4       8   0.0268   0.048339
          8         6      16   0.0442   0.025818
         10         8      34   0.0621   0.014218

         12        10      60   1.1013   0.012264

         15        13     157   5.6233  -0.004295

Checking correctness¶

Let's check that the hybrid method is correct.

In [5]:

for c2_layers in [4, 6]:
    c2 = tfim_circuit(N_QUBITS, c2_layers, seed=1)
    hybrid_value = hybrid_expectation_value(theta, c2, initial_state=0)

    full = c2.compose(c1)
    dense_value = Statevector(full).expectation_value(observable).real

    print(
        f"c2_layers={c2_layers}: hybrid={hybrid_value:.8f}  dense={dense_value:.8f}  "
        f"match={np.isclose(hybrid_value, dense_value, atol=1e-6)}"
    )

for c2_layers in [4, 6]: c2 = tfim_circuit(N_QUBITS, c2_layers, seed=1) hybrid_value = hybrid_expectation_value(theta, c2, initial_state=0) full = c2.compose(c1) dense_value = Statevector(full).expectation_value(observable).real print( f"c2_layers={c2_layers}: hybrid={hybrid_value:.8f} dense={dense_value:.8f} " f"match={np.isclose(hybrid_value, dense_value, atol=1e-6)}" )

c2_layers=4: hybrid=0.04833893  dense=0.04833893  match=True

c2_layers=6: hybrid=0.02581783  dense=0.02581804  match=True

When is hybrid simulation useful?¶

The hybrid split is only interesting when neither pure picture is affordable on its own.

In [6]:

print(f"{'n_qubits':>9} {'dense state size':>18}")
for n in [16, 20, 24, 28, 32]:
    n_bytes = (2**n) * 16  # complex128
    if n_bytes < 2**30:
        size = f"{n_bytes / 2**20:.1f} MB"
    else:
        size = f"{n_bytes / 2**30:.1f} GB"
    print(f"{n:9d} {size:>18}")

print(f"{'n_qubits':>9} {'dense state size':>18}") for n in [16, 20, 24, 28, 32]: n_bytes = (2**n) * 16 # complex128 if n_bytes < 2**30: size = f"{n_bytes / 2**20:.1f} MB" else: size = f"{n_bytes / 2**30:.1f} GB" print(f"{n:9d} {size:>18}")

 n_qubits   dense state size
       16             1.0 MB
       20            16.0 MB
       24           256.0 MB
       28             4.0 GB
       32            64.0 GB