Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1from numpy import sqrt, inner, zeros, inf, finfo 

2from numpy.linalg import norm 

3 

4from .utils import make_system 

5 

6__all__ = ['minres'] 

7 

8 

9def minres(A, b, x0=None, shift=0.0, tol=1e-5, maxiter=None, 

10 M=None, callback=None, show=False, check=False): 

11 """ 

12 Use MINimum RESidual iteration to solve Ax=b 

13 

14 MINRES minimizes norm(A*x - b) for a real symmetric matrix A. Unlike 

15 the Conjugate Gradient method, A can be indefinite or singular. 

16 

17 If shift != 0 then the method solves (A - shift*I)x = b 

18 

19 Parameters 

20 ---------- 

21 A : {sparse matrix, dense matrix, LinearOperator} 

22 The real symmetric N-by-N matrix of the linear system 

23 Alternatively, ``A`` can be a linear operator which can 

24 produce ``Ax`` using, e.g., 

25 ``scipy.sparse.linalg.LinearOperator``. 

26 b : {array, matrix} 

27 Right hand side of the linear system. Has shape (N,) or (N,1). 

28 

29 Returns 

30 ------- 

31 x : {array, matrix} 

32 The converged solution. 

33 info : integer 

34 Provides convergence information: 

35 0 : successful exit 

36 >0 : convergence to tolerance not achieved, number of iterations 

37 <0 : illegal input or breakdown 

38 

39 Other Parameters 

40 ---------------- 

41 x0 : {array, matrix} 

42 Starting guess for the solution. 

43 tol : float 

44 Tolerance to achieve. The algorithm terminates when the relative 

45 residual is below `tol`. 

46 maxiter : integer 

47 Maximum number of iterations. Iteration will stop after maxiter 

48 steps even if the specified tolerance has not been achieved. 

49 M : {sparse matrix, dense matrix, LinearOperator} 

50 Preconditioner for A. The preconditioner should approximate the 

51 inverse of A. Effective preconditioning dramatically improves the 

52 rate of convergence, which implies that fewer iterations are needed 

53 to reach a given error tolerance. 

54 callback : function 

55 User-supplied function to call after each iteration. It is called 

56 as callback(xk), where xk is the current solution vector. 

57 

58 Examples 

59 -------- 

60 >>> import numpy as np 

61 >>> from scipy.sparse import csc_matrix 

62 >>> from scipy.sparse.linalg import minres 

63 >>> A = csc_matrix([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float) 

64 >>> A = A + A.T 

65 >>> b = np.array([2, 4, -1], dtype=float) 

66 >>> x, exitCode = minres(A, b) 

67 >>> print(exitCode) # 0 indicates successful convergence 

68 0 

69 >>> np.allclose(A.dot(x), b) 

70 True 

71 

72 References 

73 ---------- 

74 Solution of sparse indefinite systems of linear equations, 

75 C. C. Paige and M. A. Saunders (1975), 

76 SIAM J. Numer. Anal. 12(4), pp. 617-629. 

77 https://web.stanford.edu/group/SOL/software/minres/ 

78 

79 This file is a translation of the following MATLAB implementation: 

80 https://web.stanford.edu/group/SOL/software/minres/minres-matlab.zip 

81 

82 """ 

83 A, M, x, b, postprocess = make_system(A, M, x0, b) 

84 

85 matvec = A.matvec 

86 psolve = M.matvec 

87 

88 first = 'Enter minres. ' 

89 last = 'Exit minres. ' 

90 

91 n = A.shape[0] 

92 

93 if maxiter is None: 

94 maxiter = 5 * n 

95 

96 msg = [' beta2 = 0. If M = I, b and x are eigenvectors ', # -1 

97 ' beta1 = 0. The exact solution is x0 ', # 0 

98 ' A solution to Ax = b was found, given rtol ', # 1 

99 ' A least-squares solution was found, given rtol ', # 2 

100 ' Reasonable accuracy achieved, given eps ', # 3 

101 ' x has converged to an eigenvector ', # 4 

102 ' acond has exceeded 0.1/eps ', # 5 

103 ' The iteration limit was reached ', # 6 

104 ' A does not define a symmetric matrix ', # 7 

105 ' M does not define a symmetric matrix ', # 8 

106 ' M does not define a pos-def preconditioner '] # 9 

107 

108 if show: 

109 print(first + 'Solution of symmetric Ax = b') 

110 print(first + 'n = %3g shift = %23.14e' % (n,shift)) 

111 print(first + 'itnlim = %3g rtol = %11.2e' % (maxiter,tol)) 

112 print() 

113 

114 istop = 0 

115 itn = 0 

116 Anorm = 0 

117 Acond = 0 

118 rnorm = 0 

119 ynorm = 0 

120 

121 xtype = x.dtype 

122 

123 eps = finfo(xtype).eps 

124 

125 # Set up y and v for the first Lanczos vector v1. 

126 # y = beta1 P' v1, where P = C**(-1). 

127 # v is really P' v1. 

128 

129 r1 = b - A*x 

130 y = psolve(r1) 

131 

132 beta1 = inner(r1, y) 

133 

134 if beta1 < 0: 

135 raise ValueError('indefinite preconditioner') 

136 elif beta1 == 0: 

137 return (postprocess(x), 0) 

138 

139 beta1 = sqrt(beta1) 

140 

141 if check: 

142 # are these too strict? 

143 

144 # see if A is symmetric 

145 w = matvec(y) 

146 r2 = matvec(w) 

147 s = inner(w,w) 

148 t = inner(y,r2) 

149 z = abs(s - t) 

150 epsa = (s + eps) * eps**(1.0/3.0) 

151 if z > epsa: 

152 raise ValueError('non-symmetric matrix') 

153 

154 # see if M is symmetric 

155 r2 = psolve(y) 

156 s = inner(y,y) 

157 t = inner(r1,r2) 

158 z = abs(s - t) 

159 epsa = (s + eps) * eps**(1.0/3.0) 

160 if z > epsa: 

161 raise ValueError('non-symmetric preconditioner') 

162 

163 # Initialize other quantities 

164 oldb = 0 

165 beta = beta1 

166 dbar = 0 

167 epsln = 0 

168 qrnorm = beta1 

169 phibar = beta1 

170 rhs1 = beta1 

171 rhs2 = 0 

172 tnorm2 = 0 

173 gmax = 0 

174 gmin = finfo(xtype).max 

175 cs = -1 

176 sn = 0 

177 w = zeros(n, dtype=xtype) 

178 w2 = zeros(n, dtype=xtype) 

179 r2 = r1 

180 

181 if show: 

182 print() 

183 print() 

184 print(' Itn x(1) Compatible LS norm(A) cond(A) gbar/|A|') 

185 

186 while itn < maxiter: 

187 itn += 1 

188 

189 s = 1.0/beta 

190 v = s*y 

191 

192 y = matvec(v) 

193 y = y - shift * v 

194 

195 if itn >= 2: 

196 y = y - (beta/oldb)*r1 

197 

198 alfa = inner(v,y) 

199 y = y - (alfa/beta)*r2 

200 r1 = r2 

201 r2 = y 

202 y = psolve(r2) 

203 oldb = beta 

204 beta = inner(r2,y) 

205 if beta < 0: 

206 raise ValueError('non-symmetric matrix') 

207 beta = sqrt(beta) 

208 tnorm2 += alfa**2 + oldb**2 + beta**2 

209 

210 if itn == 1: 

211 if beta/beta1 <= 10*eps: 

212 istop = -1 # Terminate later 

213 

214 # Apply previous rotation Qk-1 to get 

215 # [deltak epslnk+1] = [cs sn][dbark 0 ] 

216 # [gbar k dbar k+1] [sn -cs][alfak betak+1]. 

217 

218 oldeps = epsln 

219 delta = cs * dbar + sn * alfa # delta1 = 0 deltak 

220 gbar = sn * dbar - cs * alfa # gbar 1 = alfa1 gbar k 

221 epsln = sn * beta # epsln2 = 0 epslnk+1 

222 dbar = - cs * beta # dbar 2 = beta2 dbar k+1 

223 root = norm([gbar, dbar]) 

224 Arnorm = phibar * root 

225 

226 # Compute the next plane rotation Qk 

227 

228 gamma = norm([gbar, beta]) # gammak 

229 gamma = max(gamma, eps) 

230 cs = gbar / gamma # ck 

231 sn = beta / gamma # sk 

232 phi = cs * phibar # phik 

233 phibar = sn * phibar # phibark+1 

234 

235 # Update x. 

236 

237 denom = 1.0/gamma 

238 w1 = w2 

239 w2 = w 

240 w = (v - oldeps*w1 - delta*w2) * denom 

241 x = x + phi*w 

242 

243 # Go round again. 

244 

245 gmax = max(gmax, gamma) 

246 gmin = min(gmin, gamma) 

247 z = rhs1 / gamma 

248 rhs1 = rhs2 - delta*z 

249 rhs2 = - epsln*z 

250 

251 # Estimate various norms and test for convergence. 

252 

253 Anorm = sqrt(tnorm2) 

254 ynorm = norm(x) 

255 epsa = Anorm * eps 

256 epsx = Anorm * ynorm * eps 

257 epsr = Anorm * ynorm * tol 

258 diag = gbar 

259 

260 if diag == 0: 

261 diag = epsa 

262 

263 qrnorm = phibar 

264 rnorm = qrnorm 

265 if ynorm == 0 or Anorm == 0: 

266 test1 = inf 

267 else: 

268 test1 = rnorm / (Anorm*ynorm) # ||r|| / (||A|| ||x||) 

269 if Anorm == 0: 

270 test2 = inf 

271 else: 

272 test2 = root / Anorm # ||Ar|| / (||A|| ||r||) 

273 

274 # Estimate cond(A). 

275 # In this version we look at the diagonals of R in the 

276 # factorization of the lower Hessenberg matrix, Q * H = R, 

277 # where H is the tridiagonal matrix from Lanczos with one 

278 # extra row, beta(k+1) e_k^T. 

279 

280 Acond = gmax/gmin 

281 

282 # See if any of the stopping criteria are satisfied. 

283 # In rare cases, istop is already -1 from above (Abar = const*I). 

284 

285 if istop == 0: 

286 t1 = 1 + test1 # These tests work if tol < eps 

287 t2 = 1 + test2 

288 if t2 <= 1: 

289 istop = 2 

290 if t1 <= 1: 

291 istop = 1 

292 

293 if itn >= maxiter: 

294 istop = 6 

295 if Acond >= 0.1/eps: 

296 istop = 4 

297 if epsx >= beta1: 

298 istop = 3 

299 # if rnorm <= epsx : istop = 2 

300 # if rnorm <= epsr : istop = 1 

301 if test2 <= tol: 

302 istop = 2 

303 if test1 <= tol: 

304 istop = 1 

305 

306 # See if it is time to print something. 

307 

308 prnt = False 

309 if n <= 40: 

310 prnt = True 

311 if itn <= 10: 

312 prnt = True 

313 if itn >= maxiter-10: 

314 prnt = True 

315 if itn % 10 == 0: 

316 prnt = True 

317 if qrnorm <= 10*epsx: 

318 prnt = True 

319 if qrnorm <= 10*epsr: 

320 prnt = True 

321 if Acond <= 1e-2/eps: 

322 prnt = True 

323 if istop != 0: 

324 prnt = True 

325 

326 if show and prnt: 

327 str1 = '%6g %12.5e %10.3e' % (itn, x[0], test1) 

328 str2 = ' %10.3e' % (test2,) 

329 str3 = ' %8.1e %8.1e %8.1e' % (Anorm, Acond, gbar/Anorm) 

330 

331 print(str1 + str2 + str3) 

332 

333 if itn % 10 == 0: 

334 print() 

335 

336 if callback is not None: 

337 callback(x) 

338 

339 if istop != 0: 

340 break # TODO check this 

341 

342 if show: 

343 print() 

344 print(last + ' istop = %3g itn =%5g' % (istop,itn)) 

345 print(last + ' Anorm = %12.4e Acond = %12.4e' % (Anorm,Acond)) 

346 print(last + ' rnorm = %12.4e ynorm = %12.4e' % (rnorm,ynorm)) 

347 print(last + ' Arnorm = %12.4e' % (Arnorm,)) 

348 print(last + msg[istop+1]) 

349 

350 if istop == 6: 

351 info = maxiter 

352 else: 

353 info = 0 

354 

355 return (postprocess(x),info) 

356 

357 

358if __name__ == '__main__': 

359 from numpy import arange 

360 from scipy.sparse import spdiags 

361 

362 n = 10 

363 

364 residuals = [] 

365 

366 def cb(x): 

367 residuals.append(norm(b - A*x)) 

368 

369 # A = poisson((10,),format='csr') 

370 A = spdiags([arange(1,n+1,dtype=float)], [0], n, n, format='csr') 

371 M = spdiags([1.0/arange(1,n+1,dtype=float)], [0], n, n, format='csr') 

372 A.psolve = M.matvec 

373 b = zeros(A.shape[0]) 

374 x = minres(A,b,tol=1e-12,maxiter=None,callback=cb) 

375 # x = cg(A,b,x0=b,tol=1e-12,maxiter=None,callback=cb)[0]