Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/sparse/linalg/_isolve/lsqr.py: 3%
205 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Sparse Equations and Least Squares.
3The original Fortran code was written by C. C. Paige and M. A. Saunders as
4described in
6C. C. Paige and M. A. Saunders, LSQR: An algorithm for sparse linear
7equations and sparse least squares, TOMS 8(1), 43--71 (1982).
9C. C. Paige and M. A. Saunders, Algorithm 583; LSQR: Sparse linear
10equations and least-squares problems, TOMS 8(2), 195--209 (1982).
12It is licensed under the following BSD license:
14Copyright (c) 2006, Systems Optimization Laboratory
15All rights reserved.
17Redistribution and use in source and binary forms, with or without
18modification, are permitted provided that the following conditions are
19met:
21 * Redistributions of source code must retain the above copyright
22 notice, this list of conditions and the following disclaimer.
24 * Redistributions in binary form must reproduce the above
25 copyright notice, this list of conditions and the following
26 disclaimer in the documentation and/or other materials provided
27 with the distribution.
29 * Neither the name of Stanford University nor the names of its
30 contributors may be used to endorse or promote products derived
31 from this software without specific prior written permission.
33THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
34"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
35LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
36A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
37OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
38SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
39LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
40DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
41THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
42(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
43OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
45The Fortran code was translated to Python for use in CVXOPT by Jeffery
46Kline with contributions by Mridul Aanjaneya and Bob Myhill.
48Adapted for SciPy by Stefan van der Walt.
50"""
52__all__ = ['lsqr']
54import numpy as np
55from math import sqrt
56from scipy.sparse.linalg._interface import aslinearoperator
58eps = np.finfo(np.float64).eps
61def _sym_ortho(a, b):
62 """
63 Stable implementation of Givens rotation.
65 Notes
66 -----
67 The routine 'SymOrtho' was added for numerical stability. This is
68 recommended by S.-C. Choi in [1]_. It removes the unpleasant potential of
69 ``1/eps`` in some important places (see, for example text following
70 "Compute the next plane rotation Qk" in minres.py).
72 References
73 ----------
74 .. [1] S.-C. Choi, "Iterative Methods for Singular Linear Equations
75 and Least-Squares Problems", Dissertation,
76 http://www.stanford.edu/group/SOL/dissertations/sou-cheng-choi-thesis.pdf
78 """
79 if b == 0:
80 return np.sign(a), 0, abs(a)
81 elif a == 0:
82 return 0, np.sign(b), abs(b)
83 elif abs(b) > abs(a):
84 tau = a / b
85 s = np.sign(b) / sqrt(1 + tau * tau)
86 c = s * tau
87 r = b / s
88 else:
89 tau = b / a
90 c = np.sign(a) / sqrt(1+tau*tau)
91 s = c * tau
92 r = a / c
93 return c, s, r
96def lsqr(A, b, damp=0.0, atol=1e-6, btol=1e-6, conlim=1e8,
97 iter_lim=None, show=False, calc_var=False, x0=None):
98 """Find the least-squares solution to a large, sparse, linear system
99 of equations.
101 The function solves ``Ax = b`` or ``min ||Ax - b||^2`` or
102 ``min ||Ax - b||^2 + d^2 ||x - x0||^2``.
104 The matrix A may be square or rectangular (over-determined or
105 under-determined), and may have any rank.
107 ::
109 1. Unsymmetric equations -- solve Ax = b
111 2. Linear least squares -- solve Ax = b
112 in the least-squares sense
114 3. Damped least squares -- solve ( A )*x = ( b )
115 ( damp*I ) ( damp*x0 )
116 in the least-squares sense
118 Parameters
119 ----------
120 A : {sparse matrix, ndarray, LinearOperator}
121 Representation of an m-by-n matrix.
122 Alternatively, ``A`` can be a linear operator which can
123 produce ``Ax`` and ``A^T x`` using, e.g.,
124 ``scipy.sparse.linalg.LinearOperator``.
125 b : array_like, shape (m,)
126 Right-hand side vector ``b``.
127 damp : float
128 Damping coefficient. Default is 0.
129 atol, btol : float, optional
130 Stopping tolerances. `lsqr` continues iterations until a
131 certain backward error estimate is smaller than some quantity
132 depending on atol and btol. Let ``r = b - Ax`` be the
133 residual vector for the current approximate solution ``x``.
134 If ``Ax = b`` seems to be consistent, `lsqr` terminates
135 when ``norm(r) <= atol * norm(A) * norm(x) + btol * norm(b)``.
136 Otherwise, `lsqr` terminates when ``norm(A^H r) <=
137 atol * norm(A) * norm(r)``. If both tolerances are 1.0e-6 (default),
138 the final ``norm(r)`` should be accurate to about 6
139 digits. (The final ``x`` will usually have fewer correct digits,
140 depending on ``cond(A)`` and the size of LAMBDA.) If `atol`
141 or `btol` is None, a default value of 1.0e-6 will be used.
142 Ideally, they should be estimates of the relative error in the
143 entries of ``A`` and ``b`` respectively. For example, if the entries
144 of ``A`` have 7 correct digits, set ``atol = 1e-7``. This prevents
145 the algorithm from doing unnecessary work beyond the
146 uncertainty of the input data.
147 conlim : float, optional
148 Another stopping tolerance. lsqr terminates if an estimate of
149 ``cond(A)`` exceeds `conlim`. For compatible systems ``Ax =
150 b``, `conlim` could be as large as 1.0e+12 (say). For
151 least-squares problems, conlim should be less than 1.0e+8.
152 Maximum precision can be obtained by setting ``atol = btol =
153 conlim = zero``, but the number of iterations may then be
154 excessive. Default is 1e8.
155 iter_lim : int, optional
156 Explicit limitation on number of iterations (for safety).
157 show : bool, optional
158 Display an iteration log. Default is False.
159 calc_var : bool, optional
160 Whether to estimate diagonals of ``(A'A + damp^2*I)^{-1}``.
161 x0 : array_like, shape (n,), optional
162 Initial guess of x, if None zeros are used. Default is None.
164 .. versionadded:: 1.0.0
166 Returns
167 -------
168 x : ndarray of float
169 The final solution.
170 istop : int
171 Gives the reason for termination.
172 1 means x is an approximate solution to Ax = b.
173 2 means x approximately solves the least-squares problem.
174 itn : int
175 Iteration number upon termination.
176 r1norm : float
177 ``norm(r)``, where ``r = b - Ax``.
178 r2norm : float
179 ``sqrt( norm(r)^2 + damp^2 * norm(x - x0)^2 )``. Equal to `r1norm`
180 if ``damp == 0``.
181 anorm : float
182 Estimate of Frobenius norm of ``Abar = [[A]; [damp*I]]``.
183 acond : float
184 Estimate of ``cond(Abar)``.
185 arnorm : float
186 Estimate of ``norm(A'@r - damp^2*(x - x0))``.
187 xnorm : float
188 ``norm(x)``
189 var : ndarray of float
190 If ``calc_var`` is True, estimates all diagonals of
191 ``(A'A)^{-1}`` (if ``damp == 0``) or more generally ``(A'A +
192 damp^2*I)^{-1}``. This is well defined if A has full column
193 rank or ``damp > 0``. (Not sure what var means if ``rank(A)
194 < n`` and ``damp = 0.``)
196 Notes
197 -----
198 LSQR uses an iterative method to approximate the solution. The
199 number of iterations required to reach a certain accuracy depends
200 strongly on the scaling of the problem. Poor scaling of the rows
201 or columns of A should therefore be avoided where possible.
203 For example, in problem 1 the solution is unaltered by
204 row-scaling. If a row of A is very small or large compared to
205 the other rows of A, the corresponding row of ( A b ) should be
206 scaled up or down.
208 In problems 1 and 2, the solution x is easily recovered
209 following column-scaling. Unless better information is known,
210 the nonzero columns of A should be scaled so that they all have
211 the same Euclidean norm (e.g., 1.0).
213 In problem 3, there is no freedom to re-scale if damp is
214 nonzero. However, the value of damp should be assigned only
215 after attention has been paid to the scaling of A.
217 The parameter damp is intended to help regularize
218 ill-conditioned systems, by preventing the true solution from
219 being very large. Another aid to regularization is provided by
220 the parameter acond, which may be used to terminate iterations
221 before the computed solution becomes very large.
223 If some initial estimate ``x0`` is known and if ``damp == 0``,
224 one could proceed as follows:
226 1. Compute a residual vector ``r0 = b - A@x0``.
227 2. Use LSQR to solve the system ``A@dx = r0``.
228 3. Add the correction dx to obtain a final solution ``x = x0 + dx``.
230 This requires that ``x0`` be available before and after the call
231 to LSQR. To judge the benefits, suppose LSQR takes k1 iterations
232 to solve A@x = b and k2 iterations to solve A@dx = r0.
233 If x0 is "good", norm(r0) will be smaller than norm(b).
234 If the same stopping tolerances atol and btol are used for each
235 system, k1 and k2 will be similar, but the final solution x0 + dx
236 should be more accurate. The only way to reduce the total work
237 is to use a larger stopping tolerance for the second system.
238 If some value btol is suitable for A@x = b, the larger value
239 btol*norm(b)/norm(r0) should be suitable for A@dx = r0.
241 Preconditioning is another way to reduce the number of iterations.
242 If it is possible to solve a related system ``M@x = b``
243 efficiently, where M approximates A in some helpful way (e.g. M -
244 A has low rank or its elements are small relative to those of A),
245 LSQR may converge more rapidly on the system ``A@M(inverse)@z =
246 b``, after which x can be recovered by solving M@x = z.
248 If A is symmetric, LSQR should not be used!
250 Alternatives are the symmetric conjugate-gradient method (cg)
251 and/or SYMMLQ. SYMMLQ is an implementation of symmetric cg that
252 applies to any symmetric A and will converge more rapidly than
253 LSQR. If A is positive definite, there are other implementations
254 of symmetric cg that require slightly less work per iteration than
255 SYMMLQ (but will take the same number of iterations).
257 References
258 ----------
259 .. [1] C. C. Paige and M. A. Saunders (1982a).
260 "LSQR: An algorithm for sparse linear equations and
261 sparse least squares", ACM TOMS 8(1), 43-71.
262 .. [2] C. C. Paige and M. A. Saunders (1982b).
263 "Algorithm 583. LSQR: Sparse linear equations and least
264 squares problems", ACM TOMS 8(2), 195-209.
265 .. [3] M. A. Saunders (1995). "Solution of sparse rectangular
266 systems using LSQR and CRAIG", BIT 35, 588-604.
268 Examples
269 --------
270 >>> import numpy as np
271 >>> from scipy.sparse import csc_matrix
272 >>> from scipy.sparse.linalg import lsqr
273 >>> A = csc_matrix([[1., 0.], [1., 1.], [0., 1.]], dtype=float)
275 The first example has the trivial solution ``[0, 0]``
277 >>> b = np.array([0., 0., 0.], dtype=float)
278 >>> x, istop, itn, normr = lsqr(A, b)[:4]
279 >>> istop
280 0
281 >>> x
282 array([ 0., 0.])
284 The stopping code `istop=0` returned indicates that a vector of zeros was
285 found as a solution. The returned solution `x` indeed contains
286 ``[0., 0.]``. The next example has a non-trivial solution:
288 >>> b = np.array([1., 0., -1.], dtype=float)
289 >>> x, istop, itn, r1norm = lsqr(A, b)[:4]
290 >>> istop
291 1
292 >>> x
293 array([ 1., -1.])
294 >>> itn
295 1
296 >>> r1norm
297 4.440892098500627e-16
299 As indicated by `istop=1`, `lsqr` found a solution obeying the tolerance
300 limits. The given solution ``[1., -1.]`` obviously solves the equation. The
301 remaining return values include information about the number of iterations
302 (`itn=1`) and the remaining difference of left and right side of the solved
303 equation.
304 The final example demonstrates the behavior in the case where there is no
305 solution for the equation:
307 >>> b = np.array([1., 0.01, -1.], dtype=float)
308 >>> x, istop, itn, r1norm = lsqr(A, b)[:4]
309 >>> istop
310 2
311 >>> x
312 array([ 1.00333333, -0.99666667])
313 >>> A.dot(x)-b
314 array([ 0.00333333, -0.00333333, 0.00333333])
315 >>> r1norm
316 0.005773502691896255
318 `istop` indicates that the system is inconsistent and thus `x` is rather an
319 approximate solution to the corresponding least-squares problem. `r1norm`
320 contains the norm of the minimal residual that was found.
321 """
322 A = aslinearoperator(A)
323 b = np.atleast_1d(b)
324 if b.ndim > 1:
325 b = b.squeeze()
327 m, n = A.shape
328 if iter_lim is None:
329 iter_lim = 2 * n
330 var = np.zeros(n)
332 msg = ('The exact solution is x = 0 ',
333 'Ax - b is small enough, given atol, btol ',
334 'The least-squares solution is good enough, given atol ',
335 'The estimate of cond(Abar) has exceeded conlim ',
336 'Ax - b is small enough for this machine ',
337 'The least-squares solution is good enough for this machine',
338 'Cond(Abar) seems to be too large for this machine ',
339 'The iteration limit has been reached ')
341 if show:
342 print(' ')
343 print('LSQR Least-squares solution of Ax = b')
344 str1 = f'The matrix A has {m} rows and {n} columns'
345 str2 = 'damp = %20.14e calc_var = %8g' % (damp, calc_var)
346 str3 = 'atol = %8.2e conlim = %8.2e' % (atol, conlim)
347 str4 = 'btol = %8.2e iter_lim = %8g' % (btol, iter_lim)
348 print(str1)
349 print(str2)
350 print(str3)
351 print(str4)
353 itn = 0
354 istop = 0
355 ctol = 0
356 if conlim > 0:
357 ctol = 1/conlim
358 anorm = 0
359 acond = 0
360 dampsq = damp**2
361 ddnorm = 0
362 res2 = 0
363 xnorm = 0
364 xxnorm = 0
365 z = 0
366 cs2 = -1
367 sn2 = 0
369 # Set up the first vectors u and v for the bidiagonalization.
370 # These satisfy beta*u = b - A@x, alfa*v = A'@u.
371 u = b
372 bnorm = np.linalg.norm(b)
374 if x0 is None:
375 x = np.zeros(n)
376 beta = bnorm.copy()
377 else:
378 x = np.asarray(x0)
379 u = u - A.matvec(x)
380 beta = np.linalg.norm(u)
382 if beta > 0:
383 u = (1/beta) * u
384 v = A.rmatvec(u)
385 alfa = np.linalg.norm(v)
386 else:
387 v = x.copy()
388 alfa = 0
390 if alfa > 0:
391 v = (1/alfa) * v
392 w = v.copy()
394 rhobar = alfa
395 phibar = beta
396 rnorm = beta
397 r1norm = rnorm
398 r2norm = rnorm
400 # Reverse the order here from the original matlab code because
401 # there was an error on return when arnorm==0
402 arnorm = alfa * beta
403 if arnorm == 0:
404 if show:
405 print(msg[0])
406 return x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var
408 head1 = ' Itn x[0] r1norm r2norm '
409 head2 = ' Compatible LS Norm A Cond A'
411 if show:
412 print(' ')
413 print(head1, head2)
414 test1 = 1
415 test2 = alfa / beta
416 str1 = '%6g %12.5e' % (itn, x[0])
417 str2 = ' %10.3e %10.3e' % (r1norm, r2norm)
418 str3 = ' %8.1e %8.1e' % (test1, test2)
419 print(str1, str2, str3)
421 # Main iteration loop.
422 while itn < iter_lim:
423 itn = itn + 1
424 # Perform the next step of the bidiagonalization to obtain the
425 # next beta, u, alfa, v. These satisfy the relations
426 # beta*u = a@v - alfa*u,
427 # alfa*v = A'@u - beta*v.
428 u = A.matvec(v) - alfa * u
429 beta = np.linalg.norm(u)
431 if beta > 0:
432 u = (1/beta) * u
433 anorm = sqrt(anorm**2 + alfa**2 + beta**2 + dampsq)
434 v = A.rmatvec(u) - beta * v
435 alfa = np.linalg.norm(v)
436 if alfa > 0:
437 v = (1 / alfa) * v
439 # Use a plane rotation to eliminate the damping parameter.
440 # This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
441 if damp > 0:
442 rhobar1 = sqrt(rhobar**2 + dampsq)
443 cs1 = rhobar / rhobar1
444 sn1 = damp / rhobar1
445 psi = sn1 * phibar
446 phibar = cs1 * phibar
447 else:
448 # cs1 = 1 and sn1 = 0
449 rhobar1 = rhobar
450 psi = 0.
452 # Use a plane rotation to eliminate the subdiagonal element (beta)
453 # of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
454 cs, sn, rho = _sym_ortho(rhobar1, beta)
456 theta = sn * alfa
457 rhobar = -cs * alfa
458 phi = cs * phibar
459 phibar = sn * phibar
460 tau = sn * phi
462 # Update x and w.
463 t1 = phi / rho
464 t2 = -theta / rho
465 dk = (1 / rho) * w
467 x = x + t1 * w
468 w = v + t2 * w
469 ddnorm = ddnorm + np.linalg.norm(dk)**2
471 if calc_var:
472 var = var + dk**2
474 # Use a plane rotation on the right to eliminate the
475 # super-diagonal element (theta) of the upper-bidiagonal matrix.
476 # Then use the result to estimate norm(x).
477 delta = sn2 * rho
478 gambar = -cs2 * rho
479 rhs = phi - delta * z
480 zbar = rhs / gambar
481 xnorm = sqrt(xxnorm + zbar**2)
482 gamma = sqrt(gambar**2 + theta**2)
483 cs2 = gambar / gamma
484 sn2 = theta / gamma
485 z = rhs / gamma
486 xxnorm = xxnorm + z**2
488 # Test for convergence.
489 # First, estimate the condition of the matrix Abar,
490 # and the norms of rbar and Abar'rbar.
491 acond = anorm * sqrt(ddnorm)
492 res1 = phibar**2
493 res2 = res2 + psi**2
494 rnorm = sqrt(res1 + res2)
495 arnorm = alfa * abs(tau)
497 # Distinguish between
498 # r1norm = ||b - Ax|| and
499 # r2norm = rnorm in current code
500 # = sqrt(r1norm^2 + damp^2*||x - x0||^2).
501 # Estimate r1norm from
502 # r1norm = sqrt(r2norm^2 - damp^2*||x - x0||^2).
503 # Although there is cancellation, it might be accurate enough.
504 if damp > 0:
505 r1sq = rnorm**2 - dampsq * xxnorm
506 r1norm = sqrt(abs(r1sq))
507 if r1sq < 0:
508 r1norm = -r1norm
509 else:
510 r1norm = rnorm
511 r2norm = rnorm
513 # Now use these norms to estimate certain other quantities,
514 # some of which will be small near a solution.
515 test1 = rnorm / bnorm
516 test2 = arnorm / (anorm * rnorm + eps)
517 test3 = 1 / (acond + eps)
518 t1 = test1 / (1 + anorm * xnorm / bnorm)
519 rtol = btol + atol * anorm * xnorm / bnorm
521 # The following tests guard against extremely small values of
522 # atol, btol or ctol. (The user may have set any or all of
523 # the parameters atol, btol, conlim to 0.)
524 # The effect is equivalent to the normal tests using
525 # atol = eps, btol = eps, conlim = 1/eps.
526 if itn >= iter_lim:
527 istop = 7
528 if 1 + test3 <= 1:
529 istop = 6
530 if 1 + test2 <= 1:
531 istop = 5
532 if 1 + t1 <= 1:
533 istop = 4
535 # Allow for tolerances set by the user.
536 if test3 <= ctol:
537 istop = 3
538 if test2 <= atol:
539 istop = 2
540 if test1 <= rtol:
541 istop = 1
543 if show:
544 # See if it is time to print something.
545 prnt = False
546 if n <= 40:
547 prnt = True
548 if itn <= 10:
549 prnt = True
550 if itn >= iter_lim-10:
551 prnt = True
552 # if itn%10 == 0: prnt = True
553 if test3 <= 2*ctol:
554 prnt = True
555 if test2 <= 10*atol:
556 prnt = True
557 if test1 <= 10*rtol:
558 prnt = True
559 if istop != 0:
560 prnt = True
562 if prnt:
563 str1 = '%6g %12.5e' % (itn, x[0])
564 str2 = ' %10.3e %10.3e' % (r1norm, r2norm)
565 str3 = ' %8.1e %8.1e' % (test1, test2)
566 str4 = ' %8.1e %8.1e' % (anorm, acond)
567 print(str1, str2, str3, str4)
569 if istop != 0:
570 break
572 # End of iteration loop.
573 # Print the stopping condition.
574 if show:
575 print(' ')
576 print('LSQR finished')
577 print(msg[istop])
578 print(' ')
579 str1 = 'istop =%8g r1norm =%8.1e' % (istop, r1norm)
580 str2 = 'anorm =%8.1e arnorm =%8.1e' % (anorm, arnorm)
581 str3 = 'itn =%8g r2norm =%8.1e' % (itn, r2norm)
582 str4 = 'acond =%8.1e xnorm =%8.1e' % (acond, xnorm)
583 print(str1 + ' ' + str2)
584 print(str3 + ' ' + str4)
585 print(' ')
587 return x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var