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?