Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/sparse/linalg/_isolve/lsmr.py: 4%
187 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""
2Copyright (C) 2010 David Fong and Michael Saunders
4LSMR uses an iterative method.
607 Jun 2010: Documentation updated
703 Jun 2010: First release version in Python
9David Chin-lung Fong clfong@stanford.edu
10Institute for Computational and Mathematical Engineering
11Stanford University
13Michael Saunders saunders@stanford.edu
14Systems Optimization Laboratory
15Dept of MS&E, Stanford University.
17"""
19__all__ = ['lsmr']
21from numpy import zeros, infty, atleast_1d, result_type
22from numpy.linalg import norm
23from math import sqrt
24from scipy.sparse.linalg._interface import aslinearoperator
26from scipy.sparse.linalg._isolve.lsqr import _sym_ortho
29def lsmr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
30 maxiter=None, show=False, x0=None):
31 """Iterative solver for least-squares problems.
33 lsmr solves the system of linear equations ``Ax = b``. If the system
34 is inconsistent, it solves the least-squares problem ``min ||b - Ax||_2``.
35 ``A`` is a rectangular matrix of dimension m-by-n, where all cases are
36 allowed: m = n, m > n, or m < n. ``b`` is a vector of length m.
37 The matrix A may be dense or sparse (usually sparse).
39 Parameters
40 ----------
41 A : {sparse matrix, ndarray, LinearOperator}
42 Matrix A in the linear system.
43 Alternatively, ``A`` can be a linear operator which can
44 produce ``Ax`` and ``A^H x`` using, e.g.,
45 ``scipy.sparse.linalg.LinearOperator``.
46 b : array_like, shape (m,)
47 Vector ``b`` in the linear system.
48 damp : float
49 Damping factor for regularized least-squares. `lsmr` solves
50 the regularized least-squares problem::
52 min ||(b) - ( A )x||
53 ||(0) (damp*I) ||_2
55 where damp is a scalar. If damp is None or 0, the system
56 is solved without regularization. Default is 0.
57 atol, btol : float, optional
58 Stopping tolerances. `lsmr` continues iterations until a
59 certain backward error estimate is smaller than some quantity
60 depending on atol and btol. Let ``r = b - Ax`` be the
61 residual vector for the current approximate solution ``x``.
62 If ``Ax = b`` seems to be consistent, `lsmr` terminates
63 when ``norm(r) <= atol * norm(A) * norm(x) + btol * norm(b)``.
64 Otherwise, `lsmr` terminates when ``norm(A^H r) <=
65 atol * norm(A) * norm(r)``. If both tolerances are 1.0e-6 (default),
66 the final ``norm(r)`` should be accurate to about 6
67 digits. (The final ``x`` will usually have fewer correct digits,
68 depending on ``cond(A)`` and the size of LAMBDA.) If `atol`
69 or `btol` is None, a default value of 1.0e-6 will be used.
70 Ideally, they should be estimates of the relative error in the
71 entries of ``A`` and ``b`` respectively. For example, if the entries
72 of ``A`` have 7 correct digits, set ``atol = 1e-7``. This prevents
73 the algorithm from doing unnecessary work beyond the
74 uncertainty of the input data.
75 conlim : float, optional
76 `lsmr` terminates if an estimate of ``cond(A)`` exceeds
77 `conlim`. For compatible systems ``Ax = b``, conlim could be
78 as large as 1.0e+12 (say). For least-squares problems,
79 `conlim` should be less than 1.0e+8. If `conlim` is None, the
80 default value is 1e+8. Maximum precision can be obtained by
81 setting ``atol = btol = conlim = 0``, but the number of
82 iterations may then be excessive. Default is 1e8.
83 maxiter : int, optional
84 `lsmr` terminates if the number of iterations reaches
85 `maxiter`. The default is ``maxiter = min(m, n)``. For
86 ill-conditioned systems, a larger value of `maxiter` may be
87 needed. Default is False.
88 show : bool, optional
89 Print iterations logs if ``show=True``. Default is False.
90 x0 : array_like, shape (n,), optional
91 Initial guess of ``x``, if None zeros are used. Default is None.
93 .. versionadded:: 1.0.0
95 Returns
96 -------
97 x : ndarray of float
98 Least-square solution returned.
99 istop : int
100 istop gives the reason for stopping::
102 istop = 0 means x=0 is a solution. If x0 was given, then x=x0 is a
103 solution.
104 = 1 means x is an approximate solution to A@x = B,
105 according to atol and btol.
106 = 2 means x approximately solves the least-squares problem
107 according to atol.
108 = 3 means COND(A) seems to be greater than CONLIM.
109 = 4 is the same as 1 with atol = btol = eps (machine
110 precision)
111 = 5 is the same as 2 with atol = eps.
112 = 6 is the same as 3 with CONLIM = 1/eps.
113 = 7 means ITN reached maxiter before the other stopping
114 conditions were satisfied.
116 itn : int
117 Number of iterations used.
118 normr : float
119 ``norm(b-Ax)``
120 normar : float
121 ``norm(A^H (b - Ax))``
122 norma : float
123 ``norm(A)``
124 conda : float
125 Condition number of A.
126 normx : float
127 ``norm(x)``
129 Notes
130 -----
132 .. versionadded:: 0.11.0
134 References
135 ----------
136 .. [1] D. C.-L. Fong and M. A. Saunders,
137 "LSMR: An iterative algorithm for sparse least-squares problems",
138 SIAM J. Sci. Comput., vol. 33, pp. 2950-2971, 2011.
139 :arxiv:`1006.0758`
140 .. [2] LSMR Software, https://web.stanford.edu/group/SOL/software/lsmr/
142 Examples
143 --------
144 >>> import numpy as np
145 >>> from scipy.sparse import csc_matrix
146 >>> from scipy.sparse.linalg import lsmr
147 >>> A = csc_matrix([[1., 0.], [1., 1.], [0., 1.]], dtype=float)
149 The first example has the trivial solution ``[0, 0]``
151 >>> b = np.array([0., 0., 0.], dtype=float)
152 >>> x, istop, itn, normr = lsmr(A, b)[:4]
153 >>> istop
154 0
155 >>> x
156 array([0., 0.])
158 The stopping code `istop=0` returned indicates that a vector of zeros was
159 found as a solution. The returned solution `x` indeed contains
160 ``[0., 0.]``. The next example has a non-trivial solution:
162 >>> b = np.array([1., 0., -1.], dtype=float)
163 >>> x, istop, itn, normr = lsmr(A, b)[:4]
164 >>> istop
165 1
166 >>> x
167 array([ 1., -1.])
168 >>> itn
169 1
170 >>> normr
171 4.440892098500627e-16
173 As indicated by `istop=1`, `lsmr` found a solution obeying the tolerance
174 limits. The given solution ``[1., -1.]`` obviously solves the equation. The
175 remaining return values include information about the number of iterations
176 (`itn=1`) and the remaining difference of left and right side of the solved
177 equation.
178 The final example demonstrates the behavior in the case where there is no
179 solution for the equation:
181 >>> b = np.array([1., 0.01, -1.], dtype=float)
182 >>> x, istop, itn, normr = lsmr(A, b)[:4]
183 >>> istop
184 2
185 >>> x
186 array([ 1.00333333, -0.99666667])
187 >>> A.dot(x)-b
188 array([ 0.00333333, -0.00333333, 0.00333333])
189 >>> normr
190 0.005773502691896255
192 `istop` indicates that the system is inconsistent and thus `x` is rather an
193 approximate solution to the corresponding least-squares problem. `normr`
194 contains the minimal distance that was found.
195 """
197 A = aslinearoperator(A)
198 b = atleast_1d(b)
199 if b.ndim > 1:
200 b = b.squeeze()
202 msg = ('The exact solution is x = 0, or x = x0, if x0 was given ',
203 'Ax - b is small enough, given atol, btol ',
204 'The least-squares solution is good enough, given atol ',
205 'The estimate of cond(Abar) has exceeded conlim ',
206 'Ax - b is small enough for this machine ',
207 'The least-squares solution is good enough for this machine',
208 'Cond(Abar) seems to be too large for this machine ',
209 'The iteration limit has been reached ')
211 hdg1 = ' itn x(1) norm r norm Ar'
212 hdg2 = ' compatible LS norm A cond A'
213 pfreq = 20 # print frequency (for repeating the heading)
214 pcount = 0 # print counter
216 m, n = A.shape
218 # stores the num of singular values
219 minDim = min([m, n])
221 if maxiter is None:
222 maxiter = minDim
224 if x0 is None:
225 dtype = result_type(A, b, float)
226 else:
227 dtype = result_type(A, b, x0, float)
229 if show:
230 print(' ')
231 print('LSMR Least-squares solution of Ax = b\n')
232 print(f'The matrix A has {m} rows and {n} columns')
233 print('damp = %20.14e\n' % (damp))
234 print('atol = %8.2e conlim = %8.2e\n' % (atol, conlim))
235 print('btol = %8.2e maxiter = %8g\n' % (btol, maxiter))
237 u = b
238 normb = norm(b)
239 if x0 is None:
240 x = zeros(n, dtype)
241 beta = normb.copy()
242 else:
243 x = atleast_1d(x0.copy())
244 u = u - A.matvec(x)
245 beta = norm(u)
247 if beta > 0:
248 u = (1 / beta) * u
249 v = A.rmatvec(u)
250 alpha = norm(v)
251 else:
252 v = zeros(n, dtype)
253 alpha = 0
255 if alpha > 0:
256 v = (1 / alpha) * v
258 # Initialize variables for 1st iteration.
260 itn = 0
261 zetabar = alpha * beta
262 alphabar = alpha
263 rho = 1
264 rhobar = 1
265 cbar = 1
266 sbar = 0
268 h = v.copy()
269 hbar = zeros(n, dtype)
271 # Initialize variables for estimation of ||r||.
273 betadd = beta
274 betad = 0
275 rhodold = 1
276 tautildeold = 0
277 thetatilde = 0
278 zeta = 0
279 d = 0
281 # Initialize variables for estimation of ||A|| and cond(A)
283 normA2 = alpha * alpha
284 maxrbar = 0
285 minrbar = 1e+100
286 normA = sqrt(normA2)
287 condA = 1
288 normx = 0
290 # Items for use in stopping rules, normb set earlier
291 istop = 0
292 ctol = 0
293 if conlim > 0:
294 ctol = 1 / conlim
295 normr = beta
297 # Reverse the order here from the original matlab code because
298 # there was an error on return when arnorm==0
299 normar = alpha * beta
300 if normar == 0:
301 if show:
302 print(msg[0])
303 return x, istop, itn, normr, normar, normA, condA, normx
305 if normb == 0:
306 x[()] = 0
307 return x, istop, itn, normr, normar, normA, condA, normx
309 if show:
310 print(' ')
311 print(hdg1, hdg2)
312 test1 = 1
313 test2 = alpha / beta
314 str1 = '%6g %12.5e' % (itn, x[0])
315 str2 = ' %10.3e %10.3e' % (normr, normar)
316 str3 = ' %8.1e %8.1e' % (test1, test2)
317 print(''.join([str1, str2, str3]))
319 # Main iteration loop.
320 while itn < maxiter:
321 itn = itn + 1
323 # Perform the next step of the bidiagonalization to obtain the
324 # next beta, u, alpha, v. These satisfy the relations
325 # beta*u = A@v - alpha*u,
326 # alpha*v = A'@u - beta*v.
328 u *= -alpha
329 u += A.matvec(v)
330 beta = norm(u)
332 if beta > 0:
333 u *= (1 / beta)
334 v *= -beta
335 v += A.rmatvec(u)
336 alpha = norm(v)
337 if alpha > 0:
338 v *= (1 / alpha)
340 # At this point, beta = beta_{k+1}, alpha = alpha_{k+1}.
342 # Construct rotation Qhat_{k,2k+1}.
344 chat, shat, alphahat = _sym_ortho(alphabar, damp)
346 # Use a plane rotation (Q_i) to turn B_i to R_i
348 rhoold = rho
349 c, s, rho = _sym_ortho(alphahat, beta)
350 thetanew = s*alpha
351 alphabar = c*alpha
353 # Use a plane rotation (Qbar_i) to turn R_i^T to R_i^bar
355 rhobarold = rhobar
356 zetaold = zeta
357 thetabar = sbar * rho
358 rhotemp = cbar * rho
359 cbar, sbar, rhobar = _sym_ortho(cbar * rho, thetanew)
360 zeta = cbar * zetabar
361 zetabar = - sbar * zetabar
363 # Update h, h_hat, x.
365 hbar *= - (thetabar * rho / (rhoold * rhobarold))
366 hbar += h
367 x += (zeta / (rho * rhobar)) * hbar
368 h *= - (thetanew / rho)
369 h += v
371 # Estimate of ||r||.
373 # Apply rotation Qhat_{k,2k+1}.
374 betaacute = chat * betadd
375 betacheck = -shat * betadd
377 # Apply rotation Q_{k,k+1}.
378 betahat = c * betaacute
379 betadd = -s * betaacute
381 # Apply rotation Qtilde_{k-1}.
382 # betad = betad_{k-1} here.
384 thetatildeold = thetatilde
385 ctildeold, stildeold, rhotildeold = _sym_ortho(rhodold, thetabar)
386 thetatilde = stildeold * rhobar
387 rhodold = ctildeold * rhobar
388 betad = - stildeold * betad + ctildeold * betahat
390 # betad = betad_k here.
391 # rhodold = rhod_k here.
393 tautildeold = (zetaold - thetatildeold * tautildeold) / rhotildeold
394 taud = (zeta - thetatilde * tautildeold) / rhodold
395 d = d + betacheck * betacheck
396 normr = sqrt(d + (betad - taud)**2 + betadd * betadd)
398 # Estimate ||A||.
399 normA2 = normA2 + beta * beta
400 normA = sqrt(normA2)
401 normA2 = normA2 + alpha * alpha
403 # Estimate cond(A).
404 maxrbar = max(maxrbar, rhobarold)
405 if itn > 1:
406 minrbar = min(minrbar, rhobarold)
407 condA = max(maxrbar, rhotemp) / min(minrbar, rhotemp)
409 # Test for convergence.
411 # Compute norms for convergence testing.
412 normar = abs(zetabar)
413 normx = norm(x)
415 # Now use these norms to estimate certain other quantities,
416 # some of which will be small near a solution.
418 test1 = normr / normb
419 if (normA * normr) != 0:
420 test2 = normar / (normA * normr)
421 else:
422 test2 = infty
423 test3 = 1 / condA
424 t1 = test1 / (1 + normA * normx / normb)
425 rtol = btol + atol * normA * normx / normb
427 # The following tests guard against extremely small values of
428 # atol, btol or ctol. (The user may have set any or all of
429 # the parameters atol, btol, conlim to 0.)
430 # The effect is equivalent to the normAl tests using
431 # atol = eps, btol = eps, conlim = 1/eps.
433 if itn >= maxiter:
434 istop = 7
435 if 1 + test3 <= 1:
436 istop = 6
437 if 1 + test2 <= 1:
438 istop = 5
439 if 1 + t1 <= 1:
440 istop = 4
442 # Allow for tolerances set by the user.
444 if test3 <= ctol:
445 istop = 3
446 if test2 <= atol:
447 istop = 2
448 if test1 <= rtol:
449 istop = 1
451 # See if it is time to print something.
453 if show:
454 if (n <= 40) or (itn <= 10) or (itn >= maxiter - 10) or \
455 (itn % 10 == 0) or (test3 <= 1.1 * ctol) or \
456 (test2 <= 1.1 * atol) or (test1 <= 1.1 * rtol) or \
457 (istop != 0):
459 if pcount >= pfreq:
460 pcount = 0
461 print(' ')
462 print(hdg1, hdg2)
463 pcount = pcount + 1
464 str1 = '%6g %12.5e' % (itn, x[0])
465 str2 = ' %10.3e %10.3e' % (normr, normar)
466 str3 = ' %8.1e %8.1e' % (test1, test2)
467 str4 = ' %8.1e %8.1e' % (normA, condA)
468 print(''.join([str1, str2, str3, str4]))
470 if istop > 0:
471 break
473 # Print the stopping condition.
475 if show:
476 print(' ')
477 print('LSMR finished')
478 print(msg[istop])
479 print('istop =%8g normr =%8.1e' % (istop, normr))
480 print(' normA =%8.1e normAr =%8.1e' % (normA, normar))
481 print('itn =%8g condA =%8.1e' % (itn, condA))
482 print(' normx =%8.1e' % (normx))
483 print(str1, str2)
484 print(str3, str4)
486 return x, istop, itn, normr, normar, normA, condA, normx