Perf: reduced-op eigensolver

[aeddins-ibm]

In some cases, the main bottleneck in the SLC computation is getting the spectral norm (i.e. largest magnitude eigenvalue) of the commutator [error, observable]. This matrix is sparse, but still potentially really big. This PR provides an optimization so we can work with a smaller matrix (fewer effective qubits) that still captures all distinct eigenvalues of the original operator. The main work of this PR is in the _reduce_operator function.

Additionally:

• adds new logic for whether to call the usual Davidson solver, vs to use a dense solver, vs to do something intermediate, depending on operator properties.
• lowers the default tolerance on Davidson from 1e-6 to 1e-10, as some tests yielded eigenvalue error as large as ~1e-4 with the previous setting.
• adds tests for the new functions. Some old test values changed because the result converges now where it did not previously. All results in the modified pickle files have been verified against exact dense solver (numpy eigvalsh).

On a 1D 50-qubit Ising circuit with rx angle of pi/16 and rzz angle of -pi/2, this PR enabled compute_forward_bounds to obtain several more layers of bounds within a 5-minute timeout window (up to ~16k error-generators processed from an original ~14k -- significant given subsequent layers are potentially exponentially harder to process).

(This PR was developed with assistance from AI tool Claude. This PR summary write-up is not AI.)

_reduce_operator(SparsePauliOp)

For an operator with linearly-independent Pauli terms (i.e. a full rank tableau), the unique eigenvalues are set by the pairwise commutation relations between terms (and the coefficients). So one goal condition is to produce an operator with the same pairwise commutation relations, but fewer qubits. This removes degeneracies in which the same eigenvalue appears multiple times in the spectrum. For a trivial example, if your original operator Paulis were (XXX, ZZZ), you'd get the same pairwise commutation relations by the new Paulis (X, Z). The general procedure is more complicated.

We will reduce initial Paulis and coefficients {(P_i, c_i)} to fewer-qubit Paulis and coeffs {(R_i, r_i)} (where r_i is same as c_i up to a possible phase).

Steps:

1. Identify a minimal basis of "generator" Paulis {G_j}, building-blocks we can multiply together to produce any P_i. This shrinks the number of Paulis we need to deal with. If we can express these building-blocks as Paulis on fewer qubits, then we can use them to build a fewer-qubit operator that is equivalent to the original (same commutation relations between terms).
   • Do this by row-reducing the symplectic tableau (with mod-2 math since Pauli multiplication is XOR of these bits), such that all terms are linearly independent.
2. Represent each generator as a weight-1 Pauli on a new set of "reduced qubits" ("logical qubits"). An anticommuting pair of Paulis can be mapped to (X, Z) on the same qubit, reducing the qubit number.
   • Loop: Pop a generator A from list of generators.
     • If it commutes with all remaining generators, represent it as a Z operator on a new qubit. We won't map any other generators to this qubit, guaranteeing the commutation relations are preserved.
     • Otherwise, pop an anticommuting generator B from the list. Represent (A, B) as (Z, X) on a new qubit. So that anticommutation property is preserved.
       • Subtle part: update our list of remaining generators {G} so that they all commute with both A and B. If G anticommutes with A (or B), multiply G in-place by B (or A). The updated generators' tableaux still span the same space, so they still work as building-blocks for rebuilding our original operator. But since G now commutes with A and B, we are free to assign G to a new reduced qubit. (This is the "symplectic Gram Schmidt" part).
   • Repeat loop until all generators consumed.
3. We now have a list of weight-1 Paulis that isomorphically corresponds to the many-qubit generator building-blocks. Using this correspondence, we identify which generators in our list are needed to construct each Pauli in our original operator, and we accordingly multiply the corresponding weight-1 generators to produce fewer-qubit Paulis that isomorphically corresponds to our original operator. (This comparison also may yield a phase correction of the coefficient, which is why r_i can differ from c_i).

Choice of solver (secondary optimization)

The simple procedure is to pass the reduced operator to our existing Davidson solver. Since there are now fewer qubits, this can be significantly faster than without this PR.

For a small reduced operator, we do an additional optimization. In step 2 above, we identified commuting ("central") generators, corresponding to Z on a reduced qubit, that generally commute with everything else. Say the loop produced c such generators, and p pairs of anticommuting generators. We know each commuting Z will be replaced by an eigenvalue +1 or -1 everywhere it appears in the operator, but we don't know which. We can optionally plug in all 2^c combinations, then eigensolve 2^c smaller operators instead of a single larger operator. This can be faster when p+c is small.

[coveralls]

Coverage Report for CI Build 33086560019

Coverage increased (+2.5%) to 78.737%

Details

• Coverage increased (+2.5%) from the base build.
• Patch coverage: 5 uncovered changes across 1 file (121 of 126 lines covered, 96.03%).
• No coverage regressions found.

Uncovered Changes

File	Changed	Covered	%
qiskit_addon_slc/utils/reduce_op.py	96	91	94.79%
Total (2 files)	126	121	96.03%

Coverage Regressions

No coverage regressions found.

---

Coverage Stats

Relevant Lines:	950
Covered Lines:	748
Line Coverage:	78.74%
Coverage Strength:	0.79 hits per line

---

💛 - Coveralls

[aeddins-ibm]

Compared to 1 month ago, main subroutines are now much more readable, at least for someone versed in qiskit.quantum_info.

Custom function names are more descriptive; some custom code has been replaced with quantum_info methods.

Overall I think reduce_op.py (previously the more mysterious part) is in good shape. The higher-level davidson.py module still has a bit more slop-comments than ideal, but I think it's clear enough.

[mrossinek]

Thanks a lot! I think I can follow your explanation in the PR and agree with the approach taken here. The code also looks reasonable to me, but it was not super easy to account for all potential edge cases. For this, I told Claude to do a review of this PR and here are its findings:

---

HIGH — qiskit_addon_slc/utils/davidson.py:94: lowest = 0.0 initialization clamps the result

The dense per-sector loop seeds lowest = 0.0 and then only ever takes min(lowest, ...).
If every sector's smallest eigenvalue is positive, the true minimum is never recorded and the
function returns 0.0.

spo = SparsePauliOp(["XI", "ZI", "II"], [0.5, 0.5, 3.0])
get_extremal_eigenvalue(spo)  # -> (True, 0.0);  true minimum is +2.2928932188134525

Fuzzing 300 random Hermitian operators hit this 8 times — any operator with an identity term
large enough to shift the spectrum positive triggers it. The correct seed is +inf (or
float(np.linalg.eigvalsh(...)[0]) on the first sector).

get_extremal_eigenvalue is a documented public API (in qiskit_addon_slc/utils/__init__.py's
__all__, with an autofunction entry) and its docstring promises "the most-negative eigenvalue
of spo" for any Hermitian input. The addon's own caller only ever passes traceless commutators,
so production paths are unaffected — but the exported contract is broken.

MEDIUM — qiskit_addon_slc/utils/davidson.py:74: the p == 0, c == 0 case returns 0.0, discarding the operator's constant value

When spo is a multiple of the identity, _reduce_operator returns a 0-qubit SparsePauliOp.
Qiskit 2.5.2's to_matrix() on a 0-qubit operator yields [[0.+0.j]], dropping the coefficient,
so .diagonal().real.min() is 0.0:

get_extremal_eigenvalue(SparsePauliOp(["II"], [-0.73]))  # -> (True, 0.0); true minimum is -0.73

This is also the only case a 400-seed _reduce_operator spectrum fuzz flagged, confirming the
reduction itself is otherwise sound. forward.py:139-144 short-circuits identity-only commutators
before this call, so the addon's own path is safe, but the public API is not. Guarding on
reduced_op.num_qubits == 0 and returning float(sum(coeffs).real) fixes it.

LOW — tests/utils/test_davidson.py:41: the new test fixtures all have negative minima, so neither bug above is covered

All four DENSE_OPERATORS entries plus both other new tests have a strictly negative minimum
eigenvalue, so the lowest = 0.0 clamp is invisible to the suite. Adding one positive-definite
fixture (e.g. SparsePauliOp(["XI", "ZI", "II"], [0.5, 0.5, 3.0])) and one identity-only fixture
would catch both.

LOW — qiskit_addon_slc/utils/davidson.py:111: the new tol=1e-10 default is dead in the addon's own code path

forward.py:187 calls get_extremal_eigenvalue(commutator, tol=atol_eigenvalue), and
atol_eigenvalue defaults to 1e-8, so default_kwargs.update(kwargs) always overrides the new
tighter default. The docstring's rationale ("the default tol is tight because the eigenvalue
error runs well above tol") therefore does not apply to any in-repo caller; if the tightening is
meant to take effect for bound computations, atol_eigenvalue's default needs the same change.

Verified as correct (no finding)

Fuzz-tested the riskiest parts of the new reduce_op.py and found them sound:

• _get_basis maintains proper reduced-row-echelon invariants (300 seeds)
• _symplectic_gram_schmidt's recursive center resolution preserves the pair/center commutation
  structure, and _get_basis(center) re-reduction does not break it (400 seeds)
• Pauli phases are correctly recovered through the omega_phase computation, including signed and
  imaginary input labels
• The reduced operator always keeps the trailing c qubits Z-only, which the sector-sign logic in
  davidson.py depends on (500 seeds)
• c <= n always holds, so the p == 0 sparse-diagonal memory claim is valid
• bool(np.atleast_1d(converged)[0]) correctly unwraps davidson1's ndarray return
• The in-loop block.coeffs reassignment is picked up by to_matrix() with no stale caching
• The iterative path matches dense diagonalization on 60 traceless operators