jlec 14/11/14 20:34:09 Added: scipy-0.14.0-lsqr-backport.patch Log: sci-libs/scipy: Backport fix for broken lsqr (Portage version: 2.2.14/cvs/Linux x86_64, signed Manifest commit with key B9D4F231BD1558AB!)
Revision Changes Path 1.1 sci-libs/scipy/files/scipy-0.14.0-lsqr-backport.patch file : http://sources.gentoo.org/viewvc.cgi/gentoo-x86/sci-libs/scipy/files/scipy-0.14.0-lsqr-backport.patch?rev=1.1&view=markup plain: http://sources.gentoo.org/viewvc.cgi/gentoo-x86/sci-libs/scipy/files/scipy-0.14.0-lsqr-backport.patch?rev=1.1&content-type=text/plain Index: scipy-0.14.0-lsqr-backport.patch =================================================================== >From 785e3961a685a351bd9d1e8eec0a62035cc4c6aa Mon Sep 17 00:00:00 2001 From: Olivier Grisel <[email protected]> Date: Thu, 13 Nov 2014 13:39:35 +0100 Subject: [PATCH] BUG: ZeroDivisionError in scipy.sparse.linalg.lsqr --- scipy/sparse/linalg/isolve/lsqr.py | 6 ++-- scipy/sparse/linalg/isolve/tests/test_lsqr.py | 51 +++++++++++++++++++++------ 2 files changed, 45 insertions(+), 12 deletions(-) diff --git a/scipy/sparse/linalg/isolve/lsqr.py b/scipy/sparse/linalg/isolve/lsqr.py index c4a6c02..de5a6b9 100644 --- a/scipy/sparse/linalg/isolve/lsqr.py +++ b/scipy/sparse/linalg/isolve/lsqr.py @@ -57,6 +57,8 @@ from math import sqrt from scipy.sparse.linalg.interface import aslinearoperator +eps = np.finfo(np.float64).eps + def _sym_ortho(a, b): """ @@ -432,8 +434,8 @@ def lsqr(A, b, damp=0.0, atol=1e-8, btol=1e-8, conlim=1e8, # Now use these norms to estimate certain other quantities, # some of which will be small near a solution. test1 = rnorm / bnorm - test2 = arnorm / (anorm * rnorm) - test3 = 1 / acond + test2 = arnorm / (anorm * rnorm + eps) + test3 = 1 / (acond + eps) t1 = test1 / (1 + anorm * xnorm / bnorm) rtol = btol + atol * anorm * xnorm / bnorm diff --git a/scipy/sparse/linalg/isolve/tests/test_lsqr.py b/scipy/sparse/linalg/isolve/tests/test_lsqr.py index f378cca..626288b 100644 --- a/scipy/sparse/linalg/isolve/tests/test_lsqr.py +++ b/scipy/sparse/linalg/isolve/tests/test_lsqr.py @@ -1,7 +1,7 @@ from __future__ import division, print_function, absolute_import import numpy as np -from numpy.testing import assert_ +from numpy.testing import assert_, assert_equal, assert_array_almost_equal from scipy.lib.six import xrange import scipy.sparse @@ -34,6 +34,46 @@ def test_basic(): xo = X[0] assert_(norm(svx - xo) < 1e-5) + +def test_gh_2466(): + row = np.array([0, 0]) + col = np.array([0, 1]) + val = np.array([1, -1]) + A = scipy.sparse.coo_matrix((val, (row, col)), shape=(1, 2)) + b = np.asarray([4]) + lsqr(A, b) + + +def test_well_conditioned_problems(): + # Test that sparse the lsqr solver returns the right solution + # on various problems with different random seeds. + # This is a non-regression test for a potential ZeroDivisionError + # raised when computing the `test2` & `test3` convergence conditions. + n = 10 + A_sparse = scipy.sparse.eye(n, n) + A_dense = A_sparse.toarray() + + with np.errstate(invalid='raise'): + for seed in range(30): + rng = np.random.RandomState(seed + 10) + beta = rng.rand(n) + beta[beta == 0] = 0.00001 # ensure that all the betas are not null + b = A_sparse * beta[:, np.newaxis] + output = lsqr(A_sparse, b, show=show) + + # Check that the termination condition corresponds to an approximate + # solution to Ax = b + assert_equal(output[1], 1) + solution = output[0] + + # Check that we recover the ground truth solution + assert_array_almost_equal(solution, beta) + + # Sanity check: compare to the dense array solver + reference_solution = np.linalg.solve(A_dense, b).ravel() + assert_array_almost_equal(solution, reference_solution) + + if __name__ == "__main__": svx = np.linalg.solve(G, b) @@ -64,12 +104,3 @@ def test_basic(): print("") print(" || x_{direct} - x_{LSQR}|| %9.4e " % norm(svx-xo)) print("") - - -def test_gh_2466(): - row = np.array([0, 0]) - col = np.array([0, 1]) - val = np.array([1, -1]) - A = scipy.sparse.coo_matrix((val, (row, col)), shape=(1, 2)) - b = np.asarray([4]) - lsqr(A, b)
