Allow calloc-like optimisations in SparsePauliOp allocations

[jakelishman]

As of Rust 1.98, vec![0.0; 2 * len] lowers to OS/kernel-level "allocate zero-filled" functions, but using a non-primitive type in the Vec (like Complex64) doesn't permit this optimisation. We want to use calloc instead of manually writing zeros into the vector because it often gets kernel-level support to have zero impact on pages that are never written to, whereas manually issuing write calls across each row guarantees that we have to page in everything.

This has a larger performance impact for large dense matrices with high sparsity factors.

Using a benchmark script:

import timeit
import numpy as np
from qiskit.quantum_info import SparsePauliOp

def make_sparse_list(num_qubits, num_terms, locality, rng):
    terms = []
    for _ in range(num_terms):
        indices = rng.choice(num_qubits, size=locality, replace=False)
        label = "".join(rng.choice(tuple("XYZ")) for _ in indices)
        coeff = complex(rng.normal(), rng.normal())
        terms.append((label, indices, coeff))
    return terms

rng = np.random.default_rng(42)
print(f"{'n':>3} {'terms':>6} {'locality':>8} {'time (s)':>10}")
configs = [
    (6, 10, 2),
    (10, 10, 2),
    (10, 50, 2),
    (10, 100, 4),
    (14, 10, 2),
    (14, 50, 2),
    (14, 100, 4),
]
for num_qubits, num_terms, locality in configs:
    terms = make_sparse_list(num_qubits, num_terms, locality, rng)
    op = SparsePauliOp.from_sparse_list(terms, num_qubits=num_qubits)
    reps, tot = timeit.Timer(op.to_matrix).autorange()
    print(f"{num_qubits:>3} {num_terms:>6} {locality:>8} {tot / reps:>10.4f}")

On main (as of aae0d80), on a Macbook Pro M4 Max using 16 threads, the timings were:

  n  terms locality   time (s)
  6     10        2     0.0000
 10     10        2     0.0003
 10     50        2     0.0003
 10    100        4     0.0003
 14     10        2     0.2860
 14     50        2     0.2767
 14    100        4     0.2757

With both this commit and gh-16828 in place, the numbers instead become:

  n  terms locality   time (s)
  6     10        2     0.0000
 10     10        2     0.0001
 10     50        2     0.0001
 10    100        4     0.0002
 14     10        2     0.0229
 14     50        2     0.0313
 14    100        4     0.0459

AI/LLM disclosure

• No part of this submission is LLM generated.
• Some written text was generated by:
• Some submitted code was generated by: Claude Sonnet 5

The benchmark script in the commit message is a significantly modified version of a script Julien generated with Claude. Nothing else is LLM.

[qiskit-bot]

One or more of the following people are relevant to this code:

• @Qiskit/terra-core

[jakelishman]

Ha - totally forgot that bytemuck would provide an explicit interface for this via the Zeroable trait. Still happy I did the discovery about the vec![0.0; len] optimisation, though.

[jeevan0920]

Just curious, how did you discover this?

[jakelishman]

It comes very indirectly from #16809 - David and Julien had been looking at the implementation, doing some benchmarking, and Julien had discovered a weirdness where SparsePauliOp.to_matrix(sparse=True).toarray() sometimes significantly beat the direct dense output. That was super suspicious, since the conversion from sparse-matrix to dense-matrix form here should be strictly a subset of the runtime.

After that, it just seemed most likely that the weirdo ad-hoc initialisation was at fault here; the raw calloc is highly optimised and gets loads of tricks that can make zeroing appear to be much cheaper than it "should" be in certain cases.