Compensated summation for cpu_statevec_anyCtrlAnyTargDenseMatr_sub (#598) - [unitaryHACK]

[mk0dz]

Summary

Adds an optional compensated-summation path to improve numerical accuracy in CPU dense-matrix target evolution.

The implementation introduces a thinner reduction in cpu_statevec_anyCtrlAnyTargDenseMat_sub(), controlled by the compile-time flag:

-DQUEST_COMPENSATED_DENSEMAT_SUM=ON

The feature is disabled by default to preserve existing performance characteristics.

Implementation Notes

• Accumulates the Kahan compensation term into a local cpu_qcomp variable and writes it back once after the reduction.
• Applies only to the CPU dense-matrix path.
• This source file is not built with fast-math, so the compensation term is preserved by the compiler.
• Replaces the previous TODO comment with documentation describing the behaviour and trade-offs.

Files Changed

• quest/src/cpu/cpu_subroutines.cpp
• quest/src/cpu/CMakeLists.txt

Accuracy & Performance

Benchmarked against a __float128 reference using the naive vs. compensated Kahan implementation across all three floating-point precisions.

Precision	Accuracy @ 12 Targets	Runtime Cost
fp32	~58× lower error	~3.3×
fp64	~29× lower error	~2.2× (~1.76× for 10 targets on a 24-qubit state)
fp80	~25× lower error	~2.1×

Observations

• Accuracy improvements grow with target count.
• Benefits become negligible below roughly 7 targets.
• Runtime overhead remains greater than 1× because the kernel is compute-bound on matrix multiplication rather than memory-bandwidth limited.
• The Kahan recurrence prevents efficient vectorization.

Validation

All existing dense CompMatr tests pass for both compensated and uncompensated builds, including:

• applyCompMatr
• Controlled variants

Test results:

• 93,791 assertions
• 4 test cases

Rationale for Default-Off

The improvement is primarily beneficial for large or ill-conditioned matrices, while introducing a consistent ~2–3× runtime overhead.

Keeping the feature disabled by default preserves current performance expectations while allowing users to opt in when higher numerical accuracy is required.

Future Work

A possible alternative is pairwise summation, which offers:

• Error scaling of approximately O(log N)
• Better vectorization opportunities
• Potentially lower runtime overhead

This can be evaluated separately in a future PR.

I also acknowledge claude for guiding me thorough it.

/claim #598

Closes #598.

[TysonRayJones]

Hi Mukul,

Beautiful diff, and very useful and thorough benchmarking!

Super interesting that quadrupling the flops per outer iteration (by the Kahan summation) still measurably impacts runtime at 24 qubits (at which the state is ~120 - 500 MiB). You have demonstrated compensation is only worthwhile for many target qubits (otherwise of course, the sum is too small to need compensation), and the more target qubits, the exponentially greater the total number of flops per iteration. So, there might not exist a regime where compensation of the inner loop has an accuracy benefit, where the outer loop is simultaneously memory bandwidth bound to hide the slowdown. Quite a shock to me!

We will explore further to decide if exposing this hyperparameter (whether or not to compensate) is worthwhile. In any case, you have satisfied the requirements of the unitaryHACK challenge! 🎉 Please comment on issue #598 so I can assign it to you, and award the prize.