From d98c4539ea627f47ccd0edfe2187a1cbac94dcbe Mon Sep 17 00:00:00 2001 From: Vincent Gao Date: Fri, 29 May 2026 11:54:38 +0200 Subject: [PATCH 1/4] Fix elementwise variance accumulation in lsqr (closes #639) lsqr's optional `var` output estimates the diagonal of (A^H A)^-1, one variance per unknown. The update accumulated `ncp.dot(self.dk, self.dk)`, which is the scalar sum of squares of `dk`, so the same value was added to every entry of `var` and all returned variances came out identical (and wrong). Accumulate the elementwise square `ncp.square(self.dk)` instead, matching Paige & Saunders (1982, p.53) and scipy's lsqr. The solution path is unchanged; only the `var` output is corrected. Add test_lsqr_calc_var asserting the variances are not all identical and match both scipy's lsqr and the dense diag(inv(A^H A)) on full-column-rank systems. --- pylops/optimization/cls_basic.py | 7 ++++++- pytests/test_solver.py | 27 +++++++++++++++++++++++++++ 2 files changed, 33 insertions(+), 1 deletion(-) diff --git a/pylops/optimization/cls_basic.py b/pylops/optimization/cls_basic.py index 4f47565ff..0996a9f81 100644 --- a/pylops/optimization/cls_basic.py +++ b/pylops/optimization/cls_basic.py @@ -1248,8 +1248,13 @@ def step(self, x: NDArray, show: bool = False) -> NDArray: self.ncp.add(self.v, self.w, out=self.w) self.ddnorm = self.ddnorm + self.ncp.linalg.norm(self.dk) ** 2.0 if self.calc_var: + # ``var`` estimates the diagonal of (A^H A)^-1, so it must accumulate + # the *element-wise* square of ``dk`` (one variance per unknown), as in + # Paige & Saunders (1982) eq. on p.53 and scipy's lsqr. Using + # ``dot(dk, dk)`` instead adds the same scalar (sum of squares) to every + # entry, yielding identical, wrong variances (issue #639). self.var = self.var + to_numpy_conditional( - self.var, self.ncp.dot(self.dk, self.dk) + self.var, self.ncp.square(self.dk) ) # use a plane rotation on the right to eliminate the diff --git a/pytests/test_solver.py b/pytests/test_solver.py index 13a2ca0cd..1e73cd2b0 100644 --- a/pytests/test_solver.py +++ b/pytests/test_solver.py @@ -377,6 +377,33 @@ def test_lsqr_pylops_scipy(par): assert_array_almost_equal(xinv_sp, x, decimal=4) +@pytest.mark.parametrize("par", [(par3), (par4)]) +def test_lsqr_calc_var(par): + """LSQR ``var`` estimates the diagonal of (A^H A)^-1 element-wise (issue #639). + + Each unknown must get its own variance; a regression here would make all + entries identical (the previous ``dot(dk, dk)`` bug). Compare against scipy's + LSQR and the dense ``diag(inv(A^H A))`` on a full-column-rank system. + """ + np.random.seed(10) + + A = np.random.normal(0, 1, (par["ny"], par["nx"])) + Aop = MatrixMult(A, dtype=par["dtype"]) + y = Aop * (np.ones(par["nx"])) + + niter = 10 * par["nx"] + var = lsqr(Aop, y, x0=None, niter=niter, atol=1e-8, btol=1e-8)[9] + var_sp = sp_lsqr( + Aop, y, iter_lim=niter, atol=1e-8, btol=1e-8, calc_var=True + )[9] + var_dense = np.diag(np.linalg.inv(np.conj(A.T) @ A)) + + # not all identical (the bug produced a constant vector) + assert not np.allclose(var, var[0]) + assert_array_almost_equal(var, var_sp, decimal=6) + assert_array_almost_equal(var, var_dense, decimal=6) + + @pytest.mark.parametrize( "par", [(par1), (par2), (par3), (par4), (par1j), (par2j), (par3j), (par3j)] ) From 80b3eeaed12f77e10ca267ad432dd172adb2ffc6 Mon Sep 17 00:00:00 2001 From: gaoflow Date: Thu, 4 Jun 2026 23:24:58 +0200 Subject: [PATCH 2/4] Address review: trim comment, scipy-only test comparison, changelog - Condense the inline comment in LSQR to a single line. - Drop the dense diag(inv(A^H A)) comparison from the test and keep scipy's LSQR as the reference; LSQR is not guaranteed to match the dense variance. - Use niter = 2 * nx (scipy's default) and note that niter is only a ceiling. - Move the issue reference out of the docstring into the targeted comment. - Add a CHANGELOG.md entry. --- CHANGELOG.md | 3 +++ pylops/optimization/cls_basic.py | 6 +----- pytests/test_solver.py | 19 ++++++++----------- 3 files changed, 12 insertions(+), 16 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 91a14ee8c..8f8ed5194 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,6 +17,9 @@ Changelog and `pylops.waveeqprocessing.Marchenko` * Improved typing annotations across all operators (and enforced use of Literal for parameters with multiple options) +* Fixed element-wise variance accumulation in `pylops.optimization.cls_basic.LSQR` + and `pylops.lsqr`, which previously assigned the same value to every entry of + `var` (#639) # 2.6.0 * Added `pylops.medical` module and `pylops.medical.CT2D` operator diff --git a/pylops/optimization/cls_basic.py b/pylops/optimization/cls_basic.py index 0996a9f81..e166fc28a 100644 --- a/pylops/optimization/cls_basic.py +++ b/pylops/optimization/cls_basic.py @@ -1248,11 +1248,7 @@ def step(self, x: NDArray, show: bool = False) -> NDArray: self.ncp.add(self.v, self.w, out=self.w) self.ddnorm = self.ddnorm + self.ncp.linalg.norm(self.dk) ** 2.0 if self.calc_var: - # ``var`` estimates the diagonal of (A^H A)^-1, so it must accumulate - # the *element-wise* square of ``dk`` (one variance per unknown), as in - # Paige & Saunders (1982) eq. on p.53 and scipy's lsqr. Using - # ``dot(dk, dk)`` instead adds the same scalar (sum of squares) to every - # entry, yielding identical, wrong variances (issue #639). + # accumulate the element-wise square of dk (one variance per unknown) self.var = self.var + to_numpy_conditional( self.var, self.ncp.square(self.dk) ) diff --git a/pytests/test_solver.py b/pytests/test_solver.py index 1e73cd2b0..0ba0a7e2a 100644 --- a/pytests/test_solver.py +++ b/pytests/test_solver.py @@ -379,11 +379,10 @@ def test_lsqr_pylops_scipy(par): @pytest.mark.parametrize("par", [(par3), (par4)]) def test_lsqr_calc_var(par): - """LSQR ``var`` estimates the diagonal of (A^H A)^-1 element-wise (issue #639). + """LSQR ``var`` estimates the diagonal of (A^H A)^-1 element-wise. - Each unknown must get its own variance; a regression here would make all - entries identical (the previous ``dot(dk, dk)`` bug). Compare against scipy's - LSQR and the dense ``diag(inv(A^H A))`` on a full-column-rank system. + Each unknown must get its own variance; a regression would make all entries + identical. scipy's LSQR is used as the reference implementation. """ np.random.seed(10) @@ -391,17 +390,15 @@ def test_lsqr_calc_var(par): Aop = MatrixMult(A, dtype=par["dtype"]) y = Aop * (np.ones(par["nx"])) - niter = 10 * par["nx"] + # niter is only a ceiling (atol/btol stop earlier); use scipy's default of + # 2 * nx so both implementations run for the same number of iterations. + niter = 2 * par["nx"] var = lsqr(Aop, y, x0=None, niter=niter, atol=1e-8, btol=1e-8)[9] - var_sp = sp_lsqr( - Aop, y, iter_lim=niter, atol=1e-8, btol=1e-8, calc_var=True - )[9] - var_dense = np.diag(np.linalg.inv(np.conj(A.T) @ A)) + var_sp = sp_lsqr(Aop, y, iter_lim=niter, atol=1e-8, btol=1e-8, calc_var=True)[9] - # not all identical (the bug produced a constant vector) + # the dot(dk, dk) bug produced a constant vector instead (issue #639) assert not np.allclose(var, var[0]) assert_array_almost_equal(var, var_sp, decimal=6) - assert_array_almost_equal(var, var_dense, decimal=6) @pytest.mark.parametrize( From 515e51a765e20da3de414b7afa1173fa025bbbb8 Mon Sep 17 00:00:00 2001 From: gaoflow Date: Fri, 5 Jun 2026 00:01:05 +0200 Subject: [PATCH 3/4] Skip test_lsqr_calc_var under CuPy backend The test compares against scipy's lsqr as reference; scipy cannot implicitly convert cupy arrays (TypeError), so guard it with the same TEST_CUPY_PYLOPS skipif used by test_lsqr_pylops_scipy. --- pytests/test_solver.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pytests/test_solver.py b/pytests/test_solver.py index 0ba0a7e2a..3804faf5a 100644 --- a/pytests/test_solver.py +++ b/pytests/test_solver.py @@ -377,6 +377,9 @@ def test_lsqr_pylops_scipy(par): assert_array_almost_equal(xinv_sp, x, decimal=4) +@pytest.mark.skipif( + int(os.environ.get("TEST_CUPY_PYLOPS", 0)) == 1, reason="Not CuPy enabled" +) @pytest.mark.parametrize("par", [(par3), (par4)]) def test_lsqr_calc_var(par): """LSQR ``var`` estimates the diagonal of (A^H A)^-1 element-wise. From 5938feffcbc68c71a6ffab723c51aaefd3796e58 Mon Sep 17 00:00:00 2001 From: gaoflow Date: Sun, 7 Jun 2026 09:35:27 +0200 Subject: [PATCH 4/4] Address review: move changelog to 2.8.0, simplify test docstring/comments --- CHANGELOG.md | 8 +++++--- pytests/test_solver.py | 8 ++------ 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8f8ed5194..0356a57b9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,11 @@ Changelog ========= +# 2.8.0 +* Fixed element-wise variance accumulation in `pylops.optimization.cls_basic.LSQR` + and `pylops.lsqr`, which previously assigned the same value to every entry of + `var` (#639) + # 2.7.0 * Added cubic spline interpolation operator via `pylops.signalprocessing.interpspline.InterpCubicSpline` (also interfaceable via @@ -17,9 +22,6 @@ Changelog and `pylops.waveeqprocessing.Marchenko` * Improved typing annotations across all operators (and enforced use of Literal for parameters with multiple options) -* Fixed element-wise variance accumulation in `pylops.optimization.cls_basic.LSQR` - and `pylops.lsqr`, which previously assigned the same value to every entry of - `var` (#639) # 2.6.0 * Added `pylops.medical` module and `pylops.medical.CT2D` operator diff --git a/pytests/test_solver.py b/pytests/test_solver.py index 3804faf5a..96325db76 100644 --- a/pytests/test_solver.py +++ b/pytests/test_solver.py @@ -382,11 +382,8 @@ def test_lsqr_pylops_scipy(par): ) @pytest.mark.parametrize("par", [(par3), (par4)]) def test_lsqr_calc_var(par): - """LSQR ``var`` estimates the diagonal of (A^H A)^-1 element-wise. - - Each unknown must get its own variance; a regression would make all entries - identical. scipy's LSQR is used as the reference implementation. - """ + """Compare PyLops and scipy LSQR variance computation for the diagonal of + (A^H A)^-1 (issue #639).""" np.random.seed(10) A = np.random.normal(0, 1, (par["ny"], par["nx"])) @@ -399,7 +396,6 @@ def test_lsqr_calc_var(par): var = lsqr(Aop, y, x0=None, niter=niter, atol=1e-8, btol=1e-8)[9] var_sp = sp_lsqr(Aop, y, iter_lim=niter, atol=1e-8, btol=1e-8, calc_var=True)[9] - # the dot(dk, dk) bug produced a constant vector instead (issue #639) assert not np.allclose(var, var[0]) assert_array_almost_equal(var, var_sp, decimal=6)