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