Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/sparse/linalg/_isolve/minres.py: 3%
208 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1from numpy import inner, zeros, inf, finfo
2from numpy.linalg import norm
3from math import sqrt
5from .utils import make_system
7__all__ = ['minres']
10def minres(A, b, x0=None, shift=0.0, tol=1e-5, maxiter=None,
11 M=None, callback=None, show=False, check=False):
12 """
13 Use MINimum RESidual iteration to solve Ax=b
15 MINRES minimizes norm(Ax - b) for a real symmetric matrix A. Unlike
16 the Conjugate Gradient method, A can be indefinite or singular.
18 If shift != 0 then the method solves (A - shift*I)x = b
20 Parameters
21 ----------
22 A : {sparse matrix, ndarray, LinearOperator}
23 The real symmetric N-by-N matrix of the linear system
24 Alternatively, ``A`` can be a linear operator which can
25 produce ``Ax`` using, e.g.,
26 ``scipy.sparse.linalg.LinearOperator``.
27 b : ndarray
28 Right hand side of the linear system. Has shape (N,) or (N,1).
30 Returns
31 -------
32 x : ndarray
33 The converged solution.
34 info : integer
35 Provides convergence information:
36 0 : successful exit
37 >0 : convergence to tolerance not achieved, number of iterations
38 <0 : illegal input or breakdown
40 Other Parameters
41 ----------------
42 x0 : ndarray
43 Starting guess for the solution.
44 shift : float
45 Value to apply to the system ``(A - shift * I)x = b``. Default is 0.
46 tol : float
47 Tolerance to achieve. The algorithm terminates when the relative
48 residual is below `tol`.
49 maxiter : integer
50 Maximum number of iterations. Iteration will stop after maxiter
51 steps even if the specified tolerance has not been achieved.
52 M : {sparse matrix, ndarray, LinearOperator}
53 Preconditioner for A. The preconditioner should approximate the
54 inverse of A. Effective preconditioning dramatically improves the
55 rate of convergence, which implies that fewer iterations are needed
56 to reach a given error tolerance.
57 callback : function
58 User-supplied function to call after each iteration. It is called
59 as callback(xk), where xk is the current solution vector.
60 show : bool
61 If ``True``, print out a summary and metrics related to the solution
62 during iterations. Default is ``False``.
63 check : bool
64 If ``True``, run additional input validation to check that `A` and
65 `M` (if specified) are symmetric. Default is ``False``.
67 Examples
68 --------
69 >>> import numpy as np
70 >>> from scipy.sparse import csc_matrix
71 >>> from scipy.sparse.linalg import minres
72 >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)
73 >>> A = A + A.T
74 >>> b = np.array([2, 4, -1], dtype=float)
75 >>> x, exitCode = minres(A, b)
76 >>> print(exitCode) # 0 indicates successful convergence
77 0
78 >>> np.allclose(A.dot(x), b)
79 True
81 References
82 ----------
83 Solution of sparse indefinite systems of linear equations,
84 C. C. Paige and M. A. Saunders (1975),
85 SIAM J. Numer. Anal. 12(4), pp. 617-629.
86 https://web.stanford.edu/group/SOL/software/minres/
88 This file is a translation of the following MATLAB implementation:
89 https://web.stanford.edu/group/SOL/software/minres/minres-matlab.zip
91 """
92 A, M, x, b, postprocess = make_system(A, M, x0, b)
94 matvec = A.matvec
95 psolve = M.matvec
97 first = 'Enter minres. '
98 last = 'Exit minres. '
100 n = A.shape[0]
102 if maxiter is None:
103 maxiter = 5 * n
105 msg = [' beta2 = 0. If M = I, b and x are eigenvectors ', # -1
106 ' beta1 = 0. The exact solution is x0 ', # 0
107 ' A solution to Ax = b was found, given rtol ', # 1
108 ' A least-squares solution was found, given rtol ', # 2
109 ' Reasonable accuracy achieved, given eps ', # 3
110 ' x has converged to an eigenvector ', # 4
111 ' acond has exceeded 0.1/eps ', # 5
112 ' The iteration limit was reached ', # 6
113 ' A does not define a symmetric matrix ', # 7
114 ' M does not define a symmetric matrix ', # 8
115 ' M does not define a pos-def preconditioner '] # 9
117 if show:
118 print(first + 'Solution of symmetric Ax = b')
119 print(first + 'n = %3g shift = %23.14e' % (n,shift))
120 print(first + 'itnlim = %3g rtol = %11.2e' % (maxiter,tol))
121 print()
123 istop = 0
124 itn = 0
125 Anorm = 0
126 Acond = 0
127 rnorm = 0
128 ynorm = 0
130 xtype = x.dtype
132 eps = finfo(xtype).eps
134 # Set up y and v for the first Lanczos vector v1.
135 # y = beta1 P' v1, where P = C**(-1).
136 # v is really P' v1.
138 if x0 is None:
139 r1 = b.copy()
140 else:
141 r1 = b - A@x
142 y = psolve(r1)
144 beta1 = inner(r1, y)
146 if beta1 < 0:
147 raise ValueError('indefinite preconditioner')
148 elif beta1 == 0:
149 return (postprocess(x), 0)
151 bnorm = norm(b)
152 if bnorm == 0:
153 x = b
154 return (postprocess(x), 0)
156 beta1 = sqrt(beta1)
158 if check:
159 # are these too strict?
161 # see if A is symmetric
162 w = matvec(y)
163 r2 = matvec(w)
164 s = inner(w,w)
165 t = inner(y,r2)
166 z = abs(s - t)
167 epsa = (s + eps) * eps**(1.0/3.0)
168 if z > epsa:
169 raise ValueError('non-symmetric matrix')
171 # see if M is symmetric
172 r2 = psolve(y)
173 s = inner(y,y)
174 t = inner(r1,r2)
175 z = abs(s - t)
176 epsa = (s + eps) * eps**(1.0/3.0)
177 if z > epsa:
178 raise ValueError('non-symmetric preconditioner')
180 # Initialize other quantities
181 oldb = 0
182 beta = beta1
183 dbar = 0
184 epsln = 0
185 qrnorm = beta1
186 phibar = beta1
187 rhs1 = beta1
188 rhs2 = 0
189 tnorm2 = 0
190 gmax = 0
191 gmin = finfo(xtype).max
192 cs = -1
193 sn = 0
194 w = zeros(n, dtype=xtype)
195 w2 = zeros(n, dtype=xtype)
196 r2 = r1
198 if show:
199 print()
200 print()
201 print(' Itn x(1) Compatible LS norm(A) cond(A) gbar/|A|')
203 while itn < maxiter:
204 itn += 1
206 s = 1.0/beta
207 v = s*y
209 y = matvec(v)
210 y = y - shift * v
212 if itn >= 2:
213 y = y - (beta/oldb)*r1
215 alfa = inner(v,y)
216 y = y - (alfa/beta)*r2
217 r1 = r2
218 r2 = y
219 y = psolve(r2)
220 oldb = beta
221 beta = inner(r2,y)
222 if beta < 0:
223 raise ValueError('non-symmetric matrix')
224 beta = sqrt(beta)
225 tnorm2 += alfa**2 + oldb**2 + beta**2
227 if itn == 1:
228 if beta/beta1 <= 10*eps:
229 istop = -1 # Terminate later
231 # Apply previous rotation Qk-1 to get
232 # [deltak epslnk+1] = [cs sn][dbark 0 ]
233 # [gbar k dbar k+1] [sn -cs][alfak betak+1].
235 oldeps = epsln
236 delta = cs * dbar + sn * alfa # delta1 = 0 deltak
237 gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k
238 epsln = sn * beta # epsln2 = 0 epslnk+1
239 dbar = - cs * beta # dbar 2 = beta2 dbar k+1
240 root = norm([gbar, dbar])
241 Arnorm = phibar * root
243 # Compute the next plane rotation Qk
245 gamma = norm([gbar, beta]) # gammak
246 gamma = max(gamma, eps)
247 cs = gbar / gamma # ck
248 sn = beta / gamma # sk
249 phi = cs * phibar # phik
250 phibar = sn * phibar # phibark+1
252 # Update x.
254 denom = 1.0/gamma
255 w1 = w2
256 w2 = w
257 w = (v - oldeps*w1 - delta*w2) * denom
258 x = x + phi*w
260 # Go round again.
262 gmax = max(gmax, gamma)
263 gmin = min(gmin, gamma)
264 z = rhs1 / gamma
265 rhs1 = rhs2 - delta*z
266 rhs2 = - epsln*z
268 # Estimate various norms and test for convergence.
270 Anorm = sqrt(tnorm2)
271 ynorm = norm(x)
272 epsa = Anorm * eps
273 epsx = Anorm * ynorm * eps
274 epsr = Anorm * ynorm * tol
275 diag = gbar
277 if diag == 0:
278 diag = epsa
280 qrnorm = phibar
281 rnorm = qrnorm
282 if ynorm == 0 or Anorm == 0:
283 test1 = inf
284 else:
285 test1 = rnorm / (Anorm*ynorm) # ||r|| / (||A|| ||x||)
286 if Anorm == 0:
287 test2 = inf
288 else:
289 test2 = root / Anorm # ||Ar|| / (||A|| ||r||)
291 # Estimate cond(A).
292 # In this version we look at the diagonals of R in the
293 # factorization of the lower Hessenberg matrix, Q @ H = R,
294 # where H is the tridiagonal matrix from Lanczos with one
295 # extra row, beta(k+1) e_k^T.
297 Acond = gmax/gmin
299 # See if any of the stopping criteria are satisfied.
300 # In rare cases, istop is already -1 from above (Abar = const*I).
302 if istop == 0:
303 t1 = 1 + test1 # These tests work if tol < eps
304 t2 = 1 + test2
305 if t2 <= 1:
306 istop = 2
307 if t1 <= 1:
308 istop = 1
310 if itn >= maxiter:
311 istop = 6
312 if Acond >= 0.1/eps:
313 istop = 4
314 if epsx >= beta1:
315 istop = 3
316 # if rnorm <= epsx : istop = 2
317 # if rnorm <= epsr : istop = 1
318 if test2 <= tol:
319 istop = 2
320 if test1 <= tol:
321 istop = 1
323 # See if it is time to print something.
325 prnt = False
326 if n <= 40:
327 prnt = True
328 if itn <= 10:
329 prnt = True
330 if itn >= maxiter-10:
331 prnt = True
332 if itn % 10 == 0:
333 prnt = True
334 if qrnorm <= 10*epsx:
335 prnt = True
336 if qrnorm <= 10*epsr:
337 prnt = True
338 if Acond <= 1e-2/eps:
339 prnt = True
340 if istop != 0:
341 prnt = True
343 if show and prnt:
344 str1 = '%6g %12.5e %10.3e' % (itn, x[0], test1)
345 str2 = ' %10.3e' % (test2,)
346 str3 = ' %8.1e %8.1e %8.1e' % (Anorm, Acond, gbar/Anorm)
348 print(str1 + str2 + str3)
350 if itn % 10 == 0:
351 print()
353 if callback is not None:
354 callback(x)
356 if istop != 0:
357 break # TODO check this
359 if show:
360 print()
361 print(last + ' istop = %3g itn =%5g' % (istop,itn))
362 print(last + ' Anorm = %12.4e Acond = %12.4e' % (Anorm,Acond))
363 print(last + ' rnorm = %12.4e ynorm = %12.4e' % (rnorm,ynorm))
364 print(last + ' Arnorm = %12.4e' % (Arnorm,))
365 print(last + msg[istop+1])
367 if istop == 6:
368 info = maxiter
369 else:
370 info = 0
372 return (postprocess(x),info)
375if __name__ == '__main__':
376 from numpy import arange
377 from scipy.sparse import spdiags
379 n = 10
381 residuals = []
383 def cb(x):
384 residuals.append(norm(b - A@x))
386 # A = poisson((10,),format='csr')
387 A = spdiags([arange(1,n+1,dtype=float)], [0], n, n, format='csr')
388 M = spdiags([1.0/arange(1,n+1,dtype=float)], [0], n, n, format='csr')
389 A.psolve = M.matvec
390 b = zeros(A.shape[0])
391 x = minres(A,b,tol=1e-12,maxiter=None,callback=cb)
392 # x = cg(A,b,x0=b,tol=1e-12,maxiter=None,callback=cb)[0]