Sparse eigsh follow-up #10257

I would not approve [PR #10098](#10098) at current head b4ad9547f. CI is green, but I found two correctness concerns. No security issue found.

1. Blocking: repaired vectors are treated as converged without representing their residual

Comment on [_repair_locked](https://github.com/cupy/cupy/blob/b4ad9547f8f5321de89bfb00907c0ab8bc12a874/cupyx/scipy/sparse/linalg/_eigen.py#L629-L635):

Blocking: replacing V[j] with an arbitrary _restart_ortho direction and updating only alpha[j] does not make it a valid locked 1×1 Ritz block. For the new vector q, A @ q generally has components outside q, but beta[:k] remains zero and those components are never represented in T. Moreover, beta_k[j] still describes the old duplicated Ritz vector; it is subsequently used both to deflate u and as the Ritz arrowhead.

Therefore T != V A Vᴴ, and the existing residual estimate can declare the replacement converged while ||Aq - alpha[j]q|| is nonzero. For example, with A = diag(3,2,1) and duplicated locked rows [e0,e0], this restart construction produces approximately q = -0.968e1 - 0.252e2, whose Rayleigh quotient is 1.9366 but whose true residual is 0.244; the decoupled projected matrix nevertheless treats 1.9366 as exact.

I think a duplicate needs to be removed from the locked set and expanded through A again, or the complete projected operator and true residual must be rebuilt. Recomputing beta_k[j] alone is insufficient. Please also add a regression that forces this branch and checks ||Ax-xw||, not only eigenvalues.

2. The Ritz-value gate can skip the corruption it is meant to detect

Comment on the w_host gate:

The assumption that duplicate Ritz vectors require near-equal Ritz values only holds when V is already orthonormal—in which case V.T @ s cannot produce duplicates at all. In the corrupted state this function is intended to repair, V.T is not an isometry and can map Ritz vectors belonging to distinct w values to the same physical direction.

Minimal algebraic example: take V = [[1,0],[1,0]], s = I, and distinct w = [1,2]. The columns of V.T @ s are duplicates, but this gate skips the Gram check because the Ritz values are well separated.

Could _lanczos_checked return whether its orthogonality/repair path fired and use that as the cheap gate here? Otherwise the Gram check needs to be unconditional when correctness depends on it.

3. Test gap

Comment on [test_dominant_eigenvalue_hidden_from_v0](https://github.com/cupy/cupy/blob/b4ad9547f8f5321de89bfb00907c0ab8bc12a874/tests/cupyx_tests/scipy_tests/sparse_tests/test_linalg.py#L389-L401):

This test never verifies that the hidden dominant eigenvalue was discovered: returning six copies of 1 satisfies both current assertions. Please assert that max(w) is close to big, or compare against [1, 1, 1, 1, 1, big]. That directly exercises the late-discovery behavior the test and commit message claim to cover.

The first finding is the important blocker; the second explains why the current tests may not reach it reliably.

If helpful, I can set up “Re-review PR 10098 on updates” so new commits are checked automatically.

I found two blocking correctness issues in [CuPy PR #10098](#10098). I would not approve it yet.

• P1 — locked-block repair can leave the basis non-orthogonal. In [_eigen.py:634](https://github.com/cupy/cupy/blob/b4ad9547f8f5321de89bfb00907c0ab8bc12a874/cupyx/scipy/sparse/linalg/_eigen.py#L634), _repair_locked replaces V[j] using _restart_ortho(V, j, ...), which projects only against V[:j]. It does not constrain the replacement against V[j+1:k]. Thus the repair can introduce a new duplicate with a later locked row, while the Gram matrix was computed before the replacement and will not catch it. The replacement must be orthogonalized against every other locked row, or the Gram check must be recomputed until the block passes.

• P1 — replacing a locked Ritz vector leaves beta_k stale. The same repair changes V[j] and alpha[j], but the caller retains the old beta_k coupling vector. That vector is subsequently used both to deflate u and as the arrowhead in the next Ritz solve. Its j entry describes the old Ritz vector, not the fresh replacement, so the projected matrix and recurrence are inconsistent. Have _repair_locked report changed indices and recompute or clear the corresponding beta_k entries.

All reported CI checks are green, but the added tests do not exercise either of these repair-state transitions. I found no security concern.

I'll pull the PR details and diff from GitHub.The API is rate-limited; I'll grab the PR page and diff directly instead.The comment extraction mislabeled some authors; let me pull the PR description properly and then get the actual diff from the Files tab.Got the full description. Now the diff itself.The Files tab is client-rendered, so I'll clone the branch with a blob filter and diff it against the base directly.PR is open, approved by leofang with auto-merge enabled, 11 commits at head b4ad954. Now fetching the actual code changes via git.Now the full post-PR file for context, then the tests.Let me also pull the inline review threads so I don't rehash what's already been settled.That's a real finding. Let me also check the complex-dtype path for the sparse bound (abs(b) is real-valued but ones is complex).I've read the full diff at head b4ad954 (2 files, +746/−13), the description, and all six review rounds with @leofang. The PR is currently approved with auto-merge armed and CI running. Overall the design is sound and well-argued, but I found one real bug the review missed, plus a soundness gap and a coverage hole.

1. The "provable" divergence guard is silently disabled for sparse inputs (and all of svds).
_norm_upper_bound does cupy.abs(b) on a csr_matrix. cupy.abs is a bare ufunc; its argument preprocessing only accepts ndarrays / __cuda_array_interface__ objects and wraps anything else in CScalar, which raises TypeError: Unsupported type for a sparse matrix. The blanket except Exception: return None swallows that, so norm_bound is None for every sparse a — the primary input type — and the only remaining check is isfinite(anorm). The new TestEigshLateNormDiscovery uses .tocsr(), so it passes trivially and can't detect this. Fix is one token: use the builtin abs(b) (dispatches to _data_matrix.__abs__); csrmv already handles the real/complex dtype mix via _cast_common_type. I'd also drop or narrow the try/except (the # pragma: no cover is exactly how this hid) and add a unit test that _norm_upper_bound(csr) == abs(dense).sum(1).max().

Relatedly, svds passes a _ProductLinearOperator (aH @ a), so it gets no bound even after that fix. ‖AᴴA‖₂ = ‖A‖₂² ≤ ‖A‖₁·‖A‖∞ is one pass over a and provable; svds could compute it and hand it to eigsh through a private kwarg. #8009/#5000 are the inputs this guard exists for.

2. _repair_locked rests on an incorrect premise, and its repair is incomplete if it ever fires.
The comment says a degenerate Ritz basis lets "several columns of s name the same physical direction." numpy.linalg.eigh returns an orthonormal eigenvector matrix regardless of multiplicity, so V[:k] = (V.T @ s).T is orthonormal iff V was. The author's own final measurement (V[:k] off-diagonal 4.3e-16 on the ghost input) confirms this — the earlier 4.0e-01 was entirely V[k]. So the path only fires when V already lost orthogonality without the sweep-boundary gate catching it, and in that case _restart_ortho(V, j, …) orthogonalizes the new row only against V[:j], not V[j+1:k]; worst is computed once so later rows aren't rechecked; and the new V[j] is not an eigenvector, so ⟨V[i], A V[j]⟩ ≠ 0 and beta_k[j] is stale — the "each locked row is its own 1×1 block" claim doesn't hold. I'd either remove it (it costs a Gram + host sync per restart on any spectrum with legitimate repeated eigenvalues, which are common) or make it fail closed like the other three sites. Follow-up material given the approval, but the comment should not ship as written.

3. Nits

• test_clustered_large_k, _gaussian, test_null_space_start, and the two semidefinite tests ignore the class's k/return_eigenvectors params, so each runs 6× per dtype redundantly. Move them beside TestEigshTinyN.
• test_dominant_eigenvalue_hidden_from_v0 only asserts "no ghost"; the description claims the hidden eigenvalue is found. Either assert w.max() ≈ big or soften the claim.
• _lanczos_checked copies alpha/beta to host, then _eigsh_solve_ritz copies them again. Cheap, but the host copies could be threaded through.
• _norm_upper_bound ignores dia/other formats; isinstance(a, spmatrix) + tocsr() would cover them.

Given #1 is a silent no-op in a fail-closed check that the review explicitly asked for, I'd pull it into this PR before auto-merge fires rather than defer.

Reviewed it from the actual head (b4ad954, 11 commits, +746/−13 across cupyx/scipy/sparse/linalg/_eigen.py and the sparse linalg tests). Status right now: open, leofang approved and the approval was dismissed by the last commit, /test mini running.

What it does: guards the Lanczos normalization V[i+1] = u/beta[i] against lucky breakdown. The inner loop stays sync-free (_kernel_normalize emits a zero column instead of 0/0), and detection/repair happens at the sweep boundary in _lanczos_checked, plus a locked-block Gram check in _repair_locked and a deterministic reseed in _restart_ortho. The failure class is real and the regression tests are good.

Two things I'd want answered before merge:

1. The Gram branch can zero a large beta[p]. The whole design rests on "decoupling at beta[p] perturbs the spectrum by at most beta[p], so a 64·eps·‖A‖ threshold is harmless." That holds for the hits path. It doesn't hold for the orthogonality path: p = j-1 comes from the Gram row that failed, and beta[p] = 0 is then applied with no magnitude check on beta[p]. The comment asserts that a row losing orthogonality implies its predecessor was a collapsed residual, but the textbook cause of orthogonality loss in Lanczos is Ritz convergence (Paige), which happens at perfectly ordinary beta. The entry gate only requires some small beta in the range, not that the Gram-identified position is one. A cheap fix: only accept the Gram-derived q when beta_h[q] <= ortho_rtol * anorm, and fail closed otherwise. Related question: _lanczos_fast already does full reorthogonalization every step, so the residual loss was coming from a rank-deficient V. Now that _repair_locked fixes the duplicate locked rows at the source, is the sweep-boundary Gram branch still earning its place, or can it be narrowed to a fail-closed assertion?

2. _repair_locked can fabricate a converged eigenpair. When a duplicate locked row is reseeded, alpha[j] becomes the Rayleigh quotient of an arbitrary fresh direction. But beta[:k] = 0, so that entry is a decoupled 1×1 block: s[-1, j] = 0, residual identically zero, and the next Ritz solve can select it as a converged eigenvalue. It's a legitimate Rayleigh quotient in [λmin, λmax], so it won't look like garbage. Dropping the duplicated pair from the locked set for that restart seems safer than substituting a fabricated one.

Smaller, follow-up-able:

• The strong divergence guard never fires for svds. svds always calls aslinearoperator and passes aH @ a, so _norm_upper_bound returns None and only the isfinite check applies. That's exactly the truncated SVD returns null singular values. #8009 path. A _ProductLinearOperator could be bounded from its factors.
• Two redundant D2H copies per sweep. _lanczos_checked does asnumpy(alpha) + asnumpy(beta), then _eigsh_solve_ritz immediately does both again. Returning them halves the per-sweep sync count, which is the fixed cost your own ncv sweep (1.8% → 0.1% as ncv grows) fingered as the source of the healthy-path tax.
• _norm_upper_bound: the dense branch materializes a full |A|, an n² temporary on top of an already-large input. cupy.abs(b).sum(axis=1).max() works for both branches and avoids the real-abs @ complex-ones cast in the sparse path. Also except Exception swallows OutOfMemoryError, silently disabling the guard.
• svds docstring doesn't mention the new RuntimeError, though it propagates. Going from "returns NaN" to "raises" is a visible behaviour change for both entry points.
• _repair_locked's matvec isn't counted in the returned work, so maxiter is slightly under-charged.

On cost: a permanent 2–6% on every healthy eigsh/svds call is a real tax on a core solver, and the deferred early-exit-on-exhausted-Krylov-space is the thing that would make the degenerate path sane rather than 20x. Both seem acceptable given the alternative is silent NaN, but the early-exit follow-up is worth filing now so it doesn't get lost.