Fix elementwise variance accumulation in lsqr (closes #639)

[gaoflow]

Problem

lsqr's optional var output estimates the diagonal of $(A^H A)^{-1}$ — i.e. one variance per unknown. The accumulation step used

self.var = self.var + to_numpy_conditional(self.var, self.ncp.dot(self.dk, self.dk))

but ncp.dot(self.dk, self.dk) is the scalar sum of squares of dk, so the same number is added to every entry of var. The result is that all returned variances are identical (and wrong), as reported in #639.

Reproduction (from the issue)

import numpy as np
from scipy.sparse.linalg import aslinearoperator, lsqr as sp_lsqr
from pylops import LinearOperator, lsqr

np.random.seed(0)
A = np.random.randn(20, 10); b = np.random.rand(20)
Aop = LinearOperator(aslinearoperator(A))

var_pylops = lsqr(Aop, b, atol=1e-8, btol=1e-8, calc_var=True, niter=10_000)[9]
var_scipy  = sp_lsqr(aslinearoperator(A), b, atol=1e-8, btol=1e-8, calc_var=True, iter_lim=10_000)[9]
var_dense  = np.diag(np.linalg.inv(A.T @ A))

Before this PR var_pylops is a constant vector ([1.0598, 1.0598, ...]), while var_scipy and var_dense agree.

Fix

Accumulate the elementwise square instead:

self.var = self.var + to_numpy_conditional(self.var, self.ncp.square(self.dk))

This matches the recurrence in Paige & Saunders (1982, p.53) and scipy's lsqr. ncp.square is available across the NumPy / CuPy backends used elsewhere in this method, and the solution path is untouched — only the var output changes. Credit to @IruNikZe for the diagnosis in #639.

Note: like scipy, var is an iterative estimate of $\mathrm{diag}((A^H A)^{-1})$ and converges as iterations proceed; the identity-operator/one-iteration edge case raised later in the issue is inherent to the LSQR estimator (scipy behaves the same) and is not addressed here.

Tests

Added test_lsqr_calc_var (overdetermined real, full column rank) asserting the variances are not all identical (guards against the regression) and match both scipy's lsqr and the dense diag(inv(A^H A)).

Verified locally: pytest pytests/test_solver.py (98 passed), ruff check clean, flake8 (max-line 88) clean on changed files.

[codacy-production]

Up to standards ✅

🟢 Issues 0 issues

Results:
0 new issues

View in Codacy

🟢 Metrics 0 complexity · 0 duplication

Metric	Results
Complexity	0
Duplication	0

View in Codacy

_{NEW Get contextual insights on your PRs based on Codacy's metrics, along with PR and Jira context, without leaving GitHub. Enable AI reviewer
TIP This summary will be updated as you push new changes.}

[mrava87]

@IruNikZe I was wondering if you got anywhere with this in scipy (and pylops), and if you will be willing to review this PR?

[MothNik]

@mrava87 I didn't get anywhere yet, unfortunately. Seems the variances converge (way) slower than the solution and there is not so much that can be done about that.

@gaoflow @mrava87 I'm however happy to review this PR this week.

[mrava87]

@MothNik thanks a lot for the feedback and willingness to review this PR.

What you say makes somehow sense to me, as usually also in VI type of solvers the mean always converges faster than the std, so not surprised this may be also the case for a solver like LSQR

[MothNik]

Please consider also updating the Changelog.md (see here) to reflect your changes.

[gaoflow]

Thanks @MothNik, all addressed in the latest commit:

• Condensed the inline comment in LSQR to a single line.
• Dropped the dense diag(inv(A^H A)) comparison — you're right that LSQR isn't guaranteed to match it; the test now compares only against scipy's LSQR as the reference.
• Good catch on niter: switched to 2 * nx (scipy's default) and added a note that it's only a ceiling since atol/btol stop earlier.
• Moved the #639 reference out of the docstring into the specific "constant vector" comment.
• Added a CHANGELOG.md entry under 2.7.0.

[gaoflow]

Pushed 515e51a to fix the CuPy CI failure: test_lsqr_calc_var compares against scipy's lsqr as the reference, but scipy can't implicitly convert cupy arrays (TypeError: Implicit conversion to a NumPy array is not allowed). Guarded it with the same TEST_CUPY_PYLOPS skipif that test_lsqr_pylops_scipy already uses for the same reason. Passes on the numpy backend; skipped under CuPy.