EstimatorV2: O(N²) post-loop and fork-unsafe DEFAULT_EXECUTOR cause severe slowdowns

[Siddharthgolecha]

Informations

• Qiskit Aer version: 0.17.2
• Qiskit version: 2.4.1
• Python version: 3.13
• Operating system:
  • User's reported environment: Red Hat Enterprise Linux 9, x86_64, 16-core CPU
  • Reproduction environment: macOS 26.4, arm64, 16-core CPU (same Aer / Qiskit versions). Both platforms exhibit the same bugs; absolute timings differ but the mechanism is identical.

What is the current behavior?

EstimatorV2 exhibits two independent, simultaneous bugs that together turn a workload that should take ~1 minute (per-process parallel) into a ~17 minute single-threaded run on a 15-qubit, 54-parameter HHFM circuit with 100,000 parameter rows and 15 observables, MPS backend, 16-core machine:

configuration	wall time	speedup vs single-PUB
1 PUB, internal parallelism (max_parallel_threads=N_CPUS)	~17 min	1×
16 PUBs, internal parallelism	~3 min	~5.5×
16 PUBs, external ProcessPoolExecutor, max_parallel_threads=1	~60 s	~17×

Internal parallelism with max_parallel_threads=N_CPUS should approach the per-process result. It does not.

Bug #1 — O(N²) post-processing loop in EstimatorV2._run_pub

qiskit_aer/primitives/estimator_v2.py:138,143:

flat_indices = list(param_indices.ravel())          # 100k-element list of tuples
...
for index in np.ndindex(*bc_param_ind.shape):       # 100k * 15 = 1.5M iterations
    param_index = bc_param_ind[index]
    flat_index = flat_indices.index(param_index)    # O(N) linear scan per iter
    ...

For 100k rows × 15 observables this performs ~7.5 × 10¹⁰ comparisons — about 750 s of pure-Python work at ~10 ns/op, matching the ~530 s post-loop time observed in py-spy flame graphs of the 1033 s internal_sv run.

py-spy flame graph (worker thread, 6,742 samples, 67 s wall):

samples	% of worker	location
3,420	50.7 %	estimator_v2.py:143 — flat_indices.index(...)
1,998	29.6 %	cpp_execute_circuits (the actual C++ simulation)

The pure-Python post-loop (51 %) costs more wall time than the C++ simulation (30 %). Splitting the same 100k-row work into 16 PUBs drops the post-loop's share from 51 % → 7 %, confirming it scales worse than the simulation with n_lines.

Bug #2 — DEFAULT_EXECUTOR is not fork-safe

qiskit_aer/jobs/utils.py:19:

DEFAULT_EXECUTOR = ThreadPoolExecutor(max_workers=1)

This module-scope singleton is constructed at import time. In any long-lived Python process (Jupyter kernel, IPython, services), ordinary kernel activity warms the executor — its worker thread is alive and its internal queue has been touched — before the user's code runs.

os.fork() then copies the executor object's memory (including its lock state and queue) but only the calling thread — the worker thread does not survive. Every forked child inherits a DEFAULT_EXECUTOR that looks normal but has no live worker. The first call into EstimatorV2.run() enqueues to a queue nothing will ever consume, and the relevant call path raises during setup.

Direct evidence from a py-spy record --threads --idle --subprocesses capture of a live Jupyter kernel running the user's notebook (1,284,652 raw samples):

process	role	samples	share
parent kernel	doing all the work alone, single-threaded	1,257,037	97.9 %
16 forked workers	combined	27,615	2.1 %

Every one of the 16 children's leaf stack ends in:

... → ProcessPoolExecutor.submit → _Popen → _launch (popen_fork.py:74)
   ↑ os.fork() returns; child resumes below ↑
   → _process_worker (concurrent/futures/process.py:256)
   → _ExceptionWithTraceback.__init__ (concurrent/futures/process.py:137)
   → format_exception (traceback.py:154) → linecache.updatecache

That call path is reachable only from the except branch of _process_worker. Every child's run(estimator, [pub]) raised; the worker is now formatting the traceback (slow because linecache reads source files line by line). Across all 16 children:

• fewer than 100 samples reach qiskit_aer code total,
• zero samples reach cpp_execute_circuits — not a single child ever ran any quantum simulation.

So the user's "external parallelisation" cell was a no-op for parallelism. The 16 workers were ghosts; the parent did 100 % of the work alone, single-threaded, hitting the O(N²) post-loop from bug #1. From the outside it looked like "16 workers running"; in reality the workers contributed zero useful work.

The standalone-script reproduction works (~60 s) because nothing in a fresh Python process warms the executor's worker thread before fork. In Jupyter the kernel warms it implicitly before any cell runs, so the same code path silently breaks.

Steps to reproduce the problem

A complete reproduction (benchmark.py and run_bench.py) is available; minimal repro for each bug below.

Bug #1 — O(N²) post-loop

import numpy as np
from qiskit import QuantumCircuit
from qiskit.circuit import Parameter
from qiskit.quantum_info import SparsePauliOp
from qiskit.primitives.containers.estimator_pub import EstimatorPub
from qiskit_aer.primitives import EstimatorV2

n_qubits = 4
n_rows = 50_000  # already painful; user hit 100k

theta = [Parameter(f"t{i}") for i in range(n_qubits)]
qc = QuantumCircuit(n_qubits)
for i, p in enumerate(theta):
    qc.ry(p, i)

obs = [SparsePauliOp.from_list([("Z" * n_qubits, 1.0)])]
params = np.random.default_rng(0).uniform(-1, 1, size=(n_rows, n_qubits))

est = EstimatorV2()
pub = EstimatorPub.coerce((qc, obs, params))
import time; t0 = time.perf_counter()
est.run([pub]).result()
print(f"wall: {time.perf_counter() - t0:.1f} s")
# Wall time grows quadratically with n_rows. Profile with py-spy and
# you will see >50% of samples in estimator_v2.py:143.

Bug #2 — DEFAULT_EXECUTOR after fork

import multiprocessing as mp
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator

def child_work():
    qc = QuantumCircuit(2); qc.h(0); qc.cx(0, 1); qc.measure_all()
    AerSimulator().run(qc).result()
    print("child done")

if __name__ == "__main__":
    # Warm DEFAULT_EXECUTOR in the parent (any Aer .run() call does this;
    # in Jupyter it happens automatically before any user cell runs).
    qc = QuantumCircuit(2); qc.h(0); qc.cx(0, 1); qc.measure_all()
    AerSimulator().run(qc).result()

    ctx = mp.get_context("fork")
    p = ctx.Process(target=child_work)
    p.start(); p.join(timeout=30)
    print("exitcode:", p.exitcode)
    # Without the fix: child raises in EstimatorV2.run() / AerJob.submit()
    # because the inherited ThreadPoolExecutor has no live worker thread.

To see the full Jupyter manifestation, run a notebook cell that warms an AerSimulator, then a second cell that submits 16 EstimatorV2.run() calls via ProcessPoolExecutor(mp_context="fork"). The parent will silently do all the work; the children will sit in format_exception with zero useful CPU work.

What is the expected behavior?

1. For bug Error with OpenMP #1: EstimatorV2._run_pub should perform parameter-index lookup in O(1), not O(N). The post-loop should scale linearly with n_rows × n_observables, not quadratically. On the 100k × 15 workload the post-loop should consume <5 % of wall time (vs the current 51 %), and internal_sv wall should drop from ~17 minutes to ~8 minutes (with the C++ simulation now correctly identified as the bottleneck).

2. For bug Update README and build instructions on Mac #2: mp_context="fork" workers should be able to use Aer's AerJob / EstimatorV2 without silently failing. A user submitting 16 jobs via ProcessPoolExecutor(mp_context="fork") should actually get 16× parallelism — not a serial parent doing all the work while 16 children format tracebacks.

Suggested solutions

Fix #1 — Replace the O(N²) list scan with a dict lookup + 1D fast path

# qiskit_aer/primitives/estimator_v2.py

flat_index_map = {tup: i for i, tup in enumerate(param_indices.ravel())}
evs = np.zeros_like(bc_param_ind, dtype=float)
stds = np.full(bc_param_ind.shape, precision)
for index in np.ndindex(*bc_param_ind.shape):
    param_index = bc_param_ind[index]
    if bc_param_ind.shape == param_shape:  # 1d passthrough
        flat_index = param_index[0] if isinstance(param_index, tuple) else int(param_index)
    else:
        flat_index = flat_index_map[param_index]
    for pauli, coeff in bc_obs[index].items():
        expval = result.data(flat_index)[pauli]
        evs[index] += expval * coeff

Pure bugfix — semantics identical, numerical output bit-for-bit unchanged.

Fix #2 — Reset DEFAULT_EXECUTOR in every forked child

# qiskit_aer/jobs/utils.py
import os
from concurrent.futures import ThreadPoolExecutor

DEFAULT_EXECUTOR = ThreadPoolExecutor(max_workers=1)


def _reset_default_executor_in_child():
    global DEFAULT_EXECUTOR
    DEFAULT_EXECUTOR = ThreadPoolExecutor(max_workers=1)


if hasattr(os, "register_at_fork"):
    os.register_at_fork(after_in_child=_reset_default_executor_in_child)

~10 lines, no public API change, no caller changes; from qiskit_aer.jobs.utils import DEFAULT_EXECUTOR continues to work because the module-level binding is updated in place. Guarded with hasattr so the import is safe on platforms without os.register_at_fork (Windows).

I have implemented and tested both fixes locally and am preparing a PR.

Workarounds for users until the fixes land

1. Always chunk parameter rows. A 16-PUB split (np.array_split(data, N_CPUS)) cuts the post-loop O(N²) cost by 16² = 256× and explains the ~5× speedup observed for the multi-PUB internal case.
2. Prefer mp_context="spawn" over "fork" in any long-lived Python process (Jupyter, IPython). Spawn workers do not inherit the parent's DEFAULT_EXECUTOR state at all. Cost: ~1–2 s of import overhead per worker.
3. Run benchmarks in standalone scripts, not Jupyter cells, when measuring runtime. Jupyter kernels accumulate threadpool state across cells; this can silently turn parallel cells into serial ones (bug Update README and build instructions on Mac #2) and change the timing of subsequent cells in non-obvious ways.

References

• qiskit_aer/primitives/estimator_v2.py:138,143 — bug Error with OpenMP #1 location
• qiskit_aer/jobs/utils.py:19 — bug Update README and build instructions on Mac #2 location
• qiskit_aer/jobs/aerjob.py:61 — caller of DEFAULT_EXECUTOR, unchanged
• Full diagnostic write-up (profiling methodology, flame graphs, Jupyter
  py-spy capture, platform comparison data) available on request.

[cottephi]

Hi,
Thank you for creating the issue and PR.
Is there any progress on that ?