Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Loading…

ENH: Add LSMR iterative solver. #103

Merged
merged 16 commits into from

4 participants

@stefanv
Owner

LSMR is an iterative method for solving (sparse) linear and linear least-squares systems. The method is based on the Golub-Kahan bi-diagonalization process. It is analytically equivalent to the standard method of MINRES applied to the normal equation. Compared to LSQR (already included in scipy), it is safer to terminate LSMR early.

Also see http://www.stanford.edu/~clfong/lsmr.html

@davidfong

LSMR algorithm has been published at
http://dx.doi.org/10.1137/10079687X

@rgommers
Owner

Looks good based on a quick browse of the code. Here are some things I noticed:

  • the while loop should have a maxiter break condition or be changed into a for loop, the risk of hanging is too high otherwise (see the recent fmin_bfgs and changes)
  • the eight if statements at the end in the if show block can be merged to a single if statement by combining conditions with ors.
  • some PEP8ifying may be in order: the body of if-statements on new lines, single spaces around "=" operator.
  • I'm not a big fan of adding more print statements. How about using the logging module instead? Also, the start and end multiple print statements of one-line strings can be merged to a single multi-line string.
  • the indentation of lowerBidiagonalMatrix is incorrect
@rgommers
Owner

On second thought, the print statements are not a big deal. This should be cleaned up all over the codebase at once.

Updating the reference to the published version would be good.

@rgommers
Owner

Are you planning to update this sometime, Stefan? Would be nice to get it in soon, I'd rather not pile up the work to just before the next release.

scipy/sparse/linalg/isolve/lsqr.py
((6 lines not shown))
"""
- Jeffery Kline noted: I added the routine 'SymOrtho' for numerical
+ Stable implementation of Givens rotation.
+
+ Notes
+ -----
+ Jeffery Kline: I added the routine 'SymOrtho' for numerical
@rgommers Owner

Someone's name doesn't belong in a docstring.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers commented on the diff
scipy/sparse/linalg/isolve/tests/test_lsmr.py
@@ -0,0 +1,150 @@
+"""
+Copyright (C) 2010 David Fong and Michael Saunders
+Distributed under the same license as Scipy
@rgommers Owner

Wouldn't it be better to leave out this license statement, then it's automatically covered by the SciPy license.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers
Owner

The function listing in the sparse.linalg docstring should also be updated.

@stefanv
Owner

@rgommers I guess we've missed 0.10, but I'll fix up all the things you mentioned above and update.

@stefanv
Owner
  • The while loop is already limited by itnlim.

  • I left the reference as is--the final publication is pay-only, whereas the pre-print is free.

  • Multi-line strings are annoying in that they mess up the indentation. Maybe you have a good suggestion on how to format those.

  • I think I fixed everything else.

@stefanv
Owner

Rebased. Should be able to merge automatically now. (But not sure how to tell GitHub that I did?)

@pv pv commented on the diff
scipy/sparse/linalg/isolve/lsmr.py
((384 lines not shown))
+ break
+
+ # Print the stopping condition.
+
+ if show:
+ print ' '
+ print 'LSMR finished'
+ print msg[istop]
+ print 'istop =%8g normr =%8.1e' % (istop, normr)
+ print ' normA =%8.1e normAr =%8.1e' % (normA, normar)
+ print 'itn =%8g condA =%8.1e' % (itn, condA)
+ print ' normx =%8.1e' % (normx)
+ print str1, str2
+ print str3, str4
+
+ return x, istop, itn, normr, normar, normA, condA, normx
@pv Owner
pv added a note

Here, and in other places in Scipy, it would be nice to used named tuples (although we need to bundle the implemetation in scipy.lib, as we want to support Pythons < 2.6). However, this can be changed also after this is merged.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
scipy/sparse/linalg/isolve/lsmr.py
((14 lines not shown))
+Systems Optimization Laboratory
+Dept of MS&E, Stanford University.
+
+"""
+
+__all__ = ['lsmr']
+
+from numpy import zeros, infty
+from numpy.linalg import norm
+from math import sqrt
+from scipy.sparse.linalg.interface import aslinearoperator
+
+from lsqr import _sym_ortho
+
+def lsmr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
+ itnlim=None, show=False):
@pv Owner
pv added a note

I think itnlimit is called maxiter in the other iterative routines.

@stefanv Owner
stefanv added a note

Updated.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers
Owner

For multi-line strings with correct indentation in the code: http://docs.python.org/library/textwrap.html#textwrap.dedent

@rgommers
Owner

The reference should be the full one, so including SIAM J. on Sci. Comp., vol. 33, pp. 2950-2971, 2011.. There probably shouldn't be any link in it - we don't do that anywhere, and Google Scholar is quite good at locating a pdf copy.

@stefanv
Owner

We link to Wikipedia pages all over, so I don't see the problem. And are you insisting that I refer to a version of the paper that most users cannot read?

At this stage, it really feels like we're harping on technicalities.

@rgommers
Owner

I'm also fine with leaving the link in if you prefer that. I do think it is necessary to give the full reference to the published paper, both because it's how you properly cite a paper and because it's what our doc standard says.

@rgommers
Owner

For the rest this looks finished, and I'm happy to merge it.

@rgommers rgommers merged commit 1d317ef into scipy:master
@rgommers
Owner

OK, I went ahead and merged this. I'll add a note to the release notes and complete the reference, leaving the link as is.

Thanks again, Stefan and David.

@stefanv
Owner

Erk, sorry--I was meaning to do those last cleanups, still. I was just tired and grumpy the other evening :) I'll make a new PR for that.

@stefanv
Owner

@rgommers Great, looks like you made that last tweak. Thanks for all the reviews and time spent on this.

@rgommers
Owner

No worries.You can make a PR for something else if you want:)

@davidfong

Thanks a lot for getting this in!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
This page is out of date. Refresh to see the latest.
View
2  THANKS.txt
@@ -97,7 +97,7 @@ Chris Jordan-Squire for bug fixes, documentation improvements and
Christoph Gohlke for many bug fixes and help with Windows specific issues.
Jacob Silterra for cwt-based peak finding in scipy.signal.
Denis Laxalde for the unified interface to minimizers in scipy.optimize.
-
+David Fong for the sparse LSMR solver.
Institutions
------------
View
1  scipy/sparse/linalg/__init__.py
@@ -45,6 +45,7 @@
:toctree: generated/
lsqr -- Find the least-squares solution to a sparse linear equation system
+ lsmr -- Find the least-squares solution to a sparse linear equation system
Matrix factorizations
---------------------
View
1  scipy/sparse/linalg/isolve/__init__.py
@@ -5,6 +5,7 @@
from minres import minres
from lgmres import lgmres
from lsqr import lsqr
+from lsmr import lsmr
__all__ = filter(lambda s:not s.startswith('_'),dir())
from numpy.testing import Tester
View
399 scipy/sparse/linalg/isolve/lsmr.py
@@ -0,0 +1,399 @@
+"""
+Copyright (C) 2010 David Fong and Michael Saunders
+
+LSMR uses an iterative method.
+
+07 Jun 2010: Documentation updated
+03 Jun 2010: First release version in Python
+
+David Chin-lung Fong clfong@stanford.edu
+Institute for Computational and Mathematical Engineering
+Stanford University
+
+Michael Saunders saunders@stanford.edu
+Systems Optimization Laboratory
+Dept of MS&E, Stanford University.
+
+"""
+
+__all__ = ['lsmr']
+
+from numpy import zeros, infty
+from numpy.linalg import norm
+from math import sqrt
+from scipy.sparse.linalg.interface import aslinearoperator
+
+from lsqr import _sym_ortho
+
+def lsmr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
+ maxiter=None, show=False):
+ """Iterative solver for least-squares problems.
+
+ lsmr solves the system of linear equations ``Ax = b``. If the system
+ is inconsistent, it solves the least-squares problem ``min ||b - Ax||_2``.
+ A is a rectangular matrix of dimension m-by-n, where all cases are
+ allowed: m = n, m > n, or m < n. B is a vector of length m.
+ The matrix A may be dense or sparse (usually sparse).
+
+ Parameters
+ ----------
+ A : {matrix, sparse matrix, ndarray, LinearOperator}
+ Matrix A in the linear system.
+ b : (m,) ndarray
+ Vector b in the linear system.
+ damp : float
+ Damping factor for regularized least-squares. `lsmr` solves
+ the regularized least-squares problem::
+
+ min ||(b) - ( A )x||
+ ||(0) (damp*I) ||_2
+
+ where damp is a scalar. If damp is None or 0, the system
+ is solved without regularization.
+ atol, btol : float
+ Stopping tolerances. `lsmr` continues iterations until a
+ certain backward error estimate is smaller than some quantity
+ depending on atol and btol. Let ``r = b - Ax`` be the
+ residual vector for the current approximate solution ``x``.
+ If ``Ax = b`` seems to be consistent, ``lsmr`` terminates
+ when ``norm(r) <= atol * norm(A) * norm(x) + btol * norm(b)``.
+ Otherwise, lsmr terminates when ``norm(A^{T} r) <=
+ atol * norm(A) * norm(r)``. If both tolerances are 1.0e-6 (say),
+ the final ``norm(r)`` should be accurate to about 6
+ digits. (The final x will usually have fewer correct digits,
+ depending on ``cond(A)`` and the size of LAMBDA.) If `atol`
+ or `btol` is None, a default value of 1.0e-6 will be used.
+ Ideally, they should be estimates of the relative error in the
+ entries of A and B respectively. For example, if the entries
+ of `A` have 7 correct digits, set atol = 1e-7. This prevents
+ the algorithm from doing unnecessary work beyond the
+ uncertainty of the input data.
+ conlim : float
+ `lsmr` terminates if an estimate of ``cond(A)`` exceeds
+ `conlim`. For compatible systems ``Ax = b``, conlim could be
+ as large as 1.0e+12 (say). For least-squares problems,
+ `conlim` should be less than 1.0e+8. If `conlim` is None, the
+ default value is 1e+8. Maximum precision can be obtained by
+ setting ``atol = btol = conlim = 0``, but the number of
+ iterations may then be excessive.
+ maxiter : int
+ `lsmr` terminates if the number of iterations reaches
+ `maxiter`. The default is ``maxiter = min(m, n)``. For
+ ill-conditioned systems, a larger value of `maxiter` may be
+ needed.
+ show : bool
+ Print iterations logs if ``show=True``.
+
+ Returns
+ -------
+ x : ndarray of float
+ Least-square solution returned.
+ istop : int
+ istop gives the reason for stopping::
+
+ istop = 0 means x=0 is a solution.
+ = 1 means x is an approximate solution to A*x = B,
+ according to atol and btol.
+ = 2 means x approximately solves the least-squares problem
+ according to atol.
+ = 3 means COND(A) seems to be greater than CONLIM.
+ = 4 is the same as 1 with atol = btol = eps (machine
+ precision)
+ = 5 is the same as 2 with atol = eps.
+ = 6 is the same as 3 with CONLIM = 1/eps.
+ = 7 means ITN reached maxiter before the other stopping
+ conditions were satisfied.
+
+ itn : int
+ Number of iterations used.
+ normr : float
+ ``norm(b-Ax)``
+ normar : float
+ ``norm(A^T (b - Ax))``
+ norma : float
+ ``norm(A)``
+ conda : float
+ Condition number of A.
+ normx : float
+ ``norm(x)``
+
+ References
+ ----------
+ .. [1] D. C.-L. Fong and M. A. Saunders,
+ "LSMR: An iterative algorithm for least-square problems",
+ http://arxiv.org/abs/1006.0758
+ .. [2] LSMR Software, http://www.stanford.edu/~clfong/lsmr.html
+
+ """
+
+ A = aslinearoperator(A)
+ b = b.squeeze()
+
+ msg=('The exact solution is x = 0 ',
+ 'Ax - b is small enough, given atol, btol ',
+ 'The least-squares solution is good enough, given atol ',
+ 'The estimate of cond(Abar) has exceeded conlim ',
+ 'Ax - b is small enough for this machine ',
+ 'The least-squares solution is good enough for this machine',
+ 'Cond(Abar) seems to be too large for this machine ',
+ 'The iteration limit has been reached ')
+
+ hdg1 = ' itn x(1) norm r norm A''r'
+ hdg2 = ' compatible LS norm A cond A'
+ pfreq = 20 # print frequency (for repeating the heading)
+ pcount = 0 # print counter
+
+ m, n = A.shape
+
+ # stores the num of singular values
+ minDim = min([m, n])
+
+ if maxiter is None:
+ maxiter = minDim
+
+ if show:
+ print ' '
+ print 'LSMR Least-squares solution of Ax = b\n'
+ print 'The matrix A has %8g rows and %8g cols' % (m, n)
+ print 'damp = %20.14e\n' % (damp)
+ print 'atol = %8.2e conlim = %8.2e\n' % (atol, conlim)
+ print 'btol = %8.2e maxiter = %8g\n' % (btol, maxiter)
+
+ u = b
+ beta = norm(u)
+
+ v = zeros(n)
+ alpha = 0
+
+ if beta > 0:
+ u = (1 / beta) * u
+ v = A.rmatvec(u)
+ alpha = norm(v)
+
+ if alpha > 0:
+ v = (1 / alpha) * v
+
+
+ # Initialize variables for 1st iteration.
+
+ itn = 0
+ zetabar = alpha * beta
+ alphabar = alpha
+ rho = 1
+ rhobar = 1
+ cbar = 1
+ sbar = 0
+
+ h = v.copy()
+ hbar = zeros(n)
+ x = zeros(n)
+
+ # Initialize variables for estimation of ||r||.
+
+ betadd = beta
+ betad = 0
+ rhodold = 1
+ tautildeold = 0
+ thetatilde = 0
+ zeta = 0
+ d = 0
+
+ # Initialize variables for estimation of ||A|| and cond(A)
+
+ normA2 = alpha * alpha
+ maxrbar = 0
+ minrbar = 1e+100
+ normA = sqrt(normA2)
+ condA = 1
+ normx = 0
+
+ # Items for use in stopping rules.
+ normb = beta
+ istop = 0
+ ctol = 0
+ if conlim > 0:
+ ctol = 1 / conlim
+ normr = beta
+
+ # Reverse the order here from the original matlab code because
+ # there was an error on return when arnorm==0
+ normar = alpha * beta
+ if normar == 0:
+ if show:
+ print msg[0]
+ return x, istop, itn, normr, normar, normA, condA, normx
+
+ if show:
+ print ' '
+ print hdg1, hdg2
+ test1 = 1
+ test2 = alpha / beta
+ str1 = '%6g %12.5e' % (itn, x[0])
+ str2 = ' %10.3e %10.3e' % (normr, normar)
+ str3 = ' %8.1e %8.1e' % (test1, test2)
+ print ''.join([str1, str2, str3])
+
+ # Main iteration loop.
+ while itn < maxiter:
+ itn = itn + 1
+
+ # Perform the next step of the bidiagonalization to obtain the
+ # next beta, u, alpha, v. These satisfy the relations
+ # beta*u = a*v - alpha*u,
+ # alpha*v = A'*u - beta*v.
+
+ u = A.matvec(v) - alpha * u
+ beta = norm(u)
+
+ if beta > 0:
+ u = (1 / beta) * u
+ v = A.rmatvec(u) - beta * v
+ alpha = norm(v)
+ if alpha > 0:
+ v = (1 / alpha) * v
+
+ # At this point, beta = beta_{k+1}, alpha = alpha_{k+1}.
+
+ # Construct rotation Qhat_{k,2k+1}.
+
+ chat, shat, alphahat = _sym_ortho(alphabar, damp)
+
+ # Use a plane rotation (Q_i) to turn B_i to R_i
+
+ rhoold = rho
+ c, s, rho = _sym_ortho(alphahat, beta)
+ thetanew = s*alpha
+ alphabar = c*alpha
+
+ # Use a plane rotation (Qbar_i) to turn R_i^T to R_i^bar
+
+ rhobarold = rhobar
+ zetaold = zeta
+ thetabar = sbar * rho
+ rhotemp = cbar * rho
+ cbar, sbar, rhobar = _sym_ortho(cbar * rho, thetanew)
+ zeta = cbar * zetabar
+ zetabar = - sbar * zetabar
+
+ # Update h, h_hat, x.
+
+ hbar = h - (thetabar * rho / (rhoold * rhobarold)) * hbar
+ x = x + (zeta / (rho * rhobar)) * hbar
+ h = v - (thetanew / rho) * h
+
+ # Estimate of ||r||.
+
+ # Apply rotation Qhat_{k,2k+1}.
+ betaacute = chat * betadd
+ betacheck = -shat * betadd
+
+ # Apply rotation Q_{k,k+1}.
+ betahat = c * betaacute
+ betadd = -s * betaacute
+
+ # Apply rotation Qtilde_{k-1}.
+ # betad = betad_{k-1} here.
+
+ thetatildeold = thetatilde
+ ctildeold, stildeold, rhotildeold = _sym_ortho(rhodold, thetabar)
+ thetatilde = stildeold* rhobar
+ rhodold = ctildeold * rhobar
+ betad = - stildeold * betad + ctildeold * betahat
+
+ # betad = betad_k here.
+ # rhodold = rhod_k here.
+
+ tautildeold = (zetaold - thetatildeold * tautildeold) / rhotildeold
+ taud = (zeta - thetatilde * tautildeold) / rhodold
+ d = d + betacheck * betacheck
+ normr = sqrt(d + (betad - taud)**2 + betadd * betadd)
+
+ # Estimate ||A||.
+ normA2 = normA2 + beta * beta
+ normA = sqrt(normA2)
+ normA2 = normA2 + alpha * alpha
+
+ # Estimate cond(A).
+ maxrbar = max(maxrbar, rhobarold)
+ if itn > 1:
+ minrbar= min(minrbar, rhobarold)
+ condA = max(maxrbar, rhotemp) / min(minrbar, rhotemp)
+
+ # Test for convergence.
+
+ # Compute norms for convergence testing.
+ normar = abs(zetabar)
+ normx = norm(x)
+
+ # Now use these norms to estimate certain other quantities,
+ # some of which will be small near a solution.
+
+ test1 = normr / normb
+ if (normA * normr) != 0:
+ test2 = normar / (normA * normr)
+ else:
+ test2 = infty
+ test3 = 1 / condA
+ t1 = test1 / (1 + normA * normx / normb)
+ rtol = btol + atol * normA * normx / normb
+
+ # The following tests guard against extremely small values of
+ # atol, btol or ctol. (The user may have set any or all of
+ # the parameters atol, btol, conlim to 0.)
+ # The effect is equivalent to the normAl tests using
+ # atol = eps, btol = eps, conlim = 1/eps.
+
+ if itn >= maxiter:
+ istop = 7
+ if 1 + test3 <= 1:
+ istop = 6
+ if 1 + test2 <= 1:
+ istop = 5
+ if 1 + t1 <= 1:
+ istop = 4
+
+ # Allow for tolerances set by the user.
+
+ if test3 <= ctol:
+ istop = 3
+ if test2 <= atol:
+ istop = 2
+ if test1 <= rtol:
+ istop = 1
+
+ # See if it is time to print something.
+
+ if show:
+ if (n <= 40) or (itn <= 10) or (itn >= maxiter - 10) or \
+ (itn % 10 == 0) or (test3 <= 1.1 * ctol) or \
+ (test2 <= 1.1 * atol) or (test1 <= 1.1 * rtol) or \
+ (istop != 0):
+
+ if pcount >= pfreq:
+ pcount = 0
+ print ' '
+ print hdg1, hdg2
+ pcount = pcount + 1
+ str1 = '%6g %12.5e' % (itn, x[0])
+ str2 = ' %10.3e %10.3e' % (normr, normar)
+ str3 = ' %8.1e %8.1e' % (test1, test2)
+ str4 = ' %8.1e %8.1e' % (normA, condA)
+ print ''.join([str1, str2, str3, str4])
+
+ if istop > 0:
+ break
+
+ # Print the stopping condition.
+
+ if show:
+ print ' '
+ print 'LSMR finished'
+ print msg[istop]
+ print 'istop =%8g normr =%8.1e' % (istop, normr)
+ print ' normA =%8.1e normAr =%8.1e' % (normA, normar)
+ print 'itn =%8g condA =%8.1e' % (itn, condA)
+ print ' normx =%8.1e' % (normx)
+ print str1, str2
+ print str3, str4
+
+ return x, istop, itn, normr, normar, normA, condA, normx
@pv Owner
pv added a note

Here, and in other places in Scipy, it would be nice to used named tuples (although we need to bundle the implemetation in scipy.lib, as we want to support Pythons < 2.6). However, this can be changed also after this is merged.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
View
45 scipy/sparse/linalg/isolve/lsqr.py
@@ -55,13 +55,16 @@
from math import sqrt
from scipy.sparse.linalg.interface import aslinearoperator
-def _sym_ortho(a,b):
+def _sym_ortho(a, b):
"""
- Jeffery Kline noted: I added the routine 'SymOrtho' for numerical
- stability. This is recommended by S.-C. Choi in [1]_. It removes
- the unpleasant potential of ``1/eps`` in some important places
- (see, for example text following "Compute the next
- plane rotation Qk" in minres_py).
+ Stable implementation of Givens rotation.
+
+ Notes
+ -----
+ The routine 'SymOrtho' was added for numerical stability. This is
+ recommended by S.-C. Choi in [1]_. It removes the unpleasant potential of
+ ``1/eps`` in some important places (see, for example text following
+ "Compute the next plane rotation Qk" in minres_py).
References
----------
@@ -70,34 +73,18 @@ def _sym_ortho(a,b):
http://www.stanford.edu/group/SOL/dissertations/sou-cheng-choi-thesis.pdf
"""
- aa = abs(a)
- ab = abs(b)
- if b == 0.:
- s = 0.
- r = aa
- if aa == 0.:
- c = 1.
- else:
- c = a/aa
- elif a == 0.:
- c = 0.
- s = b / ab
- r = ab
- elif ab >= aa:
- sb = 1
- if b < 0: sb=-1
- tau = a/b
- s = sb * (1 + tau**2)**-0.5
+ if b == 0: return np.sign(a), 0, abs(a)
+ elif a == 0: return 0, np.sign(b), abs(b)
+ elif abs(b) > abs(a):
+ tau = a / b
+ s = np.sign(b) / sqrt(1 + tau * tau)
c = s * tau
r = b / s
- elif aa > ab:
- sa = 1
- if a < 0: sa = -1
+ else:
tau = b / a
- c = sa * (1 + tau**2)**-0.5
+ c = np.sign(a) / sqrt(1+tau*tau)
s = c * tau
r = a / c
-
return c, s, r
View
149 scipy/sparse/linalg/isolve/tests/test_lsmr.py
@@ -0,0 +1,149 @@
+"""
+Copyright (C) 2010 David Fong and Michael Saunders
+Distributed under the same license as Scipy
@rgommers Owner

Wouldn't it be better to leave out this license statement, then it's automatically covered by the SciPy license.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
+
+Testing Code for LSMR.
+
+03 Jun 2010: First version release with lsmr.py
+
+David Chin-lung Fong clfong@stanford.edu
+Institute for Computational and Mathematical Engineering
+Stanford University
+
+Michael Saunders saunders@stanford.edu
+Systems Optimization Laboratory
+Dept of MS&E, Stanford University.
+
+"""
+from numpy import arange, concatenate, identity, zeros, ones, sqrt, \
+ transpose, hstack
+from numpy.linalg import norm
+from numpy.testing import run_module_suite, assert_almost_equal
+
+from scipy.sparse import coo_matrix
+from scipy.sparse.linalg.interface import aslinearoperator
+from scipy.sparse.linalg import lsmr
+
+class TestLSMR:
+ def setUp(self):
+ self.n = 10
+ self.m = 10
+
+ def assertCompatibleSystem(self, A, xtrue):
+ Afun = aslinearoperator(A)
+ b = Afun.matvec(xtrue)
+ x = lsmr(A,b)[0]
+ assert_almost_equal(norm(x - xtrue), 0, 6)
+
+ def testIdentityACase1(self):
+ A = identity(self.n)
+ xtrue = zeros((self.n, 1))
+ self.assertCompatibleSystem(A, xtrue)
+
+ def testIdentityACase2(self):
+ A = identity(self.n)
+ xtrue = ones((self.n,1))
+ self.assertCompatibleSystem(A, xtrue)
+
+ def testIdentityACase3(self):
+ A = identity(self.n)
+ xtrue = transpose(arange(self.n,0,-1))
+ self.assertCompatibleSystem(A, xtrue)
+
+ def testBidiagonalA(self):
+ A = lowerBidiagonalMatrix(20,self.n)
+ xtrue = transpose(arange(self.n,0,-1))
+ self.assertCompatibleSystem(A,xtrue)
+
+class TestLSMRReturns:
+ def setUp(self):
+ self.n = 10
+ self.A = lowerBidiagonalMatrix(20,self.n)
+ self.xtrue = transpose(arange(self.n,0,-1))
+ self.Afun = aslinearoperator(self.A)
+ self.b = self.Afun.matvec(self.xtrue)
+ self.returnValues = lsmr(self.A,self.b)
+
+ def testNormr(self):
+ x, istop, itn, normr, normar, normA, condA, normx = self.returnValues;
+ assert_almost_equal(normr, norm(self.b - self.Afun.matvec(x)))
+
+ def testNormar(self):
+ x, istop, itn, normr, normar, normA, condA, normx = self.returnValues;
+ assert_almost_equal(normar, \
+ norm(self.Afun.rmatvec(self.b - self.Afun.matvec(x))))
+
+ def testNormx(self):
+ x, istop, itn, normr, normar, normA, condA, normx = self.returnValues;
+ assert_almost_equal(normx, norm(x))
+
+def lowerBidiagonalMatrix(m, n):
+ # This is a simple example for testing LSMR.
+ # It uses the leading m*n submatrix from
+ # A = [ 1
+ # 1 2
+ # 2 3
+ # 3 4
+ # ...
+ # n ]
+ # suitably padded by zeros.
+ #
+ # 04 Jun 2010: First version for distribution with lsmr.py
+ if m <= n:
+ row = hstack((arange(m, dtype=int), \
+ arange(1, m, dtype=int)))
+ col = hstack((arange(m, dtype=int), \
+ arange(m-1, dtype=int)))
+ data = hstack((arange(1, m+1, dtype=float), \
+ arange(1,m, dtype=float)))
+ return coo_matrix((data, (row, col)), shape=(m,n))
+ else:
+ row = hstack((arange(n, dtype=int), \
+ arange(1, n+1, dtype=int)))
+ col = hstack((arange(n, dtype=int), \
+ arange(n, dtype=int)))
+ data = hstack((arange(1, n+1, dtype=float), \
+ arange(1,n+1, dtype=float)))
+ return coo_matrix((data,(row, col)), shape=(m,n))
+
+def lsmrtest(m, n, damp):
+ """Verbose testing of lsmr"""
+
+ A = lowerBidiagonalMatrix(m,n)
+ xtrue = arange(n,0,-1, dtype=float)
+ Afun = aslinearoperator(A)
+
+ b = Afun.matvec(xtrue)
+
+ atol = 1.0e-7;
+ btol = 1.0e-7;
+ conlim = 1.0e+10;
+ itnlim = 10*n;
+ show = 1;
+
+ x, istop, itn, normr, normar, norma, conda, normx \
+ = lsmr(A, b, damp, atol, btol, conlim, itnlim, show )
+
+ j1 = min(n,5); j2 = max(n-4,1);
+ print ' '
+ print 'First elements of x:'
+ str = [ '%10.4f' %(xi) for xi in x[0:j1] ]
+ print ''.join(str)
+ print ' '
+ print 'Last elements of x:'
+ str = [ '%10.4f' %(xi) for xi in x[j2-1:] ]
+ print ''.join(str)
+
+ r = b - Afun.matvec(x);
+ r2 = sqrt(norm(r)**2 + (damp*norm(x))**2)
+ print ' '
+ str = 'normr (est.) %17.10e' %(normr )
+ str2 = 'normr (true) %17.10e' %(r2 )
+ print str
+ print str2
+ print ' '
+
+if __name__ == "__main__":
+ # Comment out the next line to run unit tests only
+ lsmrtest(20,10,0)
+ run_module_suite()
Something went wrong with that request. Please try again.