Skip to content
Browse files

STY: Clean up LSMR patch.

  • Loading branch information...
1 parent 8971e94 commit ee7caa113faebfa1b90a75b3a5591082adb57dac @stefanv stefanv committed
Showing with 150 additions and 177 deletions.
  1. +1 −1 THANKS.txt
  2. +136 −151 scipy/sparse/linalg/isolve/lsmr.py
  3. +13 −25 scipy/sparse/linalg/isolve/lsqr.py
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
287 scipy/sparse/linalg/isolve/lsmr.py
@@ -24,6 +24,8 @@
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):
"""Iterative solver for least-squares problems.
@@ -37,83 +39,91 @@ def lsmr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
Parameters
----------
A : {matrix, sparse matrix, ndarray, LinearOperator}
- matrix A in the linear system
+ Matrix A in the linear system.
b : (m,) ndarray
- vector b in the linear system
+ Vector b in the linear system.
damp : float
- Damping factor for regularized least-squares. lsmr solves
- the regularized least-squares problem
+ 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 - A*x be the residual vector for
- the current approximate solution x. If A*x = 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.
+ Stopping tolerances. `lsmr` continues iterations until a
+ certain backward error estimate is smaller than some quantity
+ depending on atol and btol. Let ``r = b - A*x`` be the
+ residual vector for the current approximate solution ``x``.
+ If ``A*x = 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 CONLIM = 1e+8.
- Maximum precision can be obtained by setting atol = btol =
- conlim = 0, but the number of iterations may then be excessive.
+ `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.
itnlim : int
- lsmr terminates if the number of iterations reaches itnlim.
- The default is itnlim = min(m,n). For ill-conditioned systems,
- a larger value of itnlim may be needed.
+ `lsmr` terminates if the number of iterations reaches
+ `itnlim`. The default is ``itnlim = min(m,n)``. For
+ ill-conditioned systems, a larger value of `itnlim` may be
+ needed.
show : bool
- print iterations logs if show=True
+ Print iterations logs if ``show=True``.
Returns
-------
x : ndarray of float
- least-square solution returned
+ 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 itnlim before the other stopping
- conditions were satisfied.
+ 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 itnlim before the other stopping
+ conditions were satisfied.
+
itn : int
- number of iterations used
+ Number of iterations used.
normr : float
- norm(b-A*x)
+ ``norm(b-A*x)``
normar : float
- norm(A^T *(b-A*x))
+ ``norm(A^T *(b-A*x))``
norma : float
- norm(A)
+ ``norm(A)``
conda : float
- condition number of A
+ Condition number of A.
normx : float
- norm(x)
+ ``norm(x)``
- Reference
- ---------
- http://arxiv.org/abs/1006.0758
- D. C.-L. Fong and M. A. Saunders
- LSMR: An iterative algorithm for least-square problems
- Draft of 01 Jun 2010, submitted to SISC.
+ 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
"""
@@ -121,13 +131,13 @@ def lsmr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
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 ')
+ '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'
@@ -160,18 +170,18 @@ def lsmr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
alpha = 0
if beta > 0:
- u = (1/beta) * u
+ u = (1 / beta) * u
v = A.rmatvec(u)
alpha = norm(v)
if alpha > 0:
- v = (1/alpha) * v
+ v = (1 / alpha) * v
# Initialize variables for 1st iteration.
itn = 0
- zetabar = alpha*beta
+ zetabar = alpha * beta
alphabar = alpha
rho = 1
rhobar = 1
@@ -194,7 +204,7 @@ def lsmr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
# Initialize variables for estimation of ||A|| and cond(A)
- normA2 = alpha*alpha
+ normA2 = alpha * alpha
maxrbar = 0
minrbar = 1e+100
normA = sqrt(normA2)
@@ -205,7 +215,7 @@ def lsmr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
normb = beta
istop = 0
ctol = 0
- if conlim > 0: ctol = 1/conlim
+ if conlim > 0: ctol = 1 / conlim
normr = beta
# Reverse the order here from the original matlab code because
@@ -238,7 +248,7 @@ def lsmr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
beta = norm(u)
if beta > 0:
- u = (1/beta) * u
+ u = (1 / beta) * u
v = A.rmatvec(u) - beta * v
alpha = norm(v)
if alpha > 0:
@@ -248,68 +258,68 @@ def lsmr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
# Construct rotation Qhat_{k,2k+1}.
- chat, shat, alphahat = symOrtho(alphabar, damp)
+ 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 = symOrtho(alphahat, beta)
+ 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 = symOrtho(cbar*rho, thetanew)
- zeta = cbar*zetabar
- zetabar = - sbar*zetabar
+ 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
+ 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
+ betaacute = chat * betadd
+ betacheck = -shat * betadd
# Apply rotation Q_{k,k+1}.
- betahat = c*betaacute
- betadd = - s*betaacute
+ betahat = c * betaacute
+ betadd = -s * betaacute
# Apply rotation Qtilde_{k-1}.
# betad = betad_{k-1} here.
thetatildeold = thetatilde
- ctildeold, stildeold, rhotildeold = symOrtho(rhodold, thetabar)
- thetatilde = stildeold* rhobar
- rhodold = ctildeold* rhobar
- betad = - stildeold*betad + ctildeold*betahat
+ 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)
+ 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
+ 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)
+ maxrbar = max(maxrbar, rhobarold)
+ if itn > 1:
+ minrbar= min(minrbar, rhobarold)
+ condA = max(maxrbar, rhotemp) / min(minrbar, rhotemp)
# Test for convergence.
@@ -320,11 +330,11 @@ def lsmr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
# Now use these norms to estimate certain other quantities,
# some of which will be small near a solution.
- test1 = normr /normb
- test2 = normar/(normA*normr)
- test3 = 1/condA
- t1 = test1/(1 + normA*normx/normb)
- rtol = btol + atol*normA*normx/normb
+ test1 = normr / normb
+ test2 = normar / (normA * normr)
+ 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
@@ -332,29 +342,29 @@ def lsmr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
# The effect is equivalent to the normAl tests using
# atol = eps, btol = eps, conlim = 1/eps.
- if itn >= itnlim: istop = 7
- if 1 + test3 <= 1: istop = 6
- if 1 + test2 <= 1: istop = 5
- if 1 + t1 <= 1: istop = 4
+ if itn >= itnlim: 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
+ 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:
prnt = 0
- if n <= 40 : prnt = 1
- if itn <= 10 : prnt = 1
- if itn >= itnlim-10: prnt = 1
- if itn % 10 == 0 : prnt = 1
- if test3 <= 1.1*ctol : prnt = 1
- if test2 <= 1.1*atol : prnt = 1
- if test1 <= 1.1*rtol : prnt = 1
- if istop != 0 : prnt = 1
+ if n <= 40 : prnt = 1
+ if itn <= 10 : prnt = 1
+ if itn >= itnlim-10: prnt = 1
+ if itn % 10 == 0 : prnt = 1
+ if test3 <= 1.1 * ctol : prnt = 1
+ if test2 <= 1.1 * atol : prnt = 1
+ if test1 <= 1.1 * rtol : prnt = 1
+ if istop != 0 : prnt = 1
if prnt:
if pcount >= pfreq:
@@ -362,10 +372,10 @@ def lsmr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
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 )
+ 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
@@ -376,36 +386,11 @@ def lsmr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
print ' '
print 'LSMR finished'
print msg[istop]
- str1 = 'istop =%8g normr =%8.1e' %( istop, normr )
- str2 = ' normA =%8.1e normAr =%8.1e' %( normA, normar)
- str3 = 'itn =%8g condA =%8.1e' %( itn , condA )
- str4 = ' normx =%8.1e' %( normx)
+ str1 = 'istop =%8g normr =%8.1e' % (istop, normr)
+ str2 = ' normA =%8.1e normAr =%8.1e' % (normA, normar)
+ str3 = 'itn =%8g condA =%8.1e' % (itn, condA)
+ str4 = ' normx =%8.1e' % (normx)
print str1, str2
print str3, str4
return x, istop, itn, normr, normar, normA, condA, normx
-
-def sign(a):
- if a < 0: return -1
- return 1
-
-def symOrtho(a,b):
- """
- A stable implementation of Givens rotation according to
- S.-C. Choi, "Iterative Methods for Singular Linear Equations
- and Least-Squares Problems", Dissertation,
- http://www.stanford.edu/group/SOL/dissertations/sou-cheng-choi-thesis.pdf
- """
- if b==0: return sign(a), 0, abs(a)
- elif a==0: return 0, sign(b), abs(b)
- elif abs(b)>abs(a):
- tau = a / b
- s = sign(b) / sqrt(1+tau*tau)
- c = s * tau
- r = b / s
- else:
- tau = b / a
- c = sign(a) / sqrt(1+tau*tau)
- s = c * tau
- r = a / c
- return c, s, r
View
38 scipy/sparse/linalg/isolve/lsqr.py
@@ -55,9 +55,13 @@
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
+ Stable implementation of Givens rotation.
+
+ Notes
+ -----
+ Jeffery Kline: 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
@@ -70,34 +74,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

0 comments on commit ee7caa1

Please sign in to comment.
Something went wrong with that request. Please try again.