Performance improvement of _DiagonalEstimator

[t-imamichi]

Summary

I notice that SamplingVQE is slower than the former VQE when I update a tutorial of Qiskit optimization.
The bottleneck is _DiagonalEstimator._evaluate_sparsepuali. See qiskit-community/qiskit-optimization#448 (comment) for details.
I optimized that part with numpy.

Details and comments

microbenchmark

from timeit import timeit
from qiskit.algorithms.minimum_eigensolvers.diagonal_estimator import _evaluate_sparsepauli
from qiskit.quantum_info import SparsePauliOp
from qiskit.quantum_info.random import random_pauli_list

seed = 123
n = 100
m = 100

op = SparsePauliOp(random_pauli_list(num_qubits=n, size=m, seed=seed))

print(_evaluate_sparsepauli(2<<n-1, op))
print(f'{timeit(lambda: _evaluate_sparsepauli(2<<n-1, op), number=100)} sec')

main

(-2+4j)
0.5635605140000001 sec

this PR

(-2+4j)
0.0018880719999999629 sec

benchmark with QAOA

from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit.utils import algorithm_globals
from qiskit.algorithms.minimum_eigensolvers import QAOA
from qiskit.algorithms.optimizers import COBYLA
from qiskit.primitives import Sampler
from qiskit_optimization.applications import Knapsack

algorithm_globals.random_seed = 123

prob = Knapsack(values=[3, 4, 5, 6, 7], weights=[2, 3, 4, 5, 6], max_weight=10)
qp = prob.to_quadratic_program()

# QAOA
meo = MinimumEigenOptimizer(min_eigen_solver=QAOA(reps=1, sampler=Sampler(), optimizer=COBYLA(maxiter=100)))
result = meo.solve(qp)
print(result.prettyprint())
print("\nsolution:", prob.interpret(result))
print("\ntime:", result.min_eigen_solver_result.optimizer_time)

main

objective function value: 13.0
variable values: x_0=1.0, x_1=1.0, x_2=0.0, x_3=1.0, x_4=0.0
status: SUCCESS

solution: [0, 1, 3]

time: 13.972561120986938

this PR

objective function value: 13.0
variable values: x_0=1.0, x_1=1.0, x_2=0.0, x_3=1.0, x_4=0.0
status: SUCCESS

solution: [0, 1, 3]

time: 0.8536498546600342

[qiskit-bot]

Thank you for opening a new pull request.

Before your PR can be merged it will first need to pass continuous integration tests and be reviewed. Sometimes the review process can be slow, so please be patient.

While you're waiting, please feel free to review other open PRs. While only a subset of people are authorized to approve pull requests for merging, everyone is encouraged to review open pull requests. Doing reviews helps reduce the burden on the core team and helps make the project's code better for everyone.

One or more of the the following people are requested to review this:

• @Qiskit/terra-core
• @manoelmarques
• @woodsp-ibm

[coveralls]

Pull Request Test Coverage Report for Build 3618715227

• 5 of 5 (100.0%) changed or added relevant lines in 1 file are covered.
• 2 unchanged lines in 1 file lost coverage.
• Overall coverage increased (+0.04%) to 84.598%

---

Files with Coverage Reduction	New Missed Lines	%
src/sabre_swap/layer.rs	2	98.95%

Totals
Change from base Build 3618561726:	0.04%
Covered Lines:	63148
Relevant Lines:	74645

---

💛 - Coveralls

[t-imamichi]

It would be great if this PR is backport stable because Qiskit optimization relies on SamplingVQE.

[Cryoris]

LGTM, but I'll let @mtreinish or @jakelishman comment on the backporting part 🙂 From my side backporting would be fine considering this a "performance bug".... 😁

[jakelishman]

I think we can include this in a backport safely enough, as long as there's direct test coverage that the change produces correct results.

Every time I get tagged in one of these I always end up getting performance nerd-sniped. Here's my offering:

_PARITY = np.array([-1 if bin(i).count("1") % 2 else 1 for i in range(256)], dtype=np.complex128)

def _evaluate_sparsepauli(state: int, observable: SparsePauliOp) -> complex:
    packed_uint8 = np.packbits(observable.paulis.z, axis=1, bitorder="little")
    state_bytes = np.frombuffer(state.to_bytes(packed_uint8.shape[1], "little"), dtype=np.uint8)
    reduced = np.bitwise_xor.reduce(packed_uint8 & state_bytes, axis=1)
    return np.sum(observable.coeffs * _PARITY[reduced])

For small numbers of qubits and operators in the list, it's about ~10-20% faster, but say for an observable with 100q and 100 Paulis, it's nearly 10x faster due to avoiding the Python object conversions. For the optimisation code block given in the top comment, on my machine the evaluation went from 11.5s on main to 1.1s with the version in this PR, and about 0.8s with my version of this function (I didn't time very accurately).

[t-imamichi]

Thank you, Jake. Your code is much faster than mine. Cool!
I added a short reno too.

[Cryoris]

LGTM thanks for the improvement (and nerd-sniping 😉)!