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

1""" 

2Find a few eigenvectors and eigenvalues of a matrix. 

3 

4 

5Uses ARPACK: http://www.caam.rice.edu/software/ARPACK/ 

6 

7""" 

8# Wrapper implementation notes 

9# 

10# ARPACK Entry Points 

11# ------------------- 

12# The entry points to ARPACK are 

13# - (s,d)seupd : single and double precision symmetric matrix 

14# - (s,d,c,z)neupd: single,double,complex,double complex general matrix 

15# This wrapper puts the *neupd (general matrix) interfaces in eigs() 

16# and the *seupd (symmetric matrix) in eigsh(). 

17# There is no Hermitian complex/double complex interface. 

18# To find eigenvalues of a Hermitian matrix you 

19# must use eigs() and not eigsh() 

20# It might be desirable to handle the Hermitian case differently 

21# and, for example, return real eigenvalues. 

22 

23# Number of eigenvalues returned and complex eigenvalues 

24# ------------------------------------------------------ 

25# The ARPACK nonsymmetric real and double interface (s,d)naupd return 

26# eigenvalues and eigenvectors in real (float,double) arrays. 

27# Since the eigenvalues and eigenvectors are, in general, complex 

28# ARPACK puts the real and imaginary parts in consecutive entries 

29# in real-valued arrays. This wrapper puts the real entries 

30# into complex data types and attempts to return the requested eigenvalues 

31# and eigenvectors. 

32 

33 

34# Solver modes 

35# ------------ 

36# ARPACK and handle shifted and shift-inverse computations 

37# for eigenvalues by providing a shift (sigma) and a solver. 

38 

39__docformat__ = "restructuredtext en" 

40 

41__all__ = ['eigs', 'eigsh', 'svds', 'ArpackError', 'ArpackNoConvergence'] 

42 

43from . import _arpack 

44arpack_int = _arpack.timing.nbx.dtype 

45 

46import numpy as np 

47import warnings 

48from scipy.sparse.linalg.interface import aslinearoperator, LinearOperator 

49from scipy.sparse import eye, issparse, isspmatrix, isspmatrix_csr 

50from scipy.linalg import eig, eigh, lu_factor, lu_solve 

51from scipy.sparse.sputils import isdense, is_pydata_spmatrix 

52from scipy.sparse.linalg import gmres, splu 

53from scipy.sparse.linalg.eigen.lobpcg import lobpcg 

54from scipy._lib._util import _aligned_zeros 

55from scipy._lib._threadsafety import ReentrancyLock 

56 

57 

58_type_conv = {'f': 's', 'd': 'd', 'F': 'c', 'D': 'z'} 

59_ndigits = {'f': 5, 'd': 12, 'F': 5, 'D': 12} 

60 

61DNAUPD_ERRORS = { 

62 0: "Normal exit.", 

63 1: "Maximum number of iterations taken. " 

64 "All possible eigenvalues of OP has been found. IPARAM(5) " 

65 "returns the number of wanted converged Ritz values.", 

66 2: "No longer an informational error. Deprecated starting " 

67 "with release 2 of ARPACK.", 

68 3: "No shifts could be applied during a cycle of the " 

69 "Implicitly restarted Arnoldi iteration. One possibility " 

70 "is to increase the size of NCV relative to NEV. ", 

71 -1: "N must be positive.", 

72 -2: "NEV must be positive.", 

73 -3: "NCV-NEV >= 2 and less than or equal to N.", 

74 -4: "The maximum number of Arnoldi update iterations allowed " 

75 "must be greater than zero.", 

76 -5: " WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'", 

77 -6: "BMAT must be one of 'I' or 'G'.", 

78 -7: "Length of private work array WORKL is not sufficient.", 

79 -8: "Error return from LAPACK eigenvalue calculation;", 

80 -9: "Starting vector is zero.", 

81 -10: "IPARAM(7) must be 1,2,3,4.", 

82 -11: "IPARAM(7) = 1 and BMAT = 'G' are incompatible.", 

83 -12: "IPARAM(1) must be equal to 0 or 1.", 

84 -13: "NEV and WHICH = 'BE' are incompatible.", 

85 -9999: "Could not build an Arnoldi factorization. " 

86 "IPARAM(5) returns the size of the current Arnoldi " 

87 "factorization. The user is advised to check that " 

88 "enough workspace and array storage has been allocated." 

89} 

90 

91SNAUPD_ERRORS = DNAUPD_ERRORS 

92 

93ZNAUPD_ERRORS = DNAUPD_ERRORS.copy() 

94ZNAUPD_ERRORS[-10] = "IPARAM(7) must be 1,2,3." 

95 

96CNAUPD_ERRORS = ZNAUPD_ERRORS 

97 

98DSAUPD_ERRORS = { 

99 0: "Normal exit.", 

100 1: "Maximum number of iterations taken. " 

101 "All possible eigenvalues of OP has been found.", 

102 2: "No longer an informational error. Deprecated starting with " 

103 "release 2 of ARPACK.", 

104 3: "No shifts could be applied during a cycle of the Implicitly " 

105 "restarted Arnoldi iteration. One possibility is to increase " 

106 "the size of NCV relative to NEV. ", 

107 -1: "N must be positive.", 

108 -2: "NEV must be positive.", 

109 -3: "NCV must be greater than NEV and less than or equal to N.", 

110 -4: "The maximum number of Arnoldi update iterations allowed " 

111 "must be greater than zero.", 

112 -5: "WHICH must be one of 'LM', 'SM', 'LA', 'SA' or 'BE'.", 

113 -6: "BMAT must be one of 'I' or 'G'.", 

114 -7: "Length of private work array WORKL is not sufficient.", 

115 -8: "Error return from trid. eigenvalue calculation; " 

116 "Informational error from LAPACK routine dsteqr .", 

117 -9: "Starting vector is zero.", 

118 -10: "IPARAM(7) must be 1,2,3,4,5.", 

119 -11: "IPARAM(7) = 1 and BMAT = 'G' are incompatible.", 

120 -12: "IPARAM(1) must be equal to 0 or 1.", 

121 -13: "NEV and WHICH = 'BE' are incompatible. ", 

122 -9999: "Could not build an Arnoldi factorization. " 

123 "IPARAM(5) returns the size of the current Arnoldi " 

124 "factorization. The user is advised to check that " 

125 "enough workspace and array storage has been allocated.", 

126} 

127 

128SSAUPD_ERRORS = DSAUPD_ERRORS 

129 

130DNEUPD_ERRORS = { 

131 0: "Normal exit.", 

132 1: "The Schur form computed by LAPACK routine dlahqr " 

133 "could not be reordered by LAPACK routine dtrsen. " 

134 "Re-enter subroutine dneupd with IPARAM(5)NCV and " 

135 "increase the size of the arrays DR and DI to have " 

136 "dimension at least dimension NCV and allocate at least NCV " 

137 "columns for Z. NOTE: Not necessary if Z and V share " 

138 "the same space. Please notify the authors if this error" 

139 "occurs.", 

140 -1: "N must be positive.", 

141 -2: "NEV must be positive.", 

142 -3: "NCV-NEV >= 2 and less than or equal to N.", 

143 -5: "WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'", 

144 -6: "BMAT must be one of 'I' or 'G'.", 

145 -7: "Length of private work WORKL array is not sufficient.", 

146 -8: "Error return from calculation of a real Schur form. " 

147 "Informational error from LAPACK routine dlahqr .", 

148 -9: "Error return from calculation of eigenvectors. " 

149 "Informational error from LAPACK routine dtrevc.", 

150 -10: "IPARAM(7) must be 1,2,3,4.", 

151 -11: "IPARAM(7) = 1 and BMAT = 'G' are incompatible.", 

152 -12: "HOWMNY = 'S' not yet implemented", 

153 -13: "HOWMNY must be one of 'A' or 'P' if RVEC = .true.", 

154 -14: "DNAUPD did not find any eigenvalues to sufficient " 

155 "accuracy.", 

156 -15: "DNEUPD got a different count of the number of converged " 

157 "Ritz values than DNAUPD got. This indicates the user " 

158 "probably made an error in passing data from DNAUPD to " 

159 "DNEUPD or that the data was modified before entering " 

160 "DNEUPD", 

161} 

162 

163SNEUPD_ERRORS = DNEUPD_ERRORS.copy() 

164SNEUPD_ERRORS[1] = ("The Schur form computed by LAPACK routine slahqr " 

165 "could not be reordered by LAPACK routine strsen . " 

166 "Re-enter subroutine dneupd with IPARAM(5)=NCV and " 

167 "increase the size of the arrays DR and DI to have " 

168 "dimension at least dimension NCV and allocate at least " 

169 "NCV columns for Z. NOTE: Not necessary if Z and V share " 

170 "the same space. Please notify the authors if this error " 

171 "occurs.") 

172SNEUPD_ERRORS[-14] = ("SNAUPD did not find any eigenvalues to sufficient " 

173 "accuracy.") 

174SNEUPD_ERRORS[-15] = ("SNEUPD got a different count of the number of " 

175 "converged Ritz values than SNAUPD got. This indicates " 

176 "the user probably made an error in passing data from " 

177 "SNAUPD to SNEUPD or that the data was modified before " 

178 "entering SNEUPD") 

179 

180ZNEUPD_ERRORS = {0: "Normal exit.", 

181 1: "The Schur form computed by LAPACK routine csheqr " 

182 "could not be reordered by LAPACK routine ztrsen. " 

183 "Re-enter subroutine zneupd with IPARAM(5)=NCV and " 

184 "increase the size of the array D to have " 

185 "dimension at least dimension NCV and allocate at least " 

186 "NCV columns for Z. NOTE: Not necessary if Z and V share " 

187 "the same space. Please notify the authors if this error " 

188 "occurs.", 

189 -1: "N must be positive.", 

190 -2: "NEV must be positive.", 

191 -3: "NCV-NEV >= 1 and less than or equal to N.", 

192 -5: "WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'", 

193 -6: "BMAT must be one of 'I' or 'G'.", 

194 -7: "Length of private work WORKL array is not sufficient.", 

195 -8: "Error return from LAPACK eigenvalue calculation. " 

196 "This should never happened.", 

197 -9: "Error return from calculation of eigenvectors. " 

198 "Informational error from LAPACK routine ztrevc.", 

199 -10: "IPARAM(7) must be 1,2,3", 

200 -11: "IPARAM(7) = 1 and BMAT = 'G' are incompatible.", 

201 -12: "HOWMNY = 'S' not yet implemented", 

202 -13: "HOWMNY must be one of 'A' or 'P' if RVEC = .true.", 

203 -14: "ZNAUPD did not find any eigenvalues to sufficient " 

204 "accuracy.", 

205 -15: "ZNEUPD got a different count of the number of " 

206 "converged Ritz values than ZNAUPD got. This " 

207 "indicates the user probably made an error in passing " 

208 "data from ZNAUPD to ZNEUPD or that the data was " 

209 "modified before entering ZNEUPD" 

210 } 

211 

212CNEUPD_ERRORS = ZNEUPD_ERRORS.copy() 

213CNEUPD_ERRORS[-14] = ("CNAUPD did not find any eigenvalues to sufficient " 

214 "accuracy.") 

215CNEUPD_ERRORS[-15] = ("CNEUPD got a different count of the number of " 

216 "converged Ritz values than CNAUPD got. This indicates " 

217 "the user probably made an error in passing data from " 

218 "CNAUPD to CNEUPD or that the data was modified before " 

219 "entering CNEUPD") 

220 

221DSEUPD_ERRORS = { 

222 0: "Normal exit.", 

223 -1: "N must be positive.", 

224 -2: "NEV must be positive.", 

225 -3: "NCV must be greater than NEV and less than or equal to N.", 

226 -5: "WHICH must be one of 'LM', 'SM', 'LA', 'SA' or 'BE'.", 

227 -6: "BMAT must be one of 'I' or 'G'.", 

228 -7: "Length of private work WORKL array is not sufficient.", 

229 -8: ("Error return from trid. eigenvalue calculation; " 

230 "Information error from LAPACK routine dsteqr."), 

231 -9: "Starting vector is zero.", 

232 -10: "IPARAM(7) must be 1,2,3,4,5.", 

233 -11: "IPARAM(7) = 1 and BMAT = 'G' are incompatible.", 

234 -12: "NEV and WHICH = 'BE' are incompatible.", 

235 -14: "DSAUPD did not find any eigenvalues to sufficient accuracy.", 

236 -15: "HOWMNY must be one of 'A' or 'S' if RVEC = .true.", 

237 -16: "HOWMNY = 'S' not yet implemented", 

238 -17: ("DSEUPD got a different count of the number of converged " 

239 "Ritz values than DSAUPD got. This indicates the user " 

240 "probably made an error in passing data from DSAUPD to " 

241 "DSEUPD or that the data was modified before entering " 

242 "DSEUPD.") 

243} 

244 

245SSEUPD_ERRORS = DSEUPD_ERRORS.copy() 

246SSEUPD_ERRORS[-14] = ("SSAUPD did not find any eigenvalues " 

247 "to sufficient accuracy.") 

248SSEUPD_ERRORS[-17] = ("SSEUPD got a different count of the number of " 

249 "converged " 

250 "Ritz values than SSAUPD got. This indicates the user " 

251 "probably made an error in passing data from SSAUPD to " 

252 "SSEUPD or that the data was modified before entering " 

253 "SSEUPD.") 

254 

255_SAUPD_ERRORS = {'d': DSAUPD_ERRORS, 

256 's': SSAUPD_ERRORS} 

257_NAUPD_ERRORS = {'d': DNAUPD_ERRORS, 

258 's': SNAUPD_ERRORS, 

259 'z': ZNAUPD_ERRORS, 

260 'c': CNAUPD_ERRORS} 

261_SEUPD_ERRORS = {'d': DSEUPD_ERRORS, 

262 's': SSEUPD_ERRORS} 

263_NEUPD_ERRORS = {'d': DNEUPD_ERRORS, 

264 's': SNEUPD_ERRORS, 

265 'z': ZNEUPD_ERRORS, 

266 'c': CNEUPD_ERRORS} 

267 

268# accepted values of parameter WHICH in _SEUPD 

269_SEUPD_WHICH = ['LM', 'SM', 'LA', 'SA', 'BE'] 

270 

271# accepted values of parameter WHICH in _NAUPD 

272_NEUPD_WHICH = ['LM', 'SM', 'LR', 'SR', 'LI', 'SI'] 

273 

274 

275class ArpackError(RuntimeError): 

276 """ 

277 ARPACK error 

278 """ 

279 def __init__(self, info, infodict=_NAUPD_ERRORS): 

280 msg = infodict.get(info, "Unknown error") 

281 RuntimeError.__init__(self, "ARPACK error %d: %s" % (info, msg)) 

282 

283 

284class ArpackNoConvergence(ArpackError): 

285 """ 

286 ARPACK iteration did not converge 

287 

288 Attributes 

289 ---------- 

290 eigenvalues : ndarray 

291 Partial result. Converged eigenvalues. 

292 eigenvectors : ndarray 

293 Partial result. Converged eigenvectors. 

294 

295 """ 

296 def __init__(self, msg, eigenvalues, eigenvectors): 

297 ArpackError.__init__(self, -1, {-1: msg}) 

298 self.eigenvalues = eigenvalues 

299 self.eigenvectors = eigenvectors 

300 

301 

302def choose_ncv(k): 

303 """ 

304 Choose number of lanczos vectors based on target number 

305 of singular/eigen values and vectors to compute, k. 

306 """ 

307 return max(2 * k + 1, 20) 

308 

309 

310class _ArpackParams(object): 

311 def __init__(self, n, k, tp, mode=1, sigma=None, 

312 ncv=None, v0=None, maxiter=None, which="LM", tol=0): 

313 if k <= 0: 

314 raise ValueError("k must be positive, k=%d" % k) 

315 

316 if maxiter is None: 

317 maxiter = n * 10 

318 if maxiter <= 0: 

319 raise ValueError("maxiter must be positive, maxiter=%d" % maxiter) 

320 

321 if tp not in 'fdFD': 

322 raise ValueError("matrix type must be 'f', 'd', 'F', or 'D'") 

323 

324 if v0 is not None: 

325 # ARPACK overwrites its initial resid, make a copy 

326 self.resid = np.array(v0, copy=True) 

327 info = 1 

328 else: 

329 # ARPACK will use a random initial vector. 

330 self.resid = np.zeros(n, tp) 

331 info = 0 

332 

333 if sigma is None: 

334 #sigma not used 

335 self.sigma = 0 

336 else: 

337 self.sigma = sigma 

338 

339 if ncv is None: 

340 ncv = choose_ncv(k) 

341 ncv = min(ncv, n) 

342 

343 self.v = np.zeros((n, ncv), tp) # holds Ritz vectors 

344 self.iparam = np.zeros(11, arpack_int) 

345 

346 # set solver mode and parameters 

347 ishfts = 1 

348 self.mode = mode 

349 self.iparam[0] = ishfts 

350 self.iparam[2] = maxiter 

351 self.iparam[3] = 1 

352 self.iparam[6] = mode 

353 

354 self.n = n 

355 self.tol = tol 

356 self.k = k 

357 self.maxiter = maxiter 

358 self.ncv = ncv 

359 self.which = which 

360 self.tp = tp 

361 self.info = info 

362 

363 self.converged = False 

364 self.ido = 0 

365 

366 def _raise_no_convergence(self): 

367 msg = "No convergence (%d iterations, %d/%d eigenvectors converged)" 

368 k_ok = self.iparam[4] 

369 num_iter = self.iparam[2] 

370 try: 

371 ev, vec = self.extract(True) 

372 except ArpackError as err: 

373 msg = "%s [%s]" % (msg, err) 

374 ev = np.zeros((0,)) 

375 vec = np.zeros((self.n, 0)) 

376 k_ok = 0 

377 raise ArpackNoConvergence(msg % (num_iter, k_ok, self.k), ev, vec) 

378 

379 

380class _SymmetricArpackParams(_ArpackParams): 

381 def __init__(self, n, k, tp, matvec, mode=1, M_matvec=None, 

382 Minv_matvec=None, sigma=None, 

383 ncv=None, v0=None, maxiter=None, which="LM", tol=0): 

384 # The following modes are supported: 

385 # mode = 1: 

386 # Solve the standard eigenvalue problem: 

387 # A*x = lambda*x : 

388 # A - symmetric 

389 # Arguments should be 

390 # matvec = left multiplication by A 

391 # M_matvec = None [not used] 

392 # Minv_matvec = None [not used] 

393 # 

394 # mode = 2: 

395 # Solve the general eigenvalue problem: 

396 # A*x = lambda*M*x 

397 # A - symmetric 

398 # M - symmetric positive definite 

399 # Arguments should be 

400 # matvec = left multiplication by A 

401 # M_matvec = left multiplication by M 

402 # Minv_matvec = left multiplication by M^-1 

403 # 

404 # mode = 3: 

405 # Solve the general eigenvalue problem in shift-invert mode: 

406 # A*x = lambda*M*x 

407 # A - symmetric 

408 # M - symmetric positive semi-definite 

409 # Arguments should be 

410 # matvec = None [not used] 

411 # M_matvec = left multiplication by M 

412 # or None, if M is the identity 

413 # Minv_matvec = left multiplication by [A-sigma*M]^-1 

414 # 

415 # mode = 4: 

416 # Solve the general eigenvalue problem in Buckling mode: 

417 # A*x = lambda*AG*x 

418 # A - symmetric positive semi-definite 

419 # AG - symmetric indefinite 

420 # Arguments should be 

421 # matvec = left multiplication by A 

422 # M_matvec = None [not used] 

423 # Minv_matvec = left multiplication by [A-sigma*AG]^-1 

424 # 

425 # mode = 5: 

426 # Solve the general eigenvalue problem in Cayley-transformed mode: 

427 # A*x = lambda*M*x 

428 # A - symmetric 

429 # M - symmetric positive semi-definite 

430 # Arguments should be 

431 # matvec = left multiplication by A 

432 # M_matvec = left multiplication by M 

433 # or None, if M is the identity 

434 # Minv_matvec = left multiplication by [A-sigma*M]^-1 

435 if mode == 1: 

436 if matvec is None: 

437 raise ValueError("matvec must be specified for mode=1") 

438 if M_matvec is not None: 

439 raise ValueError("M_matvec cannot be specified for mode=1") 

440 if Minv_matvec is not None: 

441 raise ValueError("Minv_matvec cannot be specified for mode=1") 

442 

443 self.OP = matvec 

444 self.B = lambda x: x 

445 self.bmat = 'I' 

446 elif mode == 2: 

447 if matvec is None: 

448 raise ValueError("matvec must be specified for mode=2") 

449 if M_matvec is None: 

450 raise ValueError("M_matvec must be specified for mode=2") 

451 if Minv_matvec is None: 

452 raise ValueError("Minv_matvec must be specified for mode=2") 

453 

454 self.OP = lambda x: Minv_matvec(matvec(x)) 

455 self.OPa = Minv_matvec 

456 self.OPb = matvec 

457 self.B = M_matvec 

458 self.bmat = 'G' 

459 elif mode == 3: 

460 if matvec is not None: 

461 raise ValueError("matvec must not be specified for mode=3") 

462 if Minv_matvec is None: 

463 raise ValueError("Minv_matvec must be specified for mode=3") 

464 

465 if M_matvec is None: 

466 self.OP = Minv_matvec 

467 self.OPa = Minv_matvec 

468 self.B = lambda x: x 

469 self.bmat = 'I' 

470 else: 

471 self.OP = lambda x: Minv_matvec(M_matvec(x)) 

472 self.OPa = Minv_matvec 

473 self.B = M_matvec 

474 self.bmat = 'G' 

475 elif mode == 4: 

476 if matvec is None: 

477 raise ValueError("matvec must be specified for mode=4") 

478 if M_matvec is not None: 

479 raise ValueError("M_matvec must not be specified for mode=4") 

480 if Minv_matvec is None: 

481 raise ValueError("Minv_matvec must be specified for mode=4") 

482 self.OPa = Minv_matvec 

483 self.OP = lambda x: self.OPa(matvec(x)) 

484 self.B = matvec 

485 self.bmat = 'G' 

486 elif mode == 5: 

487 if matvec is None: 

488 raise ValueError("matvec must be specified for mode=5") 

489 if Minv_matvec is None: 

490 raise ValueError("Minv_matvec must be specified for mode=5") 

491 

492 self.OPa = Minv_matvec 

493 self.A_matvec = matvec 

494 

495 if M_matvec is None: 

496 self.OP = lambda x: Minv_matvec(matvec(x) + sigma * x) 

497 self.B = lambda x: x 

498 self.bmat = 'I' 

499 else: 

500 self.OP = lambda x: Minv_matvec(matvec(x) 

501 + sigma * M_matvec(x)) 

502 self.B = M_matvec 

503 self.bmat = 'G' 

504 else: 

505 raise ValueError("mode=%i not implemented" % mode) 

506 

507 if which not in _SEUPD_WHICH: 

508 raise ValueError("which must be one of %s" 

509 % ' '.join(_SEUPD_WHICH)) 

510 if k >= n: 

511 raise ValueError("k must be less than ndim(A), k=%d" % k) 

512 

513 _ArpackParams.__init__(self, n, k, tp, mode, sigma, 

514 ncv, v0, maxiter, which, tol) 

515 

516 if self.ncv > n or self.ncv <= k: 

517 raise ValueError("ncv must be k<ncv<=n, ncv=%s" % self.ncv) 

518 

519 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1 

520 self.workd = _aligned_zeros(3 * n, self.tp) 

521 self.workl = _aligned_zeros(self.ncv * (self.ncv + 8), self.tp) 

522 

523 ltr = _type_conv[self.tp] 

524 if ltr not in ["s", "d"]: 

525 raise ValueError("Input matrix is not real-valued.") 

526 

527 self._arpack_solver = _arpack.__dict__[ltr + 'saupd'] 

528 self._arpack_extract = _arpack.__dict__[ltr + 'seupd'] 

529 

530 self.iterate_infodict = _SAUPD_ERRORS[ltr] 

531 self.extract_infodict = _SEUPD_ERRORS[ltr] 

532 

533 self.ipntr = np.zeros(11, arpack_int) 

534 

535 def iterate(self): 

536 self.ido, self.tol, self.resid, self.v, self.iparam, self.ipntr, self.info = \ 

537 self._arpack_solver(self.ido, self.bmat, self.which, self.k, 

538 self.tol, self.resid, self.v, self.iparam, 

539 self.ipntr, self.workd, self.workl, self.info) 

540 

541 xslice = slice(self.ipntr[0] - 1, self.ipntr[0] - 1 + self.n) 

542 yslice = slice(self.ipntr[1] - 1, self.ipntr[1] - 1 + self.n) 

543 if self.ido == -1: 

544 # initialization 

545 self.workd[yslice] = self.OP(self.workd[xslice]) 

546 elif self.ido == 1: 

547 # compute y = Op*x 

548 if self.mode == 1: 

549 self.workd[yslice] = self.OP(self.workd[xslice]) 

550 elif self.mode == 2: 

551 self.workd[xslice] = self.OPb(self.workd[xslice]) 

552 self.workd[yslice] = self.OPa(self.workd[xslice]) 

553 elif self.mode == 5: 

554 Bxslice = slice(self.ipntr[2] - 1, self.ipntr[2] - 1 + self.n) 

555 Ax = self.A_matvec(self.workd[xslice]) 

556 self.workd[yslice] = self.OPa(Ax + (self.sigma * 

557 self.workd[Bxslice])) 

558 else: 

559 Bxslice = slice(self.ipntr[2] - 1, self.ipntr[2] - 1 + self.n) 

560 self.workd[yslice] = self.OPa(self.workd[Bxslice]) 

561 elif self.ido == 2: 

562 self.workd[yslice] = self.B(self.workd[xslice]) 

563 elif self.ido == 3: 

564 raise ValueError("ARPACK requested user shifts. Assure ISHIFT==0") 

565 else: 

566 self.converged = True 

567 

568 if self.info == 0: 

569 pass 

570 elif self.info == 1: 

571 self._raise_no_convergence() 

572 else: 

573 raise ArpackError(self.info, infodict=self.iterate_infodict) 

574 

575 def extract(self, return_eigenvectors): 

576 rvec = return_eigenvectors 

577 ierr = 0 

578 howmny = 'A' # return all eigenvectors 

579 sselect = np.zeros(self.ncv, 'int') # unused 

580 d, z, ierr = self._arpack_extract(rvec, howmny, sselect, self.sigma, 

581 self.bmat, self.which, self.k, 

582 self.tol, self.resid, self.v, 

583 self.iparam[0:7], self.ipntr, 

584 self.workd[0:2 * self.n], 

585 self.workl, ierr) 

586 if ierr != 0: 

587 raise ArpackError(ierr, infodict=self.extract_infodict) 

588 k_ok = self.iparam[4] 

589 d = d[:k_ok] 

590 z = z[:, :k_ok] 

591 

592 if return_eigenvectors: 

593 return d, z 

594 else: 

595 return d 

596 

597 

598class _UnsymmetricArpackParams(_ArpackParams): 

599 def __init__(self, n, k, tp, matvec, mode=1, M_matvec=None, 

600 Minv_matvec=None, sigma=None, 

601 ncv=None, v0=None, maxiter=None, which="LM", tol=0): 

602 # The following modes are supported: 

603 # mode = 1: 

604 # Solve the standard eigenvalue problem: 

605 # A*x = lambda*x 

606 # A - square matrix 

607 # Arguments should be 

608 # matvec = left multiplication by A 

609 # M_matvec = None [not used] 

610 # Minv_matvec = None [not used] 

611 # 

612 # mode = 2: 

613 # Solve the generalized eigenvalue problem: 

614 # A*x = lambda*M*x 

615 # A - square matrix 

616 # M - symmetric, positive semi-definite 

617 # Arguments should be 

618 # matvec = left multiplication by A 

619 # M_matvec = left multiplication by M 

620 # Minv_matvec = left multiplication by M^-1 

621 # 

622 # mode = 3,4: 

623 # Solve the general eigenvalue problem in shift-invert mode: 

624 # A*x = lambda*M*x 

625 # A - square matrix 

626 # M - symmetric, positive semi-definite 

627 # Arguments should be 

628 # matvec = None [not used] 

629 # M_matvec = left multiplication by M 

630 # or None, if M is the identity 

631 # Minv_matvec = left multiplication by [A-sigma*M]^-1 

632 # if A is real and mode==3, use the real part of Minv_matvec 

633 # if A is real and mode==4, use the imag part of Minv_matvec 

634 # if A is complex and mode==3, 

635 # use real and imag parts of Minv_matvec 

636 if mode == 1: 

637 if matvec is None: 

638 raise ValueError("matvec must be specified for mode=1") 

639 if M_matvec is not None: 

640 raise ValueError("M_matvec cannot be specified for mode=1") 

641 if Minv_matvec is not None: 

642 raise ValueError("Minv_matvec cannot be specified for mode=1") 

643 

644 self.OP = matvec 

645 self.B = lambda x: x 

646 self.bmat = 'I' 

647 elif mode == 2: 

648 if matvec is None: 

649 raise ValueError("matvec must be specified for mode=2") 

650 if M_matvec is None: 

651 raise ValueError("M_matvec must be specified for mode=2") 

652 if Minv_matvec is None: 

653 raise ValueError("Minv_matvec must be specified for mode=2") 

654 

655 self.OP = lambda x: Minv_matvec(matvec(x)) 

656 self.OPa = Minv_matvec 

657 self.OPb = matvec 

658 self.B = M_matvec 

659 self.bmat = 'G' 

660 elif mode in (3, 4): 

661 if matvec is None: 

662 raise ValueError("matvec must be specified " 

663 "for mode in (3,4)") 

664 if Minv_matvec is None: 

665 raise ValueError("Minv_matvec must be specified " 

666 "for mode in (3,4)") 

667 

668 self.matvec = matvec 

669 if tp in 'DF': # complex type 

670 if mode == 3: 

671 self.OPa = Minv_matvec 

672 else: 

673 raise ValueError("mode=4 invalid for complex A") 

674 else: # real type 

675 if mode == 3: 

676 self.OPa = lambda x: np.real(Minv_matvec(x)) 

677 else: 

678 self.OPa = lambda x: np.imag(Minv_matvec(x)) 

679 if M_matvec is None: 

680 self.B = lambda x: x 

681 self.bmat = 'I' 

682 self.OP = self.OPa 

683 else: 

684 self.B = M_matvec 

685 self.bmat = 'G' 

686 self.OP = lambda x: self.OPa(M_matvec(x)) 

687 else: 

688 raise ValueError("mode=%i not implemented" % mode) 

689 

690 if which not in _NEUPD_WHICH: 

691 raise ValueError("Parameter which must be one of %s" 

692 % ' '.join(_NEUPD_WHICH)) 

693 if k >= n - 1: 

694 raise ValueError("k must be less than ndim(A)-1, k=%d" % k) 

695 

696 _ArpackParams.__init__(self, n, k, tp, mode, sigma, 

697 ncv, v0, maxiter, which, tol) 

698 

699 if self.ncv > n or self.ncv <= k + 1: 

700 raise ValueError("ncv must be k+1<ncv<=n, ncv=%s" % self.ncv) 

701 

702 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1 

703 self.workd = _aligned_zeros(3 * n, self.tp) 

704 self.workl = _aligned_zeros(3 * self.ncv * (self.ncv + 2), self.tp) 

705 

706 ltr = _type_conv[self.tp] 

707 self._arpack_solver = _arpack.__dict__[ltr + 'naupd'] 

708 self._arpack_extract = _arpack.__dict__[ltr + 'neupd'] 

709 

710 self.iterate_infodict = _NAUPD_ERRORS[ltr] 

711 self.extract_infodict = _NEUPD_ERRORS[ltr] 

712 

713 self.ipntr = np.zeros(14, arpack_int) 

714 

715 if self.tp in 'FD': 

716 # Use _aligned_zeros to work around a f2py bug in Numpy 1.9.1 

717 self.rwork = _aligned_zeros(self.ncv, self.tp.lower()) 

718 else: 

719 self.rwork = None 

720 

721 def iterate(self): 

722 if self.tp in 'fd': 

723 self.ido, self.tol, self.resid, self.v, self.iparam, self.ipntr, self.info =\ 

724 self._arpack_solver(self.ido, self.bmat, self.which, self.k, 

725 self.tol, self.resid, self.v, self.iparam, 

726 self.ipntr, self.workd, self.workl, 

727 self.info) 

728 else: 

729 self.ido, self.tol, self.resid, self.v, self.iparam, self.ipntr, self.info =\ 

730 self._arpack_solver(self.ido, self.bmat, self.which, self.k, 

731 self.tol, self.resid, self.v, self.iparam, 

732 self.ipntr, self.workd, self.workl, 

733 self.rwork, self.info) 

734 

735 xslice = slice(self.ipntr[0] - 1, self.ipntr[0] - 1 + self.n) 

736 yslice = slice(self.ipntr[1] - 1, self.ipntr[1] - 1 + self.n) 

737 if self.ido == -1: 

738 # initialization 

739 self.workd[yslice] = self.OP(self.workd[xslice]) 

740 elif self.ido == 1: 

741 # compute y = Op*x 

742 if self.mode in (1, 2): 

743 self.workd[yslice] = self.OP(self.workd[xslice]) 

744 else: 

745 Bxslice = slice(self.ipntr[2] - 1, self.ipntr[2] - 1 + self.n) 

746 self.workd[yslice] = self.OPa(self.workd[Bxslice]) 

747 elif self.ido == 2: 

748 self.workd[yslice] = self.B(self.workd[xslice]) 

749 elif self.ido == 3: 

750 raise ValueError("ARPACK requested user shifts. Assure ISHIFT==0") 

751 else: 

752 self.converged = True 

753 

754 if self.info == 0: 

755 pass 

756 elif self.info == 1: 

757 self._raise_no_convergence() 

758 else: 

759 raise ArpackError(self.info, infodict=self.iterate_infodict) 

760 

761 def extract(self, return_eigenvectors): 

762 k, n = self.k, self.n 

763 

764 ierr = 0 

765 howmny = 'A' # return all eigenvectors 

766 sselect = np.zeros(self.ncv, 'int') # unused 

767 sigmar = np.real(self.sigma) 

768 sigmai = np.imag(self.sigma) 

769 workev = np.zeros(3 * self.ncv, self.tp) 

770 

771 if self.tp in 'fd': 

772 dr = np.zeros(k + 1, self.tp) 

773 di = np.zeros(k + 1, self.tp) 

774 zr = np.zeros((n, k + 1), self.tp) 

775 dr, di, zr, ierr = \ 

776 self._arpack_extract(return_eigenvectors, 

777 howmny, sselect, sigmar, sigmai, workev, 

778 self.bmat, self.which, k, self.tol, self.resid, 

779 self.v, self.iparam, self.ipntr, 

780 self.workd, self.workl, self.info) 

781 if ierr != 0: 

782 raise ArpackError(ierr, infodict=self.extract_infodict) 

783 nreturned = self.iparam[4] # number of good eigenvalues returned 

784 

785 # Build complex eigenvalues from real and imaginary parts 

786 d = dr + 1.0j * di 

787 

788 # Arrange the eigenvectors: complex eigenvectors are stored as 

789 # real,imaginary in consecutive columns 

790 z = zr.astype(self.tp.upper()) 

791 

792 # The ARPACK nonsymmetric real and double interface (s,d)naupd 

793 # return eigenvalues and eigenvectors in real (float,double) 

794 # arrays. 

795 

796 # Efficiency: this should check that return_eigenvectors == True 

797 # before going through this construction. 

798 if sigmai == 0: 

799 i = 0 

800 while i <= k: 

801 # check if complex 

802 if abs(d[i].imag) != 0: 

803 # this is a complex conjugate pair with eigenvalues 

804 # in consecutive columns 

805 if i < k: 

806 z[:, i] = zr[:, i] + 1.0j * zr[:, i + 1] 

807 z[:, i + 1] = z[:, i].conjugate() 

808 i += 1 

809 else: 

810 #last eigenvalue is complex: the imaginary part of 

811 # the eigenvector has not been returned 

812 #this can only happen if nreturned > k, so we'll 

813 # throw out this case. 

814 nreturned -= 1 

815 i += 1 

816 

817 else: 

818 # real matrix, mode 3 or 4, imag(sigma) is nonzero: 

819 # see remark 3 in <s,d>neupd.f 

820 # Build complex eigenvalues from real and imaginary parts 

821 i = 0 

822 while i <= k: 

823 if abs(d[i].imag) == 0: 

824 d[i] = np.dot(zr[:, i], self.matvec(zr[:, i])) 

825 else: 

826 if i < k: 

827 z[:, i] = zr[:, i] + 1.0j * zr[:, i + 1] 

828 z[:, i + 1] = z[:, i].conjugate() 

829 d[i] = ((np.dot(zr[:, i], 

830 self.matvec(zr[:, i])) 

831 + np.dot(zr[:, i + 1], 

832 self.matvec(zr[:, i + 1]))) 

833 + 1j * (np.dot(zr[:, i], 

834 self.matvec(zr[:, i + 1])) 

835 - np.dot(zr[:, i + 1], 

836 self.matvec(zr[:, i])))) 

837 d[i + 1] = d[i].conj() 

838 i += 1 

839 else: 

840 #last eigenvalue is complex: the imaginary part of 

841 # the eigenvector has not been returned 

842 #this can only happen if nreturned > k, so we'll 

843 # throw out this case. 

844 nreturned -= 1 

845 i += 1 

846 

847 # Now we have k+1 possible eigenvalues and eigenvectors 

848 # Return the ones specified by the keyword "which" 

849 

850 if nreturned <= k: 

851 # we got less or equal as many eigenvalues we wanted 

852 d = d[:nreturned] 

853 z = z[:, :nreturned] 

854 else: 

855 # we got one extra eigenvalue (likely a cc pair, but which?) 

856 if self.mode in (1, 2): 

857 rd = d 

858 elif self.mode in (3, 4): 

859 rd = 1 / (d - self.sigma) 

860 

861 if self.which in ['LR', 'SR']: 

862 ind = np.argsort(rd.real) 

863 elif self.which in ['LI', 'SI']: 

864 # for LI,SI ARPACK returns largest,smallest 

865 # abs(imaginary) (complex pairs come together) 

866 ind = np.argsort(abs(rd.imag)) 

867 else: 

868 ind = np.argsort(abs(rd)) 

869 

870 if self.which in ['LR', 'LM', 'LI']: 

871 ind = ind[-k:][::-1] 

872 elif self.which in ['SR', 'SM', 'SI']: 

873 ind = ind[:k] 

874 

875 d = d[ind] 

876 z = z[:, ind] 

877 else: 

878 # complex is so much simpler... 

879 d, z, ierr =\ 

880 self._arpack_extract(return_eigenvectors, 

881 howmny, sselect, self.sigma, workev, 

882 self.bmat, self.which, k, self.tol, self.resid, 

883 self.v, self.iparam, self.ipntr, 

884 self.workd, self.workl, self.rwork, ierr) 

885 

886 if ierr != 0: 

887 raise ArpackError(ierr, infodict=self.extract_infodict) 

888 

889 k_ok = self.iparam[4] 

890 d = d[:k_ok] 

891 z = z[:, :k_ok] 

892 

893 if return_eigenvectors: 

894 return d, z 

895 else: 

896 return d 

897 

898 

899def _aslinearoperator_with_dtype(m): 

900 m = aslinearoperator(m) 

901 if not hasattr(m, 'dtype'): 

902 x = np.zeros(m.shape[1]) 

903 m.dtype = (m * x).dtype 

904 return m 

905 

906 

907class SpLuInv(LinearOperator): 

908 """ 

909 SpLuInv: 

910 helper class to repeatedly solve M*x=b 

911 using a sparse LU-decomposition of M 

912 """ 

913 def __init__(self, M): 

914 self.M_lu = splu(M) 

915 self.shape = M.shape 

916 self.dtype = M.dtype 

917 self.isreal = not np.issubdtype(self.dtype, np.complexfloating) 

918 

919 def _matvec(self, x): 

920 # careful here: splu.solve will throw away imaginary 

921 # part of x if M is real 

922 x = np.asarray(x) 

923 if self.isreal and np.issubdtype(x.dtype, np.complexfloating): 

924 return (self.M_lu.solve(np.real(x).astype(self.dtype)) 

925 + 1j * self.M_lu.solve(np.imag(x).astype(self.dtype))) 

926 else: 

927 return self.M_lu.solve(x.astype(self.dtype)) 

928 

929 

930class LuInv(LinearOperator): 

931 """ 

932 LuInv: 

933 helper class to repeatedly solve M*x=b 

934 using an LU-decomposition of M 

935 """ 

936 def __init__(self, M): 

937 self.M_lu = lu_factor(M) 

938 self.shape = M.shape 

939 self.dtype = M.dtype 

940 

941 def _matvec(self, x): 

942 return lu_solve(self.M_lu, x) 

943 

944 

945def gmres_loose(A, b, tol): 

946 """ 

947 gmres with looser termination condition. 

948 """ 

949 b = np.asarray(b) 

950 min_tol = 1000 * np.sqrt(b.size) * np.finfo(b.dtype).eps 

951 return gmres(A, b, tol=max(tol, min_tol), atol=0) 

952 

953 

954class IterInv(LinearOperator): 

955 """ 

956 IterInv: 

957 helper class to repeatedly solve M*x=b 

958 using an iterative method. 

959 """ 

960 def __init__(self, M, ifunc=gmres_loose, tol=0): 

961 self.M = M 

962 if hasattr(M, 'dtype'): 

963 self.dtype = M.dtype 

964 else: 

965 x = np.zeros(M.shape[1]) 

966 self.dtype = (M * x).dtype 

967 self.shape = M.shape 

968 

969 if tol <= 0: 

970 # when tol=0, ARPACK uses machine tolerance as calculated 

971 # by LAPACK's _LAMCH function. We should match this 

972 tol = 2 * np.finfo(self.dtype).eps 

973 self.ifunc = ifunc 

974 self.tol = tol 

975 

976 def _matvec(self, x): 

977 b, info = self.ifunc(self.M, x, tol=self.tol) 

978 if info != 0: 

979 raise ValueError("Error in inverting M: function " 

980 "%s did not converge (info = %i)." 

981 % (self.ifunc.__name__, info)) 

982 return b 

983 

984 

985class IterOpInv(LinearOperator): 

986 """ 

987 IterOpInv: 

988 helper class to repeatedly solve [A-sigma*M]*x = b 

989 using an iterative method 

990 """ 

991 def __init__(self, A, M, sigma, ifunc=gmres_loose, tol=0): 

992 self.A = A 

993 self.M = M 

994 self.sigma = sigma 

995 

996 def mult_func(x): 

997 return A.matvec(x) - sigma * M.matvec(x) 

998 

999 def mult_func_M_None(x): 

1000 return A.matvec(x) - sigma * x 

1001 

1002 x = np.zeros(A.shape[1]) 

1003 if M is None: 

1004 dtype = mult_func_M_None(x).dtype 

1005 self.OP = LinearOperator(self.A.shape, 

1006 mult_func_M_None, 

1007 dtype=dtype) 

1008 else: 

1009 dtype = mult_func(x).dtype 

1010 self.OP = LinearOperator(self.A.shape, 

1011 mult_func, 

1012 dtype=dtype) 

1013 self.shape = A.shape 

1014 

1015 if tol <= 0: 

1016 # when tol=0, ARPACK uses machine tolerance as calculated 

1017 # by LAPACK's _LAMCH function. We should match this 

1018 tol = 2 * np.finfo(self.OP.dtype).eps 

1019 self.ifunc = ifunc 

1020 self.tol = tol 

1021 

1022 def _matvec(self, x): 

1023 b, info = self.ifunc(self.OP, x, tol=self.tol) 

1024 if info != 0: 

1025 raise ValueError("Error in inverting [A-sigma*M]: function " 

1026 "%s did not converge (info = %i)." 

1027 % (self.ifunc.__name__, info)) 

1028 return b 

1029 

1030 @property 

1031 def dtype(self): 

1032 return self.OP.dtype 

1033 

1034 

1035def _fast_spmatrix_to_csc(A, hermitian=False): 

1036 """Convert sparse matrix to CSC (by transposing, if possible)""" 

1037 if (isspmatrix_csr(A) and hermitian 

1038 and not np.issubdtype(A.dtype, np.complexfloating)): 

1039 return A.T 

1040 elif is_pydata_spmatrix(A): 

1041 # No need to convert 

1042 return A 

1043 else: 

1044 return A.tocsc() 

1045 

1046 

1047def get_inv_matvec(M, hermitian=False, tol=0): 

1048 if isdense(M): 

1049 return LuInv(M).matvec 

1050 elif isspmatrix(M) or is_pydata_spmatrix(M): 

1051 M = _fast_spmatrix_to_csc(M, hermitian=hermitian) 

1052 return SpLuInv(M).matvec 

1053 else: 

1054 return IterInv(M, tol=tol).matvec 

1055 

1056 

1057def get_OPinv_matvec(A, M, sigma, hermitian=False, tol=0): 

1058 if sigma == 0: 

1059 return get_inv_matvec(A, hermitian=hermitian, tol=tol) 

1060 

1061 if M is None: 

1062 #M is the identity matrix 

1063 if isdense(A): 

1064 if (np.issubdtype(A.dtype, np.complexfloating) 

1065 or np.imag(sigma) == 0): 

1066 A = np.copy(A) 

1067 else: 

1068 A = A + 0j 

1069 A.flat[::A.shape[1] + 1] -= sigma 

1070 return LuInv(A).matvec 

1071 elif isspmatrix(A) or is_pydata_spmatrix(A): 

1072 A = A - sigma * eye(A.shape[0]) 

1073 A = _fast_spmatrix_to_csc(A, hermitian=hermitian) 

1074 return SpLuInv(A).matvec 

1075 else: 

1076 return IterOpInv(_aslinearoperator_with_dtype(A), 

1077 M, sigma, tol=tol).matvec 

1078 else: 

1079 if ((not isdense(A) and not isspmatrix(A) and not is_pydata_spmatrix(A)) or 

1080 (not isdense(M) and not isspmatrix(M) and not is_pydata_spmatrix(A))): 

1081 return IterOpInv(_aslinearoperator_with_dtype(A), 

1082 _aslinearoperator_with_dtype(M), 

1083 sigma, tol=tol).matvec 

1084 elif isdense(A) or isdense(M): 

1085 return LuInv(A - sigma * M).matvec 

1086 else: 

1087 OP = A - sigma * M 

1088 OP = _fast_spmatrix_to_csc(OP, hermitian=hermitian) 

1089 return SpLuInv(OP).matvec 

1090 

1091 

1092# ARPACK is not threadsafe or reentrant (SAVE variables), so we need a 

1093# lock and a re-entering check. 

1094_ARPACK_LOCK = ReentrancyLock("Nested calls to eigs/eighs not allowed: " 

1095 "ARPACK is not re-entrant") 

1096 

1097 

1098def eigs(A, k=6, M=None, sigma=None, which='LM', v0=None, 

1099 ncv=None, maxiter=None, tol=0, return_eigenvectors=True, 

1100 Minv=None, OPinv=None, OPpart=None): 

1101 """ 

1102 Find k eigenvalues and eigenvectors of the square matrix A. 

1103 

1104 Solves ``A * x[i] = w[i] * x[i]``, the standard eigenvalue problem 

1105 for w[i] eigenvalues with corresponding eigenvectors x[i]. 

1106 

1107 If M is specified, solves ``A * x[i] = w[i] * M * x[i]``, the 

1108 generalized eigenvalue problem for w[i] eigenvalues 

1109 with corresponding eigenvectors x[i] 

1110 

1111 Parameters 

1112 ---------- 

1113 A : ndarray, sparse matrix or LinearOperator 

1114 An array, sparse matrix, or LinearOperator representing 

1115 the operation ``A * x``, where A is a real or complex square matrix. 

1116 k : int, optional 

1117 The number of eigenvalues and eigenvectors desired. 

1118 `k` must be smaller than N-1. It is not possible to compute all 

1119 eigenvectors of a matrix. 

1120 M : ndarray, sparse matrix or LinearOperator, optional 

1121 An array, sparse matrix, or LinearOperator representing 

1122 the operation M*x for the generalized eigenvalue problem 

1123 

1124 A * x = w * M * x. 

1125 

1126 M must represent a real, symmetric matrix if A is real, and must 

1127 represent a complex, hermitian matrix if A is complex. For best 

1128 results, the data type of M should be the same as that of A. 

1129 Additionally: 

1130 

1131 If `sigma` is None, M is positive definite 

1132 

1133 If sigma is specified, M is positive semi-definite 

1134 

1135 If sigma is None, eigs requires an operator to compute the solution 

1136 of the linear equation ``M * x = b``. This is done internally via a 

1137 (sparse) LU decomposition for an explicit matrix M, or via an 

1138 iterative solver for a general linear operator. Alternatively, 

1139 the user can supply the matrix or operator Minv, which gives 

1140 ``x = Minv * b = M^-1 * b``. 

1141 sigma : real or complex, optional 

1142 Find eigenvalues near sigma using shift-invert mode. This requires 

1143 an operator to compute the solution of the linear system 

1144 ``[A - sigma * M] * x = b``, where M is the identity matrix if 

1145 unspecified. This is computed internally via a (sparse) LU 

1146 decomposition for explicit matrices A & M, or via an iterative 

1147 solver if either A or M is a general linear operator. 

1148 Alternatively, the user can supply the matrix or operator OPinv, 

1149 which gives ``x = OPinv * b = [A - sigma * M]^-1 * b``. 

1150 For a real matrix A, shift-invert can either be done in imaginary 

1151 mode or real mode, specified by the parameter OPpart ('r' or 'i'). 

1152 Note that when sigma is specified, the keyword 'which' (below) 

1153 refers to the shifted eigenvalues ``w'[i]`` where: 

1154 

1155 If A is real and OPpart == 'r' (default), 

1156 ``w'[i] = 1/2 * [1/(w[i]-sigma) + 1/(w[i]-conj(sigma))]``. 

1157 

1158 If A is real and OPpart == 'i', 

1159 ``w'[i] = 1/2i * [1/(w[i]-sigma) - 1/(w[i]-conj(sigma))]``. 

1160 

1161 If A is complex, ``w'[i] = 1/(w[i]-sigma)``. 

1162 

1163 v0 : ndarray, optional 

1164 Starting vector for iteration. 

1165 Default: random 

1166 ncv : int, optional 

1167 The number of Lanczos vectors generated 

1168 `ncv` must be greater than `k`; it is recommended that ``ncv > 2*k``. 

1169 Default: ``min(n, max(2*k + 1, 20))`` 

1170 which : str, ['LM' | 'SM' | 'LR' | 'SR' | 'LI' | 'SI'], optional 

1171 Which `k` eigenvectors and eigenvalues to find: 

1172 

1173 'LM' : largest magnitude 

1174 

1175 'SM' : smallest magnitude 

1176 

1177 'LR' : largest real part 

1178 

1179 'SR' : smallest real part 

1180 

1181 'LI' : largest imaginary part 

1182 

1183 'SI' : smallest imaginary part 

1184 

1185 When sigma != None, 'which' refers to the shifted eigenvalues w'[i] 

1186 (see discussion in 'sigma', above). ARPACK is generally better 

1187 at finding large values than small values. If small eigenvalues are 

1188 desired, consider using shift-invert mode for better performance. 

1189 maxiter : int, optional 

1190 Maximum number of Arnoldi update iterations allowed 

1191 Default: ``n*10`` 

1192 tol : float, optional 

1193 Relative accuracy for eigenvalues (stopping criterion) 

1194 The default value of 0 implies machine precision. 

1195 return_eigenvectors : bool, optional 

1196 Return eigenvectors (True) in addition to eigenvalues 

1197 Minv : ndarray, sparse matrix or LinearOperator, optional 

1198 See notes in M, above. 

1199 OPinv : ndarray, sparse matrix or LinearOperator, optional 

1200 See notes in sigma, above. 

1201 OPpart : {'r' or 'i'}, optional 

1202 See notes in sigma, above 

1203 

1204 Returns 

1205 ------- 

1206 w : ndarray 

1207 Array of k eigenvalues. 

1208 v : ndarray 

1209 An array of `k` eigenvectors. 

1210 ``v[:, i]`` is the eigenvector corresponding to the eigenvalue w[i]. 

1211 

1212 Raises 

1213 ------ 

1214 ArpackNoConvergence 

1215 When the requested convergence is not obtained. 

1216 The currently converged eigenvalues and eigenvectors can be found 

1217 as ``eigenvalues`` and ``eigenvectors`` attributes of the exception 

1218 object. 

1219 

1220 See Also 

1221 -------- 

1222 eigsh : eigenvalues and eigenvectors for symmetric matrix A 

1223 svds : singular value decomposition for a matrix A 

1224 

1225 Notes 

1226 ----- 

1227 This function is a wrapper to the ARPACK [1]_ SNEUPD, DNEUPD, CNEUPD, 

1228 ZNEUPD, functions which use the Implicitly Restarted Arnoldi Method to 

1229 find the eigenvalues and eigenvectors [2]_. 

1230 

1231 References 

1232 ---------- 

1233 .. [1] ARPACK Software, http://www.caam.rice.edu/software/ARPACK/ 

1234 .. [2] R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK USERS GUIDE: 

1235 Solution of Large Scale Eigenvalue Problems by Implicitly Restarted 

1236 Arnoldi Methods. SIAM, Philadelphia, PA, 1998. 

1237 

1238 Examples 

1239 -------- 

1240 Find 6 eigenvectors of the identity matrix: 

1241 

1242 >>> from scipy.sparse.linalg import eigs 

1243 >>> id = np.eye(13) 

1244 >>> vals, vecs = eigs(id, k=6) 

1245 >>> vals 

1246 array([ 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j]) 

1247 >>> vecs.shape 

1248 (13, 6) 

1249 

1250 """ 

1251 if A.shape[0] != A.shape[1]: 

1252 raise ValueError('expected square matrix (shape=%s)' % (A.shape,)) 

1253 if M is not None: 

1254 if M.shape != A.shape: 

1255 raise ValueError('wrong M dimensions %s, should be %s' 

1256 % (M.shape, A.shape)) 

1257 if np.dtype(M.dtype).char.lower() != np.dtype(A.dtype).char.lower(): 

1258 warnings.warn('M does not have the same type precision as A. ' 

1259 'This may adversely affect ARPACK convergence') 

1260 

1261 n = A.shape[0] 

1262 

1263 if k <= 0: 

1264 raise ValueError("k=%d must be greater than 0." % k) 

1265 

1266 if k >= n - 1: 

1267 warnings.warn("k >= N - 1 for N * N square matrix. " 

1268 "Attempting to use scipy.linalg.eig instead.", 

1269 RuntimeWarning) 

1270 

1271 if issparse(A): 

1272 raise TypeError("Cannot use scipy.linalg.eig for sparse A with " 

1273 "k >= N - 1. Use scipy.linalg.eig(A.toarray()) or" 

1274 " reduce k.") 

1275 if isinstance(A, LinearOperator): 

1276 raise TypeError("Cannot use scipy.linalg.eig for LinearOperator " 

1277 "A with k >= N - 1.") 

1278 if isinstance(M, LinearOperator): 

1279 raise TypeError("Cannot use scipy.linalg.eig for LinearOperator " 

1280 "M with k >= N - 1.") 

1281 

1282 return eig(A, b=M, right=return_eigenvectors) 

1283 

1284 if sigma is None: 

1285 matvec = _aslinearoperator_with_dtype(A).matvec 

1286 

1287 if OPinv is not None: 

1288 raise ValueError("OPinv should not be specified " 

1289 "with sigma = None.") 

1290 if OPpart is not None: 

1291 raise ValueError("OPpart should not be specified with " 

1292 "sigma = None or complex A") 

1293 

1294 if M is None: 

1295 #standard eigenvalue problem 

1296 mode = 1 

1297 M_matvec = None 

1298 Minv_matvec = None 

1299 if Minv is not None: 

1300 raise ValueError("Minv should not be " 

1301 "specified with M = None.") 

1302 else: 

1303 #general eigenvalue problem 

1304 mode = 2 

1305 if Minv is None: 

1306 Minv_matvec = get_inv_matvec(M, hermitian=True, tol=tol) 

1307 else: 

1308 Minv = _aslinearoperator_with_dtype(Minv) 

1309 Minv_matvec = Minv.matvec 

1310 M_matvec = _aslinearoperator_with_dtype(M).matvec 

1311 else: 

1312 #sigma is not None: shift-invert mode 

1313 if np.issubdtype(A.dtype, np.complexfloating): 

1314 if OPpart is not None: 

1315 raise ValueError("OPpart should not be specified " 

1316 "with sigma=None or complex A") 

1317 mode = 3 

1318 elif OPpart is None or OPpart.lower() == 'r': 

1319 mode = 3 

1320 elif OPpart.lower() == 'i': 

1321 if np.imag(sigma) == 0: 

1322 raise ValueError("OPpart cannot be 'i' if sigma is real") 

1323 mode = 4 

1324 else: 

1325 raise ValueError("OPpart must be one of ('r','i')") 

1326 

1327 matvec = _aslinearoperator_with_dtype(A).matvec 

1328 if Minv is not None: 

1329 raise ValueError("Minv should not be specified when sigma is") 

1330 if OPinv is None: 

1331 Minv_matvec = get_OPinv_matvec(A, M, sigma, 

1332 hermitian=False, tol=tol) 

1333 else: 

1334 OPinv = _aslinearoperator_with_dtype(OPinv) 

1335 Minv_matvec = OPinv.matvec 

1336 if M is None: 

1337 M_matvec = None 

1338 else: 

1339 M_matvec = _aslinearoperator_with_dtype(M).matvec 

1340 

1341 params = _UnsymmetricArpackParams(n, k, A.dtype.char, matvec, mode, 

1342 M_matvec, Minv_matvec, sigma, 

1343 ncv, v0, maxiter, which, tol) 

1344 

1345 with _ARPACK_LOCK: 

1346 while not params.converged: 

1347 params.iterate() 

1348 

1349 return params.extract(return_eigenvectors) 

1350 

1351 

1352def eigsh(A, k=6, M=None, sigma=None, which='LM', v0=None, 

1353 ncv=None, maxiter=None, tol=0, return_eigenvectors=True, 

1354 Minv=None, OPinv=None, mode='normal'): 

1355 """ 

1356 Find k eigenvalues and eigenvectors of the real symmetric square matrix 

1357 or complex hermitian matrix A. 

1358 

1359 Solves ``A * x[i] = w[i] * x[i]``, the standard eigenvalue problem for 

1360 w[i] eigenvalues with corresponding eigenvectors x[i]. 

1361 

1362 If M is specified, solves ``A * x[i] = w[i] * M * x[i]``, the 

1363 generalized eigenvalue problem for w[i] eigenvalues 

1364 with corresponding eigenvectors x[i]. 

1365 

1366 Parameters 

1367 ---------- 

1368 A : ndarray, sparse matrix or LinearOperator 

1369 A square operator representing the operation ``A * x``, where ``A`` is 

1370 real symmetric or complex hermitian. For buckling mode (see below) 

1371 ``A`` must additionally be positive-definite. 

1372 k : int, optional 

1373 The number of eigenvalues and eigenvectors desired. 

1374 `k` must be smaller than N. It is not possible to compute all 

1375 eigenvectors of a matrix. 

1376 

1377 Returns 

1378 ------- 

1379 w : array 

1380 Array of k eigenvalues. 

1381 v : array 

1382 An array representing the `k` eigenvectors. The column ``v[:, i]`` is 

1383 the eigenvector corresponding to the eigenvalue ``w[i]``. 

1384 

1385 Other Parameters 

1386 ---------------- 

1387 M : An N x N matrix, array, sparse matrix, or linear operator representing 

1388 the operation ``M @ x`` for the generalized eigenvalue problem 

1389 

1390 A @ x = w * M @ x. 

1391 

1392 M must represent a real, symmetric matrix if A is real, and must 

1393 represent a complex, hermitian matrix if A is complex. For best 

1394 results, the data type of M should be the same as that of A. 

1395 Additionally: 

1396 

1397 If sigma is None, M is symmetric positive definite. 

1398 

1399 If sigma is specified, M is symmetric positive semi-definite. 

1400 

1401 In buckling mode, M is symmetric indefinite. 

1402 

1403 If sigma is None, eigsh requires an operator to compute the solution 

1404 of the linear equation ``M @ x = b``. This is done internally via a 

1405 (sparse) LU decomposition for an explicit matrix M, or via an 

1406 iterative solver for a general linear operator. Alternatively, 

1407 the user can supply the matrix or operator Minv, which gives 

1408 ``x = Minv @ b = M^-1 @ b``. 

1409 sigma : real 

1410 Find eigenvalues near sigma using shift-invert mode. This requires 

1411 an operator to compute the solution of the linear system 

1412 ``[A - sigma * M] x = b``, where M is the identity matrix if 

1413 unspecified. This is computed internally via a (sparse) LU 

1414 decomposition for explicit matrices A & M, or via an iterative 

1415 solver if either A or M is a general linear operator. 

1416 Alternatively, the user can supply the matrix or operator OPinv, 

1417 which gives ``x = OPinv @ b = [A - sigma * M]^-1 @ b``. 

1418 Note that when sigma is specified, the keyword 'which' refers to 

1419 the shifted eigenvalues ``w'[i]`` where: 

1420 

1421 if mode == 'normal', ``w'[i] = 1 / (w[i] - sigma)``. 

1422 

1423 if mode == 'cayley', ``w'[i] = (w[i] + sigma) / (w[i] - sigma)``. 

1424 

1425 if mode == 'buckling', ``w'[i] = w[i] / (w[i] - sigma)``. 

1426 

1427 (see further discussion in 'mode' below) 

1428 v0 : ndarray, optional 

1429 Starting vector for iteration. 

1430 Default: random 

1431 ncv : int, optional 

1432 The number of Lanczos vectors generated ncv must be greater than k and 

1433 smaller than n; it is recommended that ``ncv > 2*k``. 

1434 Default: ``min(n, max(2*k + 1, 20))`` 

1435 which : str ['LM' | 'SM' | 'LA' | 'SA' | 'BE'] 

1436 If A is a complex hermitian matrix, 'BE' is invalid. 

1437 Which `k` eigenvectors and eigenvalues to find: 

1438 

1439 'LM' : Largest (in magnitude) eigenvalues. 

1440 

1441 'SM' : Smallest (in magnitude) eigenvalues. 

1442 

1443 'LA' : Largest (algebraic) eigenvalues. 

1444 

1445 'SA' : Smallest (algebraic) eigenvalues. 

1446 

1447 'BE' : Half (k/2) from each end of the spectrum. 

1448 

1449 When k is odd, return one more (k/2+1) from the high end. 

1450 When sigma != None, 'which' refers to the shifted eigenvalues ``w'[i]`` 

1451 (see discussion in 'sigma', above). ARPACK is generally better 

1452 at finding large values than small values. If small eigenvalues are 

1453 desired, consider using shift-invert mode for better performance. 

1454 maxiter : int, optional 

1455 Maximum number of Arnoldi update iterations allowed. 

1456 Default: ``n*10`` 

1457 tol : float 

1458 Relative accuracy for eigenvalues (stopping criterion). 

1459 The default value of 0 implies machine precision. 

1460 Minv : N x N matrix, array, sparse matrix, or LinearOperator 

1461 See notes in M, above. 

1462 OPinv : N x N matrix, array, sparse matrix, or LinearOperator 

1463 See notes in sigma, above. 

1464 return_eigenvectors : bool 

1465 Return eigenvectors (True) in addition to eigenvalues. 

1466 This value determines the order in which eigenvalues are sorted. 

1467 The sort order is also dependent on the `which` variable. 

1468 

1469 For which = 'LM' or 'SA': 

1470 If `return_eigenvectors` is True, eigenvalues are sorted by 

1471 algebraic value. 

1472 

1473 If `return_eigenvectors` is False, eigenvalues are sorted by 

1474 absolute value. 

1475 

1476 For which = 'BE' or 'LA': 

1477 eigenvalues are always sorted by algebraic value. 

1478 

1479 For which = 'SM': 

1480 If `return_eigenvectors` is True, eigenvalues are sorted by 

1481 algebraic value. 

1482 

1483 If `return_eigenvectors` is False, eigenvalues are sorted by 

1484 decreasing absolute value. 

1485 

1486 mode : string ['normal' | 'buckling' | 'cayley'] 

1487 Specify strategy to use for shift-invert mode. This argument applies 

1488 only for real-valued A and sigma != None. For shift-invert mode, 

1489 ARPACK internally solves the eigenvalue problem 

1490 ``OP * x'[i] = w'[i] * B * x'[i]`` 

1491 and transforms the resulting Ritz vectors x'[i] and Ritz values w'[i] 

1492 into the desired eigenvectors and eigenvalues of the problem 

1493 ``A * x[i] = w[i] * M * x[i]``. 

1494 The modes are as follows: 

1495 

1496 'normal' : 

1497 OP = [A - sigma * M]^-1 @ M, 

1498 B = M, 

1499 w'[i] = 1 / (w[i] - sigma) 

1500 

1501 'buckling' : 

1502 OP = [A - sigma * M]^-1 @ A, 

1503 B = A, 

1504 w'[i] = w[i] / (w[i] - sigma) 

1505 

1506 'cayley' : 

1507 OP = [A - sigma * M]^-1 @ [A + sigma * M], 

1508 B = M, 

1509 w'[i] = (w[i] + sigma) / (w[i] - sigma) 

1510 

1511 The choice of mode will affect which eigenvalues are selected by 

1512 the keyword 'which', and can also impact the stability of 

1513 convergence (see [2] for a discussion). 

1514 

1515 Raises 

1516 ------ 

1517 ArpackNoConvergence 

1518 When the requested convergence is not obtained. 

1519 

1520 The currently converged eigenvalues and eigenvectors can be found 

1521 as ``eigenvalues`` and ``eigenvectors`` attributes of the exception 

1522 object. 

1523 

1524 See Also 

1525 -------- 

1526 eigs : eigenvalues and eigenvectors for a general (nonsymmetric) matrix A 

1527 svds : singular value decomposition for a matrix A 

1528 

1529 Notes 

1530 ----- 

1531 This function is a wrapper to the ARPACK [1]_ SSEUPD and DSEUPD 

1532 functions which use the Implicitly Restarted Lanczos Method to 

1533 find the eigenvalues and eigenvectors [2]_. 

1534 

1535 References 

1536 ---------- 

1537 .. [1] ARPACK Software, http://www.caam.rice.edu/software/ARPACK/ 

1538 .. [2] R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK USERS GUIDE: 

1539 Solution of Large Scale Eigenvalue Problems by Implicitly Restarted 

1540 Arnoldi Methods. SIAM, Philadelphia, PA, 1998. 

1541 

1542 Examples 

1543 -------- 

1544 >>> from scipy.sparse.linalg import eigsh 

1545 >>> identity = np.eye(13) 

1546 >>> eigenvalues, eigenvectors = eigsh(identity, k=6) 

1547 >>> eigenvalues 

1548 array([1., 1., 1., 1., 1., 1.]) 

1549 >>> eigenvectors.shape 

1550 (13, 6) 

1551 

1552 """ 

1553 # complex hermitian matrices should be solved with eigs 

1554 if np.issubdtype(A.dtype, np.complexfloating): 

1555 if mode != 'normal': 

1556 raise ValueError("mode=%s cannot be used with " 

1557 "complex matrix A" % mode) 

1558 if which == 'BE': 

1559 raise ValueError("which='BE' cannot be used with complex matrix A") 

1560 elif which == 'LA': 

1561 which = 'LR' 

1562 elif which == 'SA': 

1563 which = 'SR' 

1564 ret = eigs(A, k, M=M, sigma=sigma, which=which, v0=v0, 

1565 ncv=ncv, maxiter=maxiter, tol=tol, 

1566 return_eigenvectors=return_eigenvectors, Minv=Minv, 

1567 OPinv=OPinv) 

1568 

1569 if return_eigenvectors: 

1570 return ret[0].real, ret[1] 

1571 else: 

1572 return ret.real 

1573 

1574 if A.shape[0] != A.shape[1]: 

1575 raise ValueError('expected square matrix (shape=%s)' % (A.shape,)) 

1576 if M is not None: 

1577 if M.shape != A.shape: 

1578 raise ValueError('wrong M dimensions %s, should be %s' 

1579 % (M.shape, A.shape)) 

1580 if np.dtype(M.dtype).char.lower() != np.dtype(A.dtype).char.lower(): 

1581 warnings.warn('M does not have the same type precision as A. ' 

1582 'This may adversely affect ARPACK convergence') 

1583 

1584 n = A.shape[0] 

1585 

1586 if k <= 0: 

1587 raise ValueError("k must be greater than 0.") 

1588 

1589 if k >= n: 

1590 warnings.warn("k >= N for N * N square matrix. " 

1591 "Attempting to use scipy.linalg.eigh instead.", 

1592 RuntimeWarning) 

1593 

1594 if issparse(A): 

1595 raise TypeError("Cannot use scipy.linalg.eigh for sparse A with " 

1596 "k >= N. Use scipy.linalg.eigh(A.toarray()) or" 

1597 " reduce k.") 

1598 if isinstance(A, LinearOperator): 

1599 raise TypeError("Cannot use scipy.linalg.eigh for LinearOperator " 

1600 "A with k >= N.") 

1601 if isinstance(M, LinearOperator): 

1602 raise TypeError("Cannot use scipy.linalg.eigh for LinearOperator " 

1603 "M with k >= N.") 

1604 

1605 return eigh(A, b=M, eigvals_only=not return_eigenvectors) 

1606 

1607 if sigma is None: 

1608 A = _aslinearoperator_with_dtype(A) 

1609 matvec = A.matvec 

1610 

1611 if OPinv is not None: 

1612 raise ValueError("OPinv should not be specified " 

1613 "with sigma = None.") 

1614 if M is None: 

1615 #standard eigenvalue problem 

1616 mode = 1 

1617 M_matvec = None 

1618 Minv_matvec = None 

1619 if Minv is not None: 

1620 raise ValueError("Minv should not be " 

1621 "specified with M = None.") 

1622 else: 

1623 #general eigenvalue problem 

1624 mode = 2 

1625 if Minv is None: 

1626 Minv_matvec = get_inv_matvec(M, hermitian=True, tol=tol) 

1627 else: 

1628 Minv = _aslinearoperator_with_dtype(Minv) 

1629 Minv_matvec = Minv.matvec 

1630 M_matvec = _aslinearoperator_with_dtype(M).matvec 

1631 else: 

1632 # sigma is not None: shift-invert mode 

1633 if Minv is not None: 

1634 raise ValueError("Minv should not be specified when sigma is") 

1635 

1636 # normal mode 

1637 if mode == 'normal': 

1638 mode = 3 

1639 matvec = None 

1640 if OPinv is None: 

1641 Minv_matvec = get_OPinv_matvec(A, M, sigma, 

1642 hermitian=True, tol=tol) 

1643 else: 

1644 OPinv = _aslinearoperator_with_dtype(OPinv) 

1645 Minv_matvec = OPinv.matvec 

1646 if M is None: 

1647 M_matvec = None 

1648 else: 

1649 M = _aslinearoperator_with_dtype(M) 

1650 M_matvec = M.matvec 

1651 

1652 # buckling mode 

1653 elif mode == 'buckling': 

1654 mode = 4 

1655 if OPinv is None: 

1656 Minv_matvec = get_OPinv_matvec(A, M, sigma, 

1657 hermitian=True, tol=tol) 

1658 else: 

1659 Minv_matvec = _aslinearoperator_with_dtype(OPinv).matvec 

1660 matvec = _aslinearoperator_with_dtype(A).matvec 

1661 M_matvec = None 

1662 

1663 # cayley-transform mode 

1664 elif mode == 'cayley': 

1665 mode = 5 

1666 matvec = _aslinearoperator_with_dtype(A).matvec 

1667 if OPinv is None: 

1668 Minv_matvec = get_OPinv_matvec(A, M, sigma, 

1669 hermitian=True, tol=tol) 

1670 else: 

1671 Minv_matvec = _aslinearoperator_with_dtype(OPinv).matvec 

1672 if M is None: 

1673 M_matvec = None 

1674 else: 

1675 M_matvec = _aslinearoperator_with_dtype(M).matvec 

1676 

1677 # unrecognized mode 

1678 else: 

1679 raise ValueError("unrecognized mode '%s'" % mode) 

1680 

1681 params = _SymmetricArpackParams(n, k, A.dtype.char, matvec, mode, 

1682 M_matvec, Minv_matvec, sigma, 

1683 ncv, v0, maxiter, which, tol) 

1684 

1685 with _ARPACK_LOCK: 

1686 while not params.converged: 

1687 params.iterate() 

1688 

1689 return params.extract(return_eigenvectors) 

1690 

1691 

1692def _augmented_orthonormal_cols(x, k): 

1693 # extract the shape of the x array 

1694 n, m = x.shape 

1695 # create the expanded array and copy x into it 

1696 y = np.empty((n, m+k), dtype=x.dtype) 

1697 y[:, :m] = x 

1698 # do some modified gram schmidt to add k random orthonormal vectors 

1699 for i in range(k): 

1700 # sample a random initial vector 

1701 v = np.random.randn(n) 

1702 if np.iscomplexobj(x): 

1703 v = v + 1j*np.random.randn(n) 

1704 # subtract projections onto the existing unit length vectors 

1705 for j in range(m+i): 

1706 u = y[:, j] 

1707 v -= (np.dot(v, u.conj()) / np.dot(u, u.conj())) * u 

1708 # normalize v 

1709 v /= np.sqrt(np.dot(v, v.conj())) 

1710 # add v into the output array 

1711 y[:, m+i] = v 

1712 # return the expanded array 

1713 return y 

1714 

1715 

1716def _augmented_orthonormal_rows(x, k): 

1717 return _augmented_orthonormal_cols(x.T, k).T 

1718 

1719 

1720def _herm(x): 

1721 return x.T.conj() 

1722 

1723 

1724def svds(A, k=6, ncv=None, tol=0, which='LM', v0=None, 

1725 maxiter=None, return_singular_vectors=True, 

1726 solver='arpack'): 

1727 """Compute the largest or smallest k singular values/vectors for a sparse matrix. The order of the singular values is not guaranteed. 

1728 

1729 Parameters 

1730 ---------- 

1731 A : {sparse matrix, LinearOperator} 

1732 Array to compute the SVD on, of shape (M, N) 

1733 k : int, optional 

1734 Number of singular values and vectors to compute. 

1735 Must be 1 <= k < min(A.shape). 

1736 ncv : int, optional 

1737 The number of Lanczos vectors generated 

1738 ncv must be greater than k+1 and smaller than n; 

1739 it is recommended that ncv > 2*k 

1740 Default: ``min(n, max(2*k + 1, 20))`` 

1741 tol : float, optional 

1742 Tolerance for singular values. Zero (default) means machine precision. 

1743 which : str, ['LM' | 'SM'], optional 

1744 Which `k` singular values to find: 

1745 

1746 - 'LM' : largest singular values 

1747 - 'SM' : smallest singular values 

1748 

1749 .. versionadded:: 0.12.0 

1750 v0 : ndarray, optional 

1751 Starting vector for iteration, of length min(A.shape). Should be an 

1752 (approximate) left singular vector if N > M and a right singular 

1753 vector otherwise. 

1754 Default: random 

1755 

1756 .. versionadded:: 0.12.0 

1757 maxiter : int, optional 

1758 Maximum number of iterations. 

1759 

1760 .. versionadded:: 0.12.0 

1761 return_singular_vectors : bool or str, optional 

1762 - True: return singular vectors (True) in addition to singular values. 

1763 

1764 .. versionadded:: 0.12.0 

1765 

1766 - "u": only return the u matrix, without computing vh (if N > M). 

1767 - "vh": only return the vh matrix, without computing u (if N <= M). 

1768 

1769 .. versionadded:: 0.16.0 

1770 solver : str, optional 

1771 Eigenvalue solver to use. Should be 'arpack' or 'lobpcg'. 

1772 Default: 'arpack' 

1773 

1774 Returns 

1775 ------- 

1776 u : ndarray, shape=(M, k) 

1777 Unitary matrix having left singular vectors as columns. 

1778 If `return_singular_vectors` is "vh", this variable is not computed, 

1779 and None is returned instead. 

1780 s : ndarray, shape=(k,) 

1781 The singular values. 

1782 vt : ndarray, shape=(k, N) 

1783 Unitary matrix having right singular vectors as rows. 

1784 If `return_singular_vectors` is "u", this variable is not computed, 

1785 and None is returned instead. 

1786 

1787 

1788 Notes 

1789 ----- 

1790 This is a naive implementation using ARPACK or LOBPCG as an eigensolver 

1791 on A.H * A or A * A.H, depending on which one is more efficient. 

1792 

1793 Examples 

1794 -------- 

1795 >>> from scipy.sparse import csc_matrix 

1796 >>> from scipy.sparse.linalg import svds, eigs 

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

1798 >>> u, s, vt = svds(A, k=2) 

1799 >>> s 

1800 array([ 2.75193379, 5.6059665 ]) 

1801 >>> np.sqrt(eigs(A.dot(A.T), k=2)[0]).real 

1802 array([ 5.6059665 , 2.75193379]) 

1803 """ 

1804 if which == 'LM': 

1805 largest = True 

1806 elif which == 'SM': 

1807 largest = False 

1808 else: 

1809 raise ValueError("which must be either 'LM' or 'SM'.") 

1810 

1811 if not (isinstance(A, LinearOperator) or isspmatrix(A) or is_pydata_spmatrix(A)): 

1812 A = np.asarray(A) 

1813 

1814 n, m = A.shape 

1815 

1816 if k <= 0 or k >= min(n, m): 

1817 raise ValueError("k must be between 1 and min(A.shape), k=%d" % k) 

1818 

1819 if isinstance(A, LinearOperator): 

1820 if n > m: 

1821 X_dot = A.matvec 

1822 X_matmat = A.matmat 

1823 XH_dot = A.rmatvec 

1824 XH_mat = A.rmatmat 

1825 else: 

1826 X_dot = A.rmatvec 

1827 X_matmat = A.rmatmat 

1828 XH_dot = A.matvec 

1829 XH_mat = A.matmat 

1830 

1831 dtype = getattr(A, 'dtype', None) 

1832 if dtype is None: 

1833 dtype = A.dot(np.zeros([m, 1])).dtype 

1834 

1835 else: 

1836 if n > m: 

1837 X_dot = X_matmat = A.dot 

1838 XH_dot = XH_mat = _herm(A).dot 

1839 else: 

1840 XH_dot = XH_mat = A.dot 

1841 X_dot = X_matmat = _herm(A).dot 

1842 

1843 def matvec_XH_X(x): 

1844 return XH_dot(X_dot(x)) 

1845 

1846 def matmat_XH_X(x): 

1847 return XH_mat(X_matmat(x)) 

1848 

1849 XH_X = LinearOperator(matvec=matvec_XH_X, dtype=A.dtype, 

1850 matmat=matmat_XH_X, 

1851 shape=(min(A.shape), min(A.shape))) 

1852 

1853 # Get a low rank approximation of the implicitly defined gramian matrix. 

1854 # This is not a stable way to approach the problem. 

1855 if solver == 'lobpcg': 

1856 

1857 if k == 1 and v0 is not None: 

1858 X = np.reshape(v0, (-1, 1)) 

1859 else: 

1860 X = np.random.RandomState(52).randn(min(A.shape), k) 

1861 

1862 eigvals, eigvec = lobpcg(XH_X, X, tol=tol ** 2, maxiter=maxiter, 

1863 largest=largest) 

1864 

1865 elif solver == 'arpack' or solver is None: 

1866 eigvals, eigvec = eigsh(XH_X, k=k, tol=tol ** 2, maxiter=maxiter, 

1867 ncv=ncv, which=which, v0=v0) 

1868 

1869 else: 

1870 raise ValueError("solver must be either 'arpack', or 'lobpcg'.") 

1871 

1872 # Gramian matrices have real non-negative eigenvalues. 

1873 eigvals = np.maximum(eigvals.real, 0) 

1874 

1875 # Use the sophisticated detection of small eigenvalues from pinvh. 

1876 t = eigvec.dtype.char.lower() 

1877 factor = {'f': 1E3, 'd': 1E6} 

1878 cond = factor[t] * np.finfo(t).eps 

1879 cutoff = cond * np.max(eigvals) 

1880 

1881 # Get a mask indicating which eigenpairs are not degenerately tiny, 

1882 # and create the re-ordered array of thresholded singular values. 

1883 above_cutoff = (eigvals > cutoff) 

1884 nlarge = above_cutoff.sum() 

1885 nsmall = k - nlarge 

1886 slarge = np.sqrt(eigvals[above_cutoff]) 

1887 s = np.zeros_like(eigvals) 

1888 s[:nlarge] = slarge 

1889 if not return_singular_vectors: 

1890 return np.sort(s) 

1891 

1892 if n > m: 

1893 vlarge = eigvec[:, above_cutoff] 

1894 ularge = X_matmat(vlarge) / slarge if return_singular_vectors != 'vh' else None 

1895 vhlarge = _herm(vlarge) 

1896 else: 

1897 ularge = eigvec[:, above_cutoff] 

1898 vhlarge = _herm(X_matmat(ularge) / slarge) if return_singular_vectors != 'u' else None 

1899 

1900 u = _augmented_orthonormal_cols(ularge, nsmall) if ularge is not None else None 

1901 vh = _augmented_orthonormal_rows(vhlarge, nsmall) if vhlarge is not None else None 

1902 

1903 indexes_sorted = np.argsort(s) 

1904 s = s[indexes_sorted] 

1905 if u is not None: 

1906 u = u[:, indexes_sorted] 

1907 if vh is not None: 

1908 vh = vh[indexes_sorted] 

1909 

1910 return u, s, vh