feat: support Windows with a Rust Davidson eigensolver

[Mostafa-Atallah2020]

Closes #84.

qiskit-addon-slc was Linux/macOS-only because of pyscf, which has no Windows wheel and cannot be built from source there. It was used in a single place: the pyscf.lib.davidson1 call in get_extremal_eigenvalue. This replaces that with a compiled Rust Davidson eigensolver (built on nalgebra, no BLAS/LAPACK), drops pyscf entirely, and adds a windows-latest runner to the test matrix so the suite is exercised on Windows in CI.

cc @mrossinek @aeddins-ibm - this is the Rust route you suggested on #88. Compared with #85 (which made pyscf optional) it removes pyscf outright and proves the full workflow runs on Windows in CI, compared with the numpy/scipy version in #88 it keeps pyscf-level speed without any dependency.

[coveralls]

Coverage Report for CI Build 33094731853

Coverage decreased (-0.2%) to 76.029%

Details

• Coverage decreased (-0.2%) from the base build.
• Patch coverage: 1 uncovered change across 1 file (10 of 11 lines covered, 90.91%).
• No coverage regressions found.

Uncovered Changes

File	Changed	Covered	%
qiskit_addon_slc/utils/davidson.py	11	10	90.91%

Coverage Regressions

No coverage regressions found.

---

Coverage Stats

Relevant Lines:	826
Covered Lines:	628
Line Coverage:	76.03%
Coverage Strength:	0.76 hits per line

---

💛 - Coveralls

[mrossinek]

Thanks a lot, this already looks like it is in great shape! I made some adjustments mostly to the repo organization to follow other Qiskit addons more closely. I also ran a /code-review with Claude. I will have the findings of that reported in a minute 👍

[mrossinek]

Review of the Rust Davidson eigensolver. Three of the inline findings compound into what I think is a release blocker, so I'll summarize the interaction here and leave the details on the lines.

The compounding failure

lib.rs:98 (stall criterion), lib.rs:107 (preconditioner clamp), and lib.rs:124 (cnorm < lindep → converged = true) interact badly on exactly the inputs this package solves. For the Pauli commutators forward.py passes in, the matrix diagonal is entirely zero (1287/1313 sampled operators), so cycle 0 gives theta ≈ 1e-16, the clamp flips the sign of the preconditioner shift, theta stays pinned near zero, and the stall criterion then fires on cycle 1 and returns (true, ~0).

90.3% of realistic commutator inputs (1186/1313) return a grossly wrong eigenvalue flagged as converged. Because forward.py:189 trusts the converged flag and skips the triangle-inequality fallback, the resulting error bounds are under-estimated — unsafe in the wrong direction for this package's purpose. The old pyscf implementation returned wrong values on some of these too, but reported converged=False, so the fallback caught them. The flag flip is the regression.

Test fixtures bake in the regression

tests/expected_fwd_bounds.pickle (binary, so no inline comment possible): 1080 of 4632 expected values changed, moving from 2.000000 — the triangle-inequality fallback, used when pyscf correctly reported non-convergence — to smaller values such as 1.111140, produced by the new solver claiming convergence with a too-small eigenvalue. The suite therefore passes while the bounds are wrong. I'd suggest not regenerating these fixtures until the solver is fixed, then diffing against the pre-regeneration values to confirm they come back to the fallback where they should. Same concern applies to expected_fwd_tightened_bounds.pickle and expected_merged_bounds.pickle if they were regenerated in the same pass.

Suggested ordering

1. lib.rs:107 — preserve the sign of d, use a relative floor.
2. lib.rs:98 — drop the stall criterion, keep only the residual test.
3. lib.rs:124 — converged = false (or restart with a fresh random vector).
4. Regenerate fixtures and confirm the reverted values.

The remaining inline comments are independent robustness items in the new Rust code.

[mrossinek]

Critical: |eigval - prev| < tol declares false convergence.

This stall criterion fires while the residual is still large — a Davidson iteration whose Ritz value stops moving is not necessarily near an eigenvalue, it may just be stuck (see the preconditioner issue on line 107, which pins theta near zero and makes this fire on cycle 1).

It is the dominant exit path: 131 of 200 random test operators leave the loop here rather than through the residual test. Removing it and keeping only residual.norm() < tol drops the worst converged-result error from 3e-3 to 2.7e-8 on a 200-operator sweep at tol=1e-6.

Suggested change

	if (eigval - prev).abs() < tol || residual.norm() < tol {
	if residual.norm() < tol {

With the stall criterion gone, prev becomes unused and can be dropped along with line 102.

[mrossinek]

Critical: the clamp on line 108-110 corrupts zero-diagonal operators.

For the Pauli commutators forward.py actually passes in, the diagonal is entirely zero — 1287 of 1313 sampled operators. That makes cycle 0 give theta ≈ 1e-16, so:

d = diag[i] - theta ≈ -1.4e-16
d.norm() < tol          → true
d = +tol                → sign flipped, magnitude inflated ~1e8

The correction is driven in the wrong direction, theta stays pinned at ~0, and the stall criterion on line 98 then fires on cycle 1 and reports success.

Reproducer:

op = SparsePauliOp(
    ['YYZYXIZ', 'YXYZZZZ', 'YIZXIXZ', 'XYYXXIY'],
    [0.976j, -1.107j, -0.126j, -1.179j],
)
# returns (True, -2.1e-16); true minimum is -2.596

90.3% of realistic commutator inputs (1186/1313) return a grossly wrong eigenvalue flagged as converged.

The clamp needs to preserve the sign of d and use a floor relative to the operator scale rather than the convergence tolerance — tol is a residual threshold and has no business as a magnitude floor on a matrix entry. Something along these lines:

let mut d = diag[i] - C64::new(theta, 0.0);
let floor = 1e-12 * scale.max(1.0);  // scale from e.g. diag.amax() or ||A||
if d.norm() < floor {
    d = C64::new(if d.re < 0.0 { -floor } else { floor }, 0.0);
}

For an all-zero diagonal the preconditioner is degenerate no matter how it's clamped, so it may be cleaner to detect that case up front and skip preconditioning entirely (correction = residual), which is a valid, if slower, Davidson variant.

[mrossinek]

Major: cnorm < lindep should not set converged = true.

A correction vector that vanishes after orthogonalization means the subspace can no longer be expanded — the iteration is stuck, not converged. Reporting success here converts what was a safe fallback into a silently wrong bound.

Any purely diagonal operator (I/Z Paulis only) hits this:

op = SparsePauliOp(
    ['ZIIII', 'IZIII', 'IIZII', 'ZZIII', 'IIIZI'],
    [1.0, 0.7, 0.4, 0.3, 0.9],
)
# returns (True, 0.185); true minimum is -2.7

I checked the old pyscf implementation on this input: it returned the same wrong value, but with converged=False, so forward.py:189 fell back to the triangle inequality and the bound stayed valid. The flag flip is the regression, not the numerics.

Suggested change

	if cnorm < lindep {
	if cnorm < lindep {
	converged = false;

Restarting with a fresh random vector (bounded by a retry count) would recover more of these cases, but simply reporting non-convergence restores the safe behaviour.

[mrossinek]

No CSR length or bounds validation — malformed input panics instead of raising.

indptr[row + 1] is read for row up to dim - 1 with no check that indptr.len() == dim + 1, and indices[k] indexes x unchecked. Since dim arrives as a parameter independent of the arrays (davidson_smallest line 148), a caller can easily desynchronize them.

Verified failure modes, all surfacing as pyo3_runtime.PanicException rather than a Python exception:

• indptr shorter than dim + 1 → index out of bounds: the len is 2 but the index is 2
• indices[k] >= dim → out-of-bounds on x
• diag or seed shorter than dim → panic in the preconditioner / DVector construction

A panic across the FFI boundary is much worse than an exception — it's not catchable as a normal error and the message is opaque. Worth validating in davidson_smallest and returning PyValueError: indptr.len() == dim + 1, indptr non-decreasing with indptr[dim] == indices.len() == data.len(), every indices[k] < dim, and diag.len() == seed.len() == dim.

[mrossinek]

max_space <= 1 silently degenerates rather than erroring.

With max_space <= 1 the condition s.len() >= max_space is true on every cycle, so the subspace collapses back to a single vector before it can ever grow to two. A 2D subspace never forms, Rayleigh-Ritz reduces to the Rayleigh quotient of the current vector, and the result is reported as converged.

max_space is user-facing through get_extremal_eigenvalue, so this is reachable from the public API. Worth rejecting max_space < 2 with a PyValueError in davidson_smallest alongside the CSR validation.

[mrossinek]

dim = 0 returns (True, 0.0).

An empty matrix has no eigenvalues, but the loop body never executes in a way that changes eigval, so the function reports a fabricated converged eigenvalue of exactly 0.0. That is indistinguishable from a legitimate result and will propagate into a bound. Reject dim == 0 with a PyValueError — cheap to add next to the other validation.

[mrossinek]

Single classical Gram-Schmidt pass, though the comment on line 120 says "modified Gram-Schmidt".

correction -= &s_mat * (s_mat.adjoint() * &correction) projects against all basis vectors simultaneously using the original correction — that is classical Gram-Schmidt, not MGS, which subtracts sequentially and re-reads the partially-updated vector each time.

With a subspace up to max_space vectors, one classical pass loses orthogonality progressively, which degrades the Rayleigh-Ritz projection in later cycles. Two options: repeat the projection twice (classical Gram-Schmidt with re-orthogonalization, which is numerically comparable to MGS and keeps the efficient matrix form), or loop the columns for true MGS. Either way the comment should match what the code does.

Suggested change

	correction -= &s_mat * (s_mat.adjoint() * &correction);
	correction -= &s_mat * (s_mat.adjoint() * &correction);
	correction -= &s_mat * (s_mat.adjoint() * &correction);