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

1"""Sparse Equations and Least Squares. 

2 

3The original Fortran code was written by C. C. Paige and M. A. Saunders as 

4described in 

5 

6C. C. Paige and M. A. Saunders, LSQR: An algorithm for sparse linear 

7equations and sparse least squares, TOMS 8(1), 43--71 (1982). 

8 

9C. C. Paige and M. A. Saunders, Algorithm 583; LSQR: Sparse linear 

10equations and least-squares problems, TOMS 8(2), 195--209 (1982). 

11 

12It is licensed under the following BSD license: 

13 

14Copyright (c) 2006, Systems Optimization Laboratory 

15All rights reserved. 

16 

17Redistribution and use in source and binary forms, with or without 

18modification, are permitted provided that the following conditions are 

19met: 

20 

21 * Redistributions of source code must retain the above copyright 

22 notice, this list of conditions and the following disclaimer. 

23 

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. 

28 

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. 

32 

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. 

44 

45The Fortran code was translated to Python for use in CVXOPT by Jeffery 

46Kline with contributions by Mridul Aanjaneya and Bob Myhill. 

47 

48Adapted for SciPy by Stefan van der Walt. 

49 

50""" 

51 

52__all__ = ['lsqr'] 

53 

54import numpy as np 

55from math import sqrt 

56from scipy.sparse.linalg._interface import aslinearoperator 

57 

58eps = np.finfo(np.float64).eps 

59 

60 

61def _sym_ortho(a, b): 

62 """ 

63 Stable implementation of Givens rotation. 

64 

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). 

71 

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 

77 

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 

94 

95 

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. 

100 

101 The function solves ``Ax = b`` or ``min ||Ax - b||^2`` or 

102 ``min ||Ax - b||^2 + d^2 ||x - x0||^2``. 

103 

104 The matrix A may be square or rectangular (over-determined or 

105 under-determined), and may have any rank. 

106 

107 :: 

108 

109 1. Unsymmetric equations -- solve Ax = b 

110 

111 2. Linear least squares -- solve Ax = b 

112 in the least-squares sense 

113 

114 3. Damped least squares -- solve ( A )*x = ( b ) 

115 ( damp*I ) ( damp*x0 ) 

116 in the least-squares sense 

117 

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. 

163 

164 .. versionadded:: 1.0.0 

165 

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.``) 

195 

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. 

202 

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. 

207 

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). 

212 

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. 

216 

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. 

222 

223 If some initial estimate ``x0`` is known and if ``damp == 0``, 

224 one could proceed as follows: 

225 

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``. 

229 

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. 

240 

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. 

247 

248 If A is symmetric, LSQR should not be used! 

249 

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). 

256 

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. 

267 

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) 

274 

275 The first example has the trivial solution ``[0, 0]`` 

276 

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.]) 

283 

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: 

287 

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 

298 

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: 

306 

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 

317 

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() 

326 

327 m, n = A.shape 

328 if iter_lim is None: 

329 iter_lim = 2 * n 

330 var = np.zeros(n) 

331 

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 ') 

340 

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) 

352 

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 

368 

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) 

373 

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) 

381 

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 

389 

390 if alfa > 0: 

391 v = (1/alfa) * v 

392 w = v.copy() 

393 

394 rhobar = alfa 

395 phibar = beta 

396 rnorm = beta 

397 r1norm = rnorm 

398 r2norm = rnorm 

399 

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 

407 

408 head1 = ' Itn x[0] r1norm r2norm ' 

409 head2 = ' Compatible LS Norm A Cond A' 

410 

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) 

420 

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) 

430 

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 

438 

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. 

451 

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) 

455 

456 theta = sn * alfa 

457 rhobar = -cs * alfa 

458 phi = cs * phibar 

459 phibar = sn * phibar 

460 tau = sn * phi 

461 

462 # Update x and w. 

463 t1 = phi / rho 

464 t2 = -theta / rho 

465 dk = (1 / rho) * w 

466 

467 x = x + t1 * w 

468 w = v + t2 * w 

469 ddnorm = ddnorm + np.linalg.norm(dk)**2 

470 

471 if calc_var: 

472 var = var + dk**2 

473 

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 

487 

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) 

496 

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 

512 

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 

520 

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 

534 

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 

542 

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 

561 

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) 

568 

569 if istop != 0: 

570 break 

571 

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(' ') 

586 

587 return x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var