Repository navigation
LOBPCG can stop early with returned residual above tol after soft locking #10402
Description
cupyx.scipy.sparse.linalg.lobpcg in CuPy 14.2.0 can stop before maxiter while a returned eigenpair has an absolute residual above tol. This also appears in the last residual-history entry, so it is not solely a final-postprocessing discrepancy.
A 128-by-128 positive diagonal matrix with repeated eigenvalues reproduces it without a custom operator, preconditioner, generalized mass, or constraints:
import numpy as np
import cupy as cp
import cupyx.scipy.sparse as sparse
from cupyx.scipy.sparse.linalg import lobpcg
diagonal = np.repeat(np.arange(1.0, 33.0), 4)
A = sparse.diags(cp.asarray(diagonal), format="csr")
X = cp.asarray(np.random.default_rng(42).normal(size=(128, 6)))
w, v, history = lobpcg(
A, X, largest=False, tol=1e-10, maxiter=200,
retResidualNormsHistory=True,
)
w, v = w.get(), v.get()
history = cp.stack(history).get()
residual = np.linalg.norm(diagonal[:, None] * v - v * w, axis=0)
print(len(history))
print(history[-1])
print(history.min(axis=0))
print(residual)On an NVIDIA L40, Python 3.10, NumPy 2.2.6, cupy-cuda12x 14.2.0, CUDA runtime 12.9, driver 570.158.01, three restored-start trials give the same result:
history length: 76 # maxiter=200
last history: [6.463312e-11, 3.060786e-11, 1.401550e-10,
1.028023e-11, 7.569174e-11, 9.249482e-11]
minimum by history column:
[8.522374e-12, 1.202785e-11, 1.079567e-11,
8.319244e-12, 7.569174e-11, 9.249482e-11]
fresh CPU residual:
[6.463298e-11, 3.060790e-11, 1.401545e-10,
1.028019e-11, 7.569152e-11, 9.249490e-11]
Eigenvalues and orthogonality are accurate (maximum errors approximately 3.55e-15 and 8.88e-16); the discrepancy concerns the requested stopping criterion. A disconnected cycle-Laplacian example also returns residual 3.77e-10 for tolerance 1e-10 before its iteration limit.
The implementation updates activeMask = activeMask & ii, where ii identifies residuals above tolerance. Thus a column cannot reactivate after its residual grows, even though all Ritz vectors continue to rotate. This seems relevant, but I have not established that simply allowing reactivation is a safe fix across the full solver suite. Minima of history columns are not fixed physical mode identities within a repeated eigenspace.
Is this intended soft-locking behavior? It would be helpful to ensure all returned pairs satisfy the requested criterion before early termination, or warn/document when they do not. The API documentation describes tol as the stopping criterion.
Activity
seberg commentedon Oct 22, 2026
This (and similar issues), CuPy follows SciPy (and this code is probably largely a copy-paste from there), so some of these questions seem like they may need to be answered there. Ping @ev-br in case you have some quick input, if just to encourage opening an umbrella issue for these in SciPy.
ev-br commentedon Oct 22, 2026
I've to admit I don't know if this is a desired behaviour or not. (And agreed it's better to fix it in both SciPy and CuPy if it's not).
Could you weight in on whether it's a bug or a feature @lobpcg?
lobpcg commentedon Oct 22, 2026
It's not desired, but expected. Soft-locking minimizes a different function, the Raleigh quotient, so the residual may indeed increase but it doesn't mean that the output is actually less accurate. The Raleigh quotient minimizing is more relevant to the true accuracy compared to the residual. In other words, using this residual for stopping criteria is a simple traditional choice but should not be taken too seriously. This evidently creates this software maintenance trouble reported here. The easiest fix is to declare this as a feature and document accordingly.
Let me run a few examples to double check this and suggest a good comment to add.