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# 

2# Author: Pearu Peterson, March 2002 

3# 

4# w/ additions by Travis Oliphant, March 2002 

5# and Jake Vanderplas, August 2012 

6 

7from warnings import warn 

8import numpy as np 

9from numpy import atleast_1d, atleast_2d 

10from .flinalg import get_flinalg_funcs 

11from .lapack import get_lapack_funcs, _compute_lwork 

12from .misc import LinAlgError, _datacopied, LinAlgWarning 

13from .decomp import _asarray_validated 

14from . import decomp, decomp_svd 

15from ._solve_toeplitz import levinson 

16 

17__all__ = ['solve', 'solve_triangular', 'solveh_banded', 'solve_banded', 

18 'solve_toeplitz', 'solve_circulant', 'inv', 'det', 'lstsq', 

19 'pinv', 'pinv2', 'pinvh', 'matrix_balance'] 

20 

21 

22# Linear equations 

23def _solve_check(n, info, lamch=None, rcond=None): 

24 """ Check arguments during the different steps of the solution phase """ 

25 if info < 0: 

26 raise ValueError('LAPACK reported an illegal value in {}-th argument' 

27 '.'.format(-info)) 

28 elif 0 < info: 

29 raise LinAlgError('Matrix is singular.') 

30 

31 if lamch is None: 

32 return 

33 E = lamch('E') 

34 if rcond < E: 

35 warn('Ill-conditioned matrix (rcond={:.6g}): ' 

36 'result may not be accurate.'.format(rcond), 

37 LinAlgWarning, stacklevel=3) 

38 

39 

40def solve(a, b, sym_pos=False, lower=False, overwrite_a=False, 

41 overwrite_b=False, debug=None, check_finite=True, assume_a='gen', 

42 transposed=False): 

43 """ 

44 Solves the linear equation set ``a * x = b`` for the unknown ``x`` 

45 for square ``a`` matrix. 

46 

47 If the data matrix is known to be a particular type then supplying the 

48 corresponding string to ``assume_a`` key chooses the dedicated solver. 

49 The available options are 

50 

51 =================== ======== 

52 generic matrix 'gen' 

53 symmetric 'sym' 

54 hermitian 'her' 

55 positive definite 'pos' 

56 =================== ======== 

57 

58 If omitted, ``'gen'`` is the default structure. 

59 

60 The datatype of the arrays define which solver is called regardless 

61 of the values. In other words, even when the complex array entries have 

62 precisely zero imaginary parts, the complex solver will be called based 

63 on the data type of the array. 

64 

65 Parameters 

66 ---------- 

67 a : (N, N) array_like 

68 Square input data 

69 b : (N, NRHS) array_like 

70 Input data for the right hand side. 

71 sym_pos : bool, optional 

72 Assume `a` is symmetric and positive definite. This key is deprecated 

73 and assume_a = 'pos' keyword is recommended instead. The functionality 

74 is the same. It will be removed in the future. 

75 lower : bool, optional 

76 If True, only the data contained in the lower triangle of `a`. Default 

77 is to use upper triangle. (ignored for ``'gen'``) 

78 overwrite_a : bool, optional 

79 Allow overwriting data in `a` (may enhance performance). 

80 Default is False. 

81 overwrite_b : bool, optional 

82 Allow overwriting data in `b` (may enhance performance). 

83 Default is False. 

84 check_finite : bool, optional 

85 Whether to check that the input matrices contain only finite numbers. 

86 Disabling may give a performance gain, but may result in problems 

87 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

88 assume_a : str, optional 

89 Valid entries are explained above. 

90 transposed: bool, optional 

91 If True, ``a^T x = b`` for real matrices, raises `NotImplementedError` 

92 for complex matrices (only for True). 

93 

94 Returns 

95 ------- 

96 x : (N, NRHS) ndarray 

97 The solution array. 

98 

99 Raises 

100 ------ 

101 ValueError 

102 If size mismatches detected or input a is not square. 

103 LinAlgError 

104 If the matrix is singular. 

105 LinAlgWarning 

106 If an ill-conditioned input a is detected. 

107 NotImplementedError 

108 If transposed is True and input a is a complex matrix. 

109 

110 Examples 

111 -------- 

112 Given `a` and `b`, solve for `x`: 

113 

114 >>> a = np.array([[3, 2, 0], [1, -1, 0], [0, 5, 1]]) 

115 >>> b = np.array([2, 4, -1]) 

116 >>> from scipy import linalg 

117 >>> x = linalg.solve(a, b) 

118 >>> x 

119 array([ 2., -2., 9.]) 

120 >>> np.dot(a, x) == b 

121 array([ True, True, True], dtype=bool) 

122 

123 Notes 

124 ----- 

125 If the input b matrix is a 1-D array with N elements, when supplied 

126 together with an NxN input a, it is assumed as a valid column vector 

127 despite the apparent size mismatch. This is compatible with the 

128 numpy.dot() behavior and the returned result is still 1-D array. 

129 

130 The generic, symmetric, hermitian and positive definite solutions are 

131 obtained via calling ?GESV, ?SYSV, ?HESV, and ?POSV routines of 

132 LAPACK respectively. 

133 """ 

134 # Flags for 1-D or N-D right-hand side 

135 b_is_1D = False 

136 

137 a1 = atleast_2d(_asarray_validated(a, check_finite=check_finite)) 

138 b1 = atleast_1d(_asarray_validated(b, check_finite=check_finite)) 

139 n = a1.shape[0] 

140 

141 overwrite_a = overwrite_a or _datacopied(a1, a) 

142 overwrite_b = overwrite_b or _datacopied(b1, b) 

143 

144 if a1.shape[0] != a1.shape[1]: 

145 raise ValueError('Input a needs to be a square matrix.') 

146 

147 if n != b1.shape[0]: 

148 # Last chance to catch 1x1 scalar a and 1-D b arrays 

149 if not (n == 1 and b1.size != 0): 

150 raise ValueError('Input b has to have same number of rows as ' 

151 'input a') 

152 

153 # accommodate empty arrays 

154 if b1.size == 0: 

155 return np.asfortranarray(b1.copy()) 

156 

157 # regularize 1-D b arrays to 2D 

158 if b1.ndim == 1: 

159 if n == 1: 

160 b1 = b1[None, :] 

161 else: 

162 b1 = b1[:, None] 

163 b_is_1D = True 

164 

165 # Backwards compatibility - old keyword. 

166 if sym_pos: 

167 assume_a = 'pos' 

168 

169 if assume_a not in ('gen', 'sym', 'her', 'pos'): 

170 raise ValueError('{} is not a recognized matrix structure' 

171 ''.format(assume_a)) 

172 

173 # Deprecate keyword "debug" 

174 if debug is not None: 

175 warn('Use of the "debug" keyword is deprecated ' 

176 'and this keyword will be removed in future ' 

177 'versions of SciPy.', DeprecationWarning, stacklevel=2) 

178 

179 # Get the correct lamch function. 

180 # The LAMCH functions only exists for S and D 

181 # So for complex values we have to convert to real/double. 

182 if a1.dtype.char in 'fF': # single precision 

183 lamch = get_lapack_funcs('lamch', dtype='f') 

184 else: 

185 lamch = get_lapack_funcs('lamch', dtype='d') 

186 

187 # Currently we do not have the other forms of the norm calculators 

188 # lansy, lanpo, lanhe. 

189 # However, in any case they only reduce computations slightly... 

190 lange = get_lapack_funcs('lange', (a1,)) 

191 

192 # Since the I-norm and 1-norm are the same for symmetric matrices 

193 # we can collect them all in this one call 

194 # Note however, that when issuing 'gen' and form!='none', then 

195 # the I-norm should be used 

196 if transposed: 

197 trans = 1 

198 norm = 'I' 

199 if np.iscomplexobj(a1): 

200 raise NotImplementedError('scipy.linalg.solve can currently ' 

201 'not solve a^T x = b or a^H x = b ' 

202 'for complex matrices.') 

203 else: 

204 trans = 0 

205 norm = '1' 

206 

207 anorm = lange(norm, a1) 

208 

209 # Generalized case 'gesv' 

210 if assume_a == 'gen': 

211 gecon, getrf, getrs = get_lapack_funcs(('gecon', 'getrf', 'getrs'), 

212 (a1, b1)) 

213 lu, ipvt, info = getrf(a1, overwrite_a=overwrite_a) 

214 _solve_check(n, info) 

215 x, info = getrs(lu, ipvt, b1, 

216 trans=trans, overwrite_b=overwrite_b) 

217 _solve_check(n, info) 

218 rcond, info = gecon(lu, anorm, norm=norm) 

219 # Hermitian case 'hesv' 

220 elif assume_a == 'her': 

221 hecon, hesv, hesv_lw = get_lapack_funcs(('hecon', 'hesv', 

222 'hesv_lwork'), (a1, b1)) 

223 lwork = _compute_lwork(hesv_lw, n, lower) 

224 lu, ipvt, x, info = hesv(a1, b1, lwork=lwork, 

225 lower=lower, 

226 overwrite_a=overwrite_a, 

227 overwrite_b=overwrite_b) 

228 _solve_check(n, info) 

229 rcond, info = hecon(lu, ipvt, anorm) 

230 # Symmetric case 'sysv' 

231 elif assume_a == 'sym': 

232 sycon, sysv, sysv_lw = get_lapack_funcs(('sycon', 'sysv', 

233 'sysv_lwork'), (a1, b1)) 

234 lwork = _compute_lwork(sysv_lw, n, lower) 

235 lu, ipvt, x, info = sysv(a1, b1, lwork=lwork, 

236 lower=lower, 

237 overwrite_a=overwrite_a, 

238 overwrite_b=overwrite_b) 

239 _solve_check(n, info) 

240 rcond, info = sycon(lu, ipvt, anorm) 

241 # Positive definite case 'posv' 

242 else: 

243 pocon, posv = get_lapack_funcs(('pocon', 'posv'), 

244 (a1, b1)) 

245 lu, x, info = posv(a1, b1, lower=lower, 

246 overwrite_a=overwrite_a, 

247 overwrite_b=overwrite_b) 

248 _solve_check(n, info) 

249 rcond, info = pocon(lu, anorm) 

250 

251 _solve_check(n, info, lamch, rcond) 

252 

253 if b_is_1D: 

254 x = x.ravel() 

255 

256 return x 

257 

258 

259def solve_triangular(a, b, trans=0, lower=False, unit_diagonal=False, 

260 overwrite_b=False, debug=None, check_finite=True): 

261 """ 

262 Solve the equation `a x = b` for `x`, assuming a is a triangular matrix. 

263 

264 Parameters 

265 ---------- 

266 a : (M, M) array_like 

267 A triangular matrix 

268 b : (M,) or (M, N) array_like 

269 Right-hand side matrix in `a x = b` 

270 lower : bool, optional 

271 Use only data contained in the lower triangle of `a`. 

272 Default is to use upper triangle. 

273 trans : {0, 1, 2, 'N', 'T', 'C'}, optional 

274 Type of system to solve: 

275 

276 ======== ========= 

277 trans system 

278 ======== ========= 

279 0 or 'N' a x = b 

280 1 or 'T' a^T x = b 

281 2 or 'C' a^H x = b 

282 ======== ========= 

283 unit_diagonal : bool, optional 

284 If True, diagonal elements of `a` are assumed to be 1 and 

285 will not be referenced. 

286 overwrite_b : bool, optional 

287 Allow overwriting data in `b` (may enhance performance) 

288 check_finite : bool, optional 

289 Whether to check that the input matrices contain only finite numbers. 

290 Disabling may give a performance gain, but may result in problems 

291 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

292 

293 Returns 

294 ------- 

295 x : (M,) or (M, N) ndarray 

296 Solution to the system `a x = b`. Shape of return matches `b`. 

297 

298 Raises 

299 ------ 

300 LinAlgError 

301 If `a` is singular 

302 

303 Notes 

304 ----- 

305 .. versionadded:: 0.9.0 

306 

307 Examples 

308 -------- 

309 Solve the lower triangular system a x = b, where:: 

310 

311 [3 0 0 0] [4] 

312 a = [2 1 0 0] b = [2] 

313 [1 0 1 0] [4] 

314 [1 1 1 1] [2] 

315 

316 >>> from scipy.linalg import solve_triangular 

317 >>> a = np.array([[3, 0, 0, 0], [2, 1, 0, 0], [1, 0, 1, 0], [1, 1, 1, 1]]) 

318 >>> b = np.array([4, 2, 4, 2]) 

319 >>> x = solve_triangular(a, b, lower=True) 

320 >>> x 

321 array([ 1.33333333, -0.66666667, 2.66666667, -1.33333333]) 

322 >>> a.dot(x) # Check the result 

323 array([ 4., 2., 4., 2.]) 

324 

325 """ 

326 

327 # Deprecate keyword "debug" 

328 if debug is not None: 

329 warn('Use of the "debug" keyword is deprecated ' 

330 'and this keyword will be removed in the future ' 

331 'versions of SciPy.', DeprecationWarning, stacklevel=2) 

332 

333 a1 = _asarray_validated(a, check_finite=check_finite) 

334 b1 = _asarray_validated(b, check_finite=check_finite) 

335 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]: 

336 raise ValueError('expected square matrix') 

337 if a1.shape[0] != b1.shape[0]: 

338 raise ValueError('incompatible dimensions') 

339 overwrite_b = overwrite_b or _datacopied(b1, b) 

340 if debug: 

341 print('solve:overwrite_b=', overwrite_b) 

342 trans = {'N': 0, 'T': 1, 'C': 2}.get(trans, trans) 

343 trtrs, = get_lapack_funcs(('trtrs',), (a1, b1)) 

344 if a1.flags.f_contiguous or trans == 2: 

345 x, info = trtrs(a1, b1, overwrite_b=overwrite_b, lower=lower, 

346 trans=trans, unitdiag=unit_diagonal) 

347 else: 

348 # transposed system is solved since trtrs expects Fortran ordering 

349 x, info = trtrs(a1.T, b1, overwrite_b=overwrite_b, lower=not lower, 

350 trans=not trans, unitdiag=unit_diagonal) 

351 

352 if info == 0: 

353 return x 

354 if info > 0: 

355 raise LinAlgError("singular matrix: resolution failed at diagonal %d" % 

356 (info-1)) 

357 raise ValueError('illegal value in %dth argument of internal trtrs' % 

358 (-info)) 

359 

360 

361def solve_banded(l_and_u, ab, b, overwrite_ab=False, overwrite_b=False, 

362 debug=None, check_finite=True): 

363 """ 

364 Solve the equation a x = b for x, assuming a is banded matrix. 

365 

366 The matrix a is stored in `ab` using the matrix diagonal ordered form:: 

367 

368 ab[u + i - j, j] == a[i,j] 

369 

370 Example of `ab` (shape of a is (6,6), `u` =1, `l` =2):: 

371 

372 * a01 a12 a23 a34 a45 

373 a00 a11 a22 a33 a44 a55 

374 a10 a21 a32 a43 a54 * 

375 a20 a31 a42 a53 * * 

376 

377 Parameters 

378 ---------- 

379 (l, u) : (integer, integer) 

380 Number of non-zero lower and upper diagonals 

381 ab : (`l` + `u` + 1, M) array_like 

382 Banded matrix 

383 b : (M,) or (M, K) array_like 

384 Right-hand side 

385 overwrite_ab : bool, optional 

386 Discard data in `ab` (may enhance performance) 

387 overwrite_b : bool, optional 

388 Discard data in `b` (may enhance performance) 

389 check_finite : bool, optional 

390 Whether to check that the input matrices contain only finite numbers. 

391 Disabling may give a performance gain, but may result in problems 

392 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

393 

394 Returns 

395 ------- 

396 x : (M,) or (M, K) ndarray 

397 The solution to the system a x = b. Returned shape depends on the 

398 shape of `b`. 

399 

400 Examples 

401 -------- 

402 Solve the banded system a x = b, where:: 

403 

404 [5 2 -1 0 0] [0] 

405 [1 4 2 -1 0] [1] 

406 a = [0 1 3 2 -1] b = [2] 

407 [0 0 1 2 2] [2] 

408 [0 0 0 1 1] [3] 

409 

410 There is one nonzero diagonal below the main diagonal (l = 1), and 

411 two above (u = 2). The diagonal banded form of the matrix is:: 

412 

413 [* * -1 -1 -1] 

414 ab = [* 2 2 2 2] 

415 [5 4 3 2 1] 

416 [1 1 1 1 *] 

417 

418 >>> from scipy.linalg import solve_banded 

419 >>> ab = np.array([[0, 0, -1, -1, -1], 

420 ... [0, 2, 2, 2, 2], 

421 ... [5, 4, 3, 2, 1], 

422 ... [1, 1, 1, 1, 0]]) 

423 >>> b = np.array([0, 1, 2, 2, 3]) 

424 >>> x = solve_banded((1, 2), ab, b) 

425 >>> x 

426 array([-2.37288136, 3.93220339, -4. , 4.3559322 , -1.3559322 ]) 

427 

428 """ 

429 

430 # Deprecate keyword "debug" 

431 if debug is not None: 

432 warn('Use of the "debug" keyword is deprecated ' 

433 'and this keyword will be removed in the future ' 

434 'versions of SciPy.', DeprecationWarning, stacklevel=2) 

435 

436 a1 = _asarray_validated(ab, check_finite=check_finite, as_inexact=True) 

437 b1 = _asarray_validated(b, check_finite=check_finite, as_inexact=True) 

438 # Validate shapes. 

439 if a1.shape[-1] != b1.shape[0]: 

440 raise ValueError("shapes of ab and b are not compatible.") 

441 (nlower, nupper) = l_and_u 

442 if nlower + nupper + 1 != a1.shape[0]: 

443 raise ValueError("invalid values for the number of lower and upper " 

444 "diagonals: l+u+1 (%d) does not equal ab.shape[0] " 

445 "(%d)" % (nlower + nupper + 1, ab.shape[0])) 

446 

447 overwrite_b = overwrite_b or _datacopied(b1, b) 

448 if a1.shape[-1] == 1: 

449 b2 = np.array(b1, copy=(not overwrite_b)) 

450 b2 /= a1[1, 0] 

451 return b2 

452 if nlower == nupper == 1: 

453 overwrite_ab = overwrite_ab or _datacopied(a1, ab) 

454 gtsv, = get_lapack_funcs(('gtsv',), (a1, b1)) 

455 du = a1[0, 1:] 

456 d = a1[1, :] 

457 dl = a1[2, :-1] 

458 du2, d, du, x, info = gtsv(dl, d, du, b1, overwrite_ab, overwrite_ab, 

459 overwrite_ab, overwrite_b) 

460 else: 

461 gbsv, = get_lapack_funcs(('gbsv',), (a1, b1)) 

462 a2 = np.zeros((2*nlower + nupper + 1, a1.shape[1]), dtype=gbsv.dtype) 

463 a2[nlower:, :] = a1 

464 lu, piv, x, info = gbsv(nlower, nupper, a2, b1, overwrite_ab=True, 

465 overwrite_b=overwrite_b) 

466 if info == 0: 

467 return x 

468 if info > 0: 

469 raise LinAlgError("singular matrix") 

470 raise ValueError('illegal value in %d-th argument of internal ' 

471 'gbsv/gtsv' % -info) 

472 

473 

474def solveh_banded(ab, b, overwrite_ab=False, overwrite_b=False, lower=False, 

475 check_finite=True): 

476 """ 

477 Solve equation a x = b. a is Hermitian positive-definite banded matrix. 

478 

479 The matrix a is stored in `ab` either in lower diagonal or upper 

480 diagonal ordered form: 

481 

482 ab[u + i - j, j] == a[i,j] (if upper form; i <= j) 

483 ab[ i - j, j] == a[i,j] (if lower form; i >= j) 

484 

485 Example of `ab` (shape of a is (6, 6), `u` =2):: 

486 

487 upper form: 

488 * * a02 a13 a24 a35 

489 * a01 a12 a23 a34 a45 

490 a00 a11 a22 a33 a44 a55 

491 

492 lower form: 

493 a00 a11 a22 a33 a44 a55 

494 a10 a21 a32 a43 a54 * 

495 a20 a31 a42 a53 * * 

496 

497 Cells marked with * are not used. 

498 

499 Parameters 

500 ---------- 

501 ab : (`u` + 1, M) array_like 

502 Banded matrix 

503 b : (M,) or (M, K) array_like 

504 Right-hand side 

505 overwrite_ab : bool, optional 

506 Discard data in `ab` (may enhance performance) 

507 overwrite_b : bool, optional 

508 Discard data in `b` (may enhance performance) 

509 lower : bool, optional 

510 Is the matrix in the lower form. (Default is upper form) 

511 check_finite : bool, optional 

512 Whether to check that the input matrices contain only finite numbers. 

513 Disabling may give a performance gain, but may result in problems 

514 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

515 

516 Returns 

517 ------- 

518 x : (M,) or (M, K) ndarray 

519 The solution to the system a x = b. Shape of return matches shape 

520 of `b`. 

521 

522 Examples 

523 -------- 

524 Solve the banded system A x = b, where:: 

525 

526 [ 4 2 -1 0 0 0] [1] 

527 [ 2 5 2 -1 0 0] [2] 

528 A = [-1 2 6 2 -1 0] b = [2] 

529 [ 0 -1 2 7 2 -1] [3] 

530 [ 0 0 -1 2 8 2] [3] 

531 [ 0 0 0 -1 2 9] [3] 

532 

533 >>> from scipy.linalg import solveh_banded 

534 

535 `ab` contains the main diagonal and the nonzero diagonals below the 

536 main diagonal. That is, we use the lower form: 

537 

538 >>> ab = np.array([[ 4, 5, 6, 7, 8, 9], 

539 ... [ 2, 2, 2, 2, 2, 0], 

540 ... [-1, -1, -1, -1, 0, 0]]) 

541 >>> b = np.array([1, 2, 2, 3, 3, 3]) 

542 >>> x = solveh_banded(ab, b, lower=True) 

543 >>> x 

544 array([ 0.03431373, 0.45938375, 0.05602241, 0.47759104, 0.17577031, 

545 0.34733894]) 

546 

547 

548 Solve the Hermitian banded system H x = b, where:: 

549 

550 [ 8 2-1j 0 0 ] [ 1 ] 

551 H = [2+1j 5 1j 0 ] b = [1+1j] 

552 [ 0 -1j 9 -2-1j] [1-2j] 

553 [ 0 0 -2+1j 6 ] [ 0 ] 

554 

555 In this example, we put the upper diagonals in the array `hb`: 

556 

557 >>> hb = np.array([[0, 2-1j, 1j, -2-1j], 

558 ... [8, 5, 9, 6 ]]) 

559 >>> b = np.array([1, 1+1j, 1-2j, 0]) 

560 >>> x = solveh_banded(hb, b) 

561 >>> x 

562 array([ 0.07318536-0.02939412j, 0.11877624+0.17696461j, 

563 0.10077984-0.23035393j, -0.00479904-0.09358128j]) 

564 

565 """ 

566 a1 = _asarray_validated(ab, check_finite=check_finite) 

567 b1 = _asarray_validated(b, check_finite=check_finite) 

568 # Validate shapes. 

569 if a1.shape[-1] != b1.shape[0]: 

570 raise ValueError("shapes of ab and b are not compatible.") 

571 

572 overwrite_b = overwrite_b or _datacopied(b1, b) 

573 overwrite_ab = overwrite_ab or _datacopied(a1, ab) 

574 

575 if a1.shape[0] == 2: 

576 ptsv, = get_lapack_funcs(('ptsv',), (a1, b1)) 

577 if lower: 

578 d = a1[0, :].real 

579 e = a1[1, :-1] 

580 else: 

581 d = a1[1, :].real 

582 e = a1[0, 1:].conj() 

583 d, du, x, info = ptsv(d, e, b1, overwrite_ab, overwrite_ab, 

584 overwrite_b) 

585 else: 

586 pbsv, = get_lapack_funcs(('pbsv',), (a1, b1)) 

587 c, x, info = pbsv(a1, b1, lower=lower, overwrite_ab=overwrite_ab, 

588 overwrite_b=overwrite_b) 

589 if info > 0: 

590 raise LinAlgError("%dth leading minor not positive definite" % info) 

591 if info < 0: 

592 raise ValueError('illegal value in %dth argument of internal ' 

593 'pbsv' % -info) 

594 return x 

595 

596 

597def solve_toeplitz(c_or_cr, b, check_finite=True): 

598 """Solve a Toeplitz system using Levinson Recursion 

599 

600 The Toeplitz matrix has constant diagonals, with c as its first column 

601 and r as its first row. If r is not given, ``r == conjugate(c)`` is 

602 assumed. 

603 

604 Parameters 

605 ---------- 

606 c_or_cr : array_like or tuple of (array_like, array_like) 

607 The vector ``c``, or a tuple of arrays (``c``, ``r``). Whatever the 

608 actual shape of ``c``, it will be converted to a 1-D array. If not 

609 supplied, ``r = conjugate(c)`` is assumed; in this case, if c[0] is 

610 real, the Toeplitz matrix is Hermitian. r[0] is ignored; the first row 

611 of the Toeplitz matrix is ``[c[0], r[1:]]``. Whatever the actual shape 

612 of ``r``, it will be converted to a 1-D array. 

613 b : (M,) or (M, K) array_like 

614 Right-hand side in ``T x = b``. 

615 check_finite : bool, optional 

616 Whether to check that the input matrices contain only finite numbers. 

617 Disabling may give a performance gain, but may result in problems 

618 (result entirely NaNs) if the inputs do contain infinities or NaNs. 

619 

620 Returns 

621 ------- 

622 x : (M,) or (M, K) ndarray 

623 The solution to the system ``T x = b``. Shape of return matches shape 

624 of `b`. 

625 

626 See Also 

627 -------- 

628 toeplitz : Toeplitz matrix 

629 

630 Notes 

631 ----- 

632 The solution is computed using Levinson-Durbin recursion, which is faster 

633 than generic least-squares methods, but can be less numerically stable. 

634 

635 Examples 

636 -------- 

637 Solve the Toeplitz system T x = b, where:: 

638 

639 [ 1 -1 -2 -3] [1] 

640 T = [ 3 1 -1 -2] b = [2] 

641 [ 6 3 1 -1] [2] 

642 [10 6 3 1] [5] 

643 

644 To specify the Toeplitz matrix, only the first column and the first 

645 row are needed. 

646 

647 >>> c = np.array([1, 3, 6, 10]) # First column of T 

648 >>> r = np.array([1, -1, -2, -3]) # First row of T 

649 >>> b = np.array([1, 2, 2, 5]) 

650 

651 >>> from scipy.linalg import solve_toeplitz, toeplitz 

652 >>> x = solve_toeplitz((c, r), b) 

653 >>> x 

654 array([ 1.66666667, -1. , -2.66666667, 2.33333333]) 

655 

656 Check the result by creating the full Toeplitz matrix and 

657 multiplying it by `x`. We should get `b`. 

658 

659 >>> T = toeplitz(c, r) 

660 >>> T.dot(x) 

661 array([ 1., 2., 2., 5.]) 

662 

663 """ 

664 # If numerical stability of this algorithm is a problem, a future 

665 # developer might consider implementing other O(N^2) Toeplitz solvers, 

666 # such as GKO (https://www.jstor.org/stable/2153371) or Bareiss. 

667 if isinstance(c_or_cr, tuple): 

668 c, r = c_or_cr 

669 c = _asarray_validated(c, check_finite=check_finite).ravel() 

670 r = _asarray_validated(r, check_finite=check_finite).ravel() 

671 else: 

672 c = _asarray_validated(c_or_cr, check_finite=check_finite).ravel() 

673 r = c.conjugate() 

674 

675 # Form a 1-D array of values to be used in the matrix, containing a reversed 

676 # copy of r[1:], followed by c. 

677 vals = np.concatenate((r[-1:0:-1], c)) 

678 if b is None: 

679 raise ValueError('illegal value, `b` is a required argument') 

680 

681 b = _asarray_validated(b) 

682 if vals.shape[0] != (2*b.shape[0] - 1): 

683 raise ValueError('incompatible dimensions') 

684 if np.iscomplexobj(vals) or np.iscomplexobj(b): 

685 vals = np.asarray(vals, dtype=np.complex128, order='c') 

686 b = np.asarray(b, dtype=np.complex128) 

687 else: 

688 vals = np.asarray(vals, dtype=np.double, order='c') 

689 b = np.asarray(b, dtype=np.double) 

690 

691 if b.ndim == 1: 

692 x, _ = levinson(vals, np.ascontiguousarray(b)) 

693 else: 

694 b_shape = b.shape 

695 b = b.reshape(b.shape[0], -1) 

696 x = np.column_stack([levinson(vals, np.ascontiguousarray(b[:, i]))[0] 

697 for i in range(b.shape[1])]) 

698 x = x.reshape(*b_shape) 

699 

700 return x 

701 

702 

703def _get_axis_len(aname, a, axis): 

704 ax = axis 

705 if ax < 0: 

706 ax += a.ndim 

707 if 0 <= ax < a.ndim: 

708 return a.shape[ax] 

709 raise ValueError("'%saxis' entry is out of bounds" % (aname,)) 

710 

711 

712def solve_circulant(c, b, singular='raise', tol=None, 

713 caxis=-1, baxis=0, outaxis=0): 

714 """Solve C x = b for x, where C is a circulant matrix. 

715 

716 `C` is the circulant matrix associated with the vector `c`. 

717 

718 The system is solved by doing division in Fourier space. The 

719 calculation is:: 

720 

721 x = ifft(fft(b) / fft(c)) 

722 

723 where `fft` and `ifft` are the fast Fourier transform and its inverse, 

724 respectively. For a large vector `c`, this is *much* faster than 

725 solving the system with the full circulant matrix. 

726 

727 Parameters 

728 ---------- 

729 c : array_like 

730 The coefficients of the circulant matrix. 

731 b : array_like 

732 Right-hand side matrix in ``a x = b``. 

733 singular : str, optional 

734 This argument controls how a near singular circulant matrix is 

735 handled. If `singular` is "raise" and the circulant matrix is 

736 near singular, a `LinAlgError` is raised. If `singular` is 

737 "lstsq", the least squares solution is returned. Default is "raise". 

738 tol : float, optional 

739 If any eigenvalue of the circulant matrix has an absolute value 

740 that is less than or equal to `tol`, the matrix is considered to be 

741 near singular. If not given, `tol` is set to:: 

742 

743 tol = abs_eigs.max() * abs_eigs.size * np.finfo(np.float64).eps 

744 

745 where `abs_eigs` is the array of absolute values of the eigenvalues 

746 of the circulant matrix. 

747 caxis : int 

748 When `c` has dimension greater than 1, it is viewed as a collection 

749 of circulant vectors. In this case, `caxis` is the axis of `c` that 

750 holds the vectors of circulant coefficients. 

751 baxis : int 

752 When `b` has dimension greater than 1, it is viewed as a collection 

753 of vectors. In this case, `baxis` is the axis of `b` that holds the 

754 right-hand side vectors. 

755 outaxis : int 

756 When `c` or `b` are multidimensional, the value returned by 

757 `solve_circulant` is multidimensional. In this case, `outaxis` is 

758 the axis of the result that holds the solution vectors. 

759 

760 Returns 

761 ------- 

762 x : ndarray 

763 Solution to the system ``C x = b``. 

764 

765 Raises 

766 ------ 

767 LinAlgError 

768 If the circulant matrix associated with `c` is near singular. 

769 

770 See Also 

771 -------- 

772 circulant : circulant matrix 

773 

774 Notes 

775 ----- 

776 For a 1-D vector `c` with length `m`, and an array `b` 

777 with shape ``(m, ...)``, 

778 

779 solve_circulant(c, b) 

780 

781 returns the same result as 

782 

783 solve(circulant(c), b) 

784 

785 where `solve` and `circulant` are from `scipy.linalg`. 

786 

787 .. versionadded:: 0.16.0 

788 

789 Examples 

790 -------- 

791 >>> from scipy.linalg import solve_circulant, solve, circulant, lstsq 

792 

793 >>> c = np.array([2, 2, 4]) 

794 >>> b = np.array([1, 2, 3]) 

795 >>> solve_circulant(c, b) 

796 array([ 0.75, -0.25, 0.25]) 

797 

798 Compare that result to solving the system with `scipy.linalg.solve`: 

799 

800 >>> solve(circulant(c), b) 

801 array([ 0.75, -0.25, 0.25]) 

802 

803 A singular example: 

804 

805 >>> c = np.array([1, 1, 0, 0]) 

806 >>> b = np.array([1, 2, 3, 4]) 

807 

808 Calling ``solve_circulant(c, b)`` will raise a `LinAlgError`. For the 

809 least square solution, use the option ``singular='lstsq'``: 

810 

811 >>> solve_circulant(c, b, singular='lstsq') 

812 array([ 0.25, 1.25, 2.25, 1.25]) 

813 

814 Compare to `scipy.linalg.lstsq`: 

815 

816 >>> x, resid, rnk, s = lstsq(circulant(c), b) 

817 >>> x 

818 array([ 0.25, 1.25, 2.25, 1.25]) 

819 

820 A broadcasting example: 

821 

822 Suppose we have the vectors of two circulant matrices stored in an array 

823 with shape (2, 5), and three `b` vectors stored in an array with shape 

824 (3, 5). For example, 

825 

826 >>> c = np.array([[1.5, 2, 3, 0, 0], [1, 1, 4, 3, 2]]) 

827 >>> b = np.arange(15).reshape(-1, 5) 

828 

829 We want to solve all combinations of circulant matrices and `b` vectors, 

830 with the result stored in an array with shape (2, 3, 5). When we 

831 disregard the axes of `c` and `b` that hold the vectors of coefficients, 

832 the shapes of the collections are (2,) and (3,), respectively, which are 

833 not compatible for broadcasting. To have a broadcast result with shape 

834 (2, 3), we add a trivial dimension to `c`: ``c[:, np.newaxis, :]`` has 

835 shape (2, 1, 5). The last dimension holds the coefficients of the 

836 circulant matrices, so when we call `solve_circulant`, we can use the 

837 default ``caxis=-1``. The coefficients of the `b` vectors are in the last 

838 dimension of the array `b`, so we use ``baxis=-1``. If we use the 

839 default `outaxis`, the result will have shape (5, 2, 3), so we'll use 

840 ``outaxis=-1`` to put the solution vectors in the last dimension. 

841 

842 >>> x = solve_circulant(c[:, np.newaxis, :], b, baxis=-1, outaxis=-1) 

843 >>> x.shape 

844 (2, 3, 5) 

845 >>> np.set_printoptions(precision=3) # For compact output of numbers. 

846 >>> x 

847 array([[[-0.118, 0.22 , 1.277, -0.142, 0.302], 

848 [ 0.651, 0.989, 2.046, 0.627, 1.072], 

849 [ 1.42 , 1.758, 2.816, 1.396, 1.841]], 

850 [[ 0.401, 0.304, 0.694, -0.867, 0.377], 

851 [ 0.856, 0.758, 1.149, -0.412, 0.831], 

852 [ 1.31 , 1.213, 1.603, 0.042, 1.286]]]) 

853 

854 Check by solving one pair of `c` and `b` vectors (cf. ``x[1, 1, :]``): 

855 

856 >>> solve_circulant(c[1], b[1, :]) 

857 array([ 0.856, 0.758, 1.149, -0.412, 0.831]) 

858 

859 """ 

860 c = np.atleast_1d(c) 

861 nc = _get_axis_len("c", c, caxis) 

862 b = np.atleast_1d(b) 

863 nb = _get_axis_len("b", b, baxis) 

864 if nc != nb: 

865 raise ValueError('Incompatible c and b axis lengths') 

866 

867 fc = np.fft.fft(np.rollaxis(c, caxis, c.ndim), axis=-1) 

868 abs_fc = np.abs(fc) 

869 if tol is None: 

870 # This is the same tolerance as used in np.linalg.matrix_rank. 

871 tol = abs_fc.max(axis=-1) * nc * np.finfo(np.float64).eps 

872 if tol.shape != (): 

873 tol.shape = tol.shape + (1,) 

874 else: 

875 tol = np.atleast_1d(tol) 

876 

877 near_zeros = abs_fc <= tol 

878 is_near_singular = np.any(near_zeros) 

879 if is_near_singular: 

880 if singular == 'raise': 

881 raise LinAlgError("near singular circulant matrix.") 

882 else: 

883 # Replace the small values with 1 to avoid errors in the 

884 # division fb/fc below. 

885 fc[near_zeros] = 1 

886 

887 fb = np.fft.fft(np.rollaxis(b, baxis, b.ndim), axis=-1) 

888 

889 q = fb / fc 

890 

891 if is_near_singular: 

892 # `near_zeros` is a boolean array, same shape as `c`, that is 

893 # True where `fc` is (near) zero. `q` is the broadcasted result 

894 # of fb / fc, so to set the values of `q` to 0 where `fc` is near 

895 # zero, we use a mask that is the broadcast result of an array 

896 # of True values shaped like `b` with `near_zeros`. 

897 mask = np.ones_like(b, dtype=bool) & near_zeros 

898 q[mask] = 0 

899 

900 x = np.fft.ifft(q, axis=-1) 

901 if not (np.iscomplexobj(c) or np.iscomplexobj(b)): 

902 x = x.real 

903 if outaxis != -1: 

904 x = np.rollaxis(x, -1, outaxis) 

905 return x 

906 

907 

908# matrix inversion 

909def inv(a, overwrite_a=False, check_finite=True): 

910 """ 

911 Compute the inverse of a matrix. 

912 

913 Parameters 

914 ---------- 

915 a : array_like 

916 Square matrix to be inverted. 

917 overwrite_a : bool, optional 

918 Discard data in `a` (may improve performance). Default is False. 

919 check_finite : bool, optional 

920 Whether to check that the input matrix contains only finite numbers. 

921 Disabling may give a performance gain, but may result in problems 

922 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

923 

924 Returns 

925 ------- 

926 ainv : ndarray 

927 Inverse of the matrix `a`. 

928 

929 Raises 

930 ------ 

931 LinAlgError 

932 If `a` is singular. 

933 ValueError 

934 If `a` is not square, or not 2D. 

935 

936 Examples 

937 -------- 

938 >>> from scipy import linalg 

939 >>> a = np.array([[1., 2.], [3., 4.]]) 

940 >>> linalg.inv(a) 

941 array([[-2. , 1. ], 

942 [ 1.5, -0.5]]) 

943 >>> np.dot(a, linalg.inv(a)) 

944 array([[ 1., 0.], 

945 [ 0., 1.]]) 

946 

947 """ 

948 a1 = _asarray_validated(a, check_finite=check_finite) 

949 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]: 

950 raise ValueError('expected square matrix') 

951 overwrite_a = overwrite_a or _datacopied(a1, a) 

952 # XXX: I found no advantage or disadvantage of using finv. 

953# finv, = get_flinalg_funcs(('inv',),(a1,)) 

954# if finv is not None: 

955# a_inv,info = finv(a1,overwrite_a=overwrite_a) 

956# if info==0: 

957# return a_inv 

958# if info>0: raise LinAlgError, "singular matrix" 

959# if info<0: raise ValueError('illegal value in %d-th argument of ' 

960# 'internal inv.getrf|getri'%(-info)) 

961 getrf, getri, getri_lwork = get_lapack_funcs(('getrf', 'getri', 

962 'getri_lwork'), 

963 (a1,)) 

964 lu, piv, info = getrf(a1, overwrite_a=overwrite_a) 

965 if info == 0: 

966 lwork = _compute_lwork(getri_lwork, a1.shape[0]) 

967 

968 # XXX: the following line fixes curious SEGFAULT when 

969 # benchmarking 500x500 matrix inverse. This seems to 

970 # be a bug in LAPACK ?getri routine because if lwork is 

971 # minimal (when using lwork[0] instead of lwork[1]) then 

972 # all tests pass. Further investigation is required if 

973 # more such SEGFAULTs occur. 

974 lwork = int(1.01 * lwork) 

975 inv_a, info = getri(lu, piv, lwork=lwork, overwrite_lu=1) 

976 if info > 0: 

977 raise LinAlgError("singular matrix") 

978 if info < 0: 

979 raise ValueError('illegal value in %d-th argument of internal ' 

980 'getrf|getri' % -info) 

981 return inv_a 

982 

983 

984# Determinant 

985 

986def det(a, overwrite_a=False, check_finite=True): 

987 """ 

988 Compute the determinant of a matrix 

989 

990 The determinant of a square matrix is a value derived arithmetically 

991 from the coefficients of the matrix. 

992 

993 The determinant for a 3x3 matrix, for example, is computed as follows:: 

994 

995 a b c 

996 d e f = A 

997 g h i 

998 

999 det(A) = a*e*i + b*f*g + c*d*h - c*e*g - b*d*i - a*f*h 

1000 

1001 Parameters 

1002 ---------- 

1003 a : (M, M) array_like 

1004 A square matrix. 

1005 overwrite_a : bool, optional 

1006 Allow overwriting data in a (may enhance performance). 

1007 check_finite : bool, optional 

1008 Whether to check that the input matrix contains only finite numbers. 

1009 Disabling may give a performance gain, but may result in problems 

1010 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

1011 

1012 Returns 

1013 ------- 

1014 det : float or complex 

1015 Determinant of `a`. 

1016 

1017 Notes 

1018 ----- 

1019 The determinant is computed via LU factorization, LAPACK routine z/dgetrf. 

1020 

1021 Examples 

1022 -------- 

1023 >>> from scipy import linalg 

1024 >>> a = np.array([[1,2,3], [4,5,6], [7,8,9]]) 

1025 >>> linalg.det(a) 

1026 0.0 

1027 >>> a = np.array([[0,2,3], [4,5,6], [7,8,9]]) 

1028 >>> linalg.det(a) 

1029 3.0 

1030 

1031 """ 

1032 a1 = _asarray_validated(a, check_finite=check_finite) 

1033 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]: 

1034 raise ValueError('expected square matrix') 

1035 overwrite_a = overwrite_a or _datacopied(a1, a) 

1036 fdet, = get_flinalg_funcs(('det',), (a1,)) 

1037 a_det, info = fdet(a1, overwrite_a=overwrite_a) 

1038 if info < 0: 

1039 raise ValueError('illegal value in %d-th argument of internal ' 

1040 'det.getrf' % -info) 

1041 return a_det 

1042 

1043 

1044# Linear Least Squares 

1045def lstsq(a, b, cond=None, overwrite_a=False, overwrite_b=False, 

1046 check_finite=True, lapack_driver=None): 

1047 """ 

1048 Compute least-squares solution to equation Ax = b. 

1049 

1050 Compute a vector x such that the 2-norm ``|b - A x|`` is minimized. 

1051 

1052 Parameters 

1053 ---------- 

1054 a : (M, N) array_like 

1055 Left-hand side array 

1056 b : (M,) or (M, K) array_like 

1057 Right hand side array 

1058 cond : float, optional 

1059 Cutoff for 'small' singular values; used to determine effective 

1060 rank of a. Singular values smaller than 

1061 ``rcond * largest_singular_value`` are considered zero. 

1062 overwrite_a : bool, optional 

1063 Discard data in `a` (may enhance performance). Default is False. 

1064 overwrite_b : bool, optional 

1065 Discard data in `b` (may enhance performance). Default is False. 

1066 check_finite : bool, optional 

1067 Whether to check that the input matrices contain only finite numbers. 

1068 Disabling may give a performance gain, but may result in problems 

1069 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

1070 lapack_driver : str, optional 

1071 Which LAPACK driver is used to solve the least-squares problem. 

1072 Options are ``'gelsd'``, ``'gelsy'``, ``'gelss'``. Default 

1073 (``'gelsd'``) is a good choice. However, ``'gelsy'`` can be slightly 

1074 faster on many problems. ``'gelss'`` was used historically. It is 

1075 generally slow but uses less memory. 

1076 

1077 .. versionadded:: 0.17.0 

1078 

1079 Returns 

1080 ------- 

1081 x : (N,) or (N, K) ndarray 

1082 Least-squares solution. Return shape matches shape of `b`. 

1083 residues : (K,) ndarray or float 

1084 Square of the 2-norm for each column in ``b - a x``, if ``M > N`` and 

1085 ``ndim(A) == n`` (returns a scalar if b is 1-D). Otherwise a 

1086 (0,)-shaped array is returned. 

1087 rank : int 

1088 Effective rank of `a`. 

1089 s : (min(M, N),) ndarray or None 

1090 Singular values of `a`. The condition number of a is 

1091 ``abs(s[0] / s[-1])``. 

1092 

1093 Raises 

1094 ------ 

1095 LinAlgError 

1096 If computation does not converge. 

1097 

1098 ValueError 

1099 When parameters are not compatible. 

1100 

1101 See Also 

1102 -------- 

1103 scipy.optimize.nnls : linear least squares with non-negativity constraint 

1104 

1105 Notes 

1106 ----- 

1107 When ``'gelsy'`` is used as a driver, `residues` is set to a (0,)-shaped 

1108 array and `s` is always ``None``. 

1109 

1110 Examples 

1111 -------- 

1112 >>> from scipy.linalg import lstsq 

1113 >>> import matplotlib.pyplot as plt 

1114 

1115 Suppose we have the following data: 

1116 

1117 >>> x = np.array([1, 2.5, 3.5, 4, 5, 7, 8.5]) 

1118 >>> y = np.array([0.3, 1.1, 1.5, 2.0, 3.2, 6.6, 8.6]) 

1119 

1120 We want to fit a quadratic polynomial of the form ``y = a + b*x**2`` 

1121 to this data. We first form the "design matrix" M, with a constant 

1122 column of 1s and a column containing ``x**2``: 

1123 

1124 >>> M = x[:, np.newaxis]**[0, 2] 

1125 >>> M 

1126 array([[ 1. , 1. ], 

1127 [ 1. , 6.25], 

1128 [ 1. , 12.25], 

1129 [ 1. , 16. ], 

1130 [ 1. , 25. ], 

1131 [ 1. , 49. ], 

1132 [ 1. , 72.25]]) 

1133 

1134 We want to find the least-squares solution to ``M.dot(p) = y``, 

1135 where ``p`` is a vector with length 2 that holds the parameters 

1136 ``a`` and ``b``. 

1137 

1138 >>> p, res, rnk, s = lstsq(M, y) 

1139 >>> p 

1140 array([ 0.20925829, 0.12013861]) 

1141 

1142 Plot the data and the fitted curve. 

1143 

1144 >>> plt.plot(x, y, 'o', label='data') 

1145 >>> xx = np.linspace(0, 9, 101) 

1146 >>> yy = p[0] + p[1]*xx**2 

1147 >>> plt.plot(xx, yy, label='least squares fit, $y = a + bx^2$') 

1148 >>> plt.xlabel('x') 

1149 >>> plt.ylabel('y') 

1150 >>> plt.legend(framealpha=1, shadow=True) 

1151 >>> plt.grid(alpha=0.25) 

1152 >>> plt.show() 

1153 

1154 """ 

1155 a1 = _asarray_validated(a, check_finite=check_finite) 

1156 b1 = _asarray_validated(b, check_finite=check_finite) 

1157 if len(a1.shape) != 2: 

1158 raise ValueError('Input array a should be 2D') 

1159 m, n = a1.shape 

1160 if len(b1.shape) == 2: 

1161 nrhs = b1.shape[1] 

1162 else: 

1163 nrhs = 1 

1164 if m != b1.shape[0]: 

1165 raise ValueError('Shape mismatch: a and b should have the same number' 

1166 ' of rows ({} != {}).'.format(m, b1.shape[0])) 

1167 if m == 0 or n == 0: # Zero-sized problem, confuses LAPACK 

1168 x = np.zeros((n,) + b1.shape[1:], dtype=np.common_type(a1, b1)) 

1169 if n == 0: 

1170 residues = np.linalg.norm(b1, axis=0)**2 

1171 else: 

1172 residues = np.empty((0,)) 

1173 return x, residues, 0, np.empty((0,)) 

1174 

1175 driver = lapack_driver 

1176 if driver is None: 

1177 driver = lstsq.default_lapack_driver 

1178 if driver not in ('gelsd', 'gelsy', 'gelss'): 

1179 raise ValueError('LAPACK driver "%s" is not found' % driver) 

1180 

1181 lapack_func, lapack_lwork = get_lapack_funcs((driver, 

1182 '%s_lwork' % driver), 

1183 (a1, b1)) 

1184 real_data = True if (lapack_func.dtype.kind == 'f') else False 

1185 

1186 if m < n: 

1187 # need to extend b matrix as it will be filled with 

1188 # a larger solution matrix 

1189 if len(b1.shape) == 2: 

1190 b2 = np.zeros((n, nrhs), dtype=lapack_func.dtype) 

1191 b2[:m, :] = b1 

1192 else: 

1193 b2 = np.zeros(n, dtype=lapack_func.dtype) 

1194 b2[:m] = b1 

1195 b1 = b2 

1196 

1197 overwrite_a = overwrite_a or _datacopied(a1, a) 

1198 overwrite_b = overwrite_b or _datacopied(b1, b) 

1199 

1200 if cond is None: 

1201 cond = np.finfo(lapack_func.dtype).eps 

1202 

1203 if driver in ('gelss', 'gelsd'): 

1204 if driver == 'gelss': 

1205 lwork = _compute_lwork(lapack_lwork, m, n, nrhs, cond) 

1206 v, x, s, rank, work, info = lapack_func(a1, b1, cond, lwork, 

1207 overwrite_a=overwrite_a, 

1208 overwrite_b=overwrite_b) 

1209 

1210 elif driver == 'gelsd': 

1211 if real_data: 

1212 lwork, iwork = _compute_lwork(lapack_lwork, m, n, nrhs, cond) 

1213 x, s, rank, info = lapack_func(a1, b1, lwork, 

1214 iwork, cond, False, False) 

1215 else: # complex data 

1216 lwork, rwork, iwork = _compute_lwork(lapack_lwork, m, n, 

1217 nrhs, cond) 

1218 x, s, rank, info = lapack_func(a1, b1, lwork, rwork, iwork, 

1219 cond, False, False) 

1220 if info > 0: 

1221 raise LinAlgError("SVD did not converge in Linear Least Squares") 

1222 if info < 0: 

1223 raise ValueError('illegal value in %d-th argument of internal %s' 

1224 % (-info, lapack_driver)) 

1225 resids = np.asarray([], dtype=x.dtype) 

1226 if m > n: 

1227 x1 = x[:n] 

1228 if rank == n: 

1229 resids = np.sum(np.abs(x[n:])**2, axis=0) 

1230 x = x1 

1231 return x, resids, rank, s 

1232 

1233 elif driver == 'gelsy': 

1234 lwork = _compute_lwork(lapack_lwork, m, n, nrhs, cond) 

1235 jptv = np.zeros((a1.shape[1], 1), dtype=np.int32) 

1236 v, x, j, rank, info = lapack_func(a1, b1, jptv, cond, 

1237 lwork, False, False) 

1238 if info < 0: 

1239 raise ValueError("illegal value in %d-th argument of internal " 

1240 "gelsy" % -info) 

1241 if m > n: 

1242 x1 = x[:n] 

1243 x = x1 

1244 return x, np.array([], x.dtype), rank, None 

1245 

1246 

1247lstsq.default_lapack_driver = 'gelsd' 

1248 

1249 

1250def pinv(a, cond=None, rcond=None, return_rank=False, check_finite=True): 

1251 """ 

1252 Compute the (Moore-Penrose) pseudo-inverse of a matrix. 

1253 

1254 Calculate a generalized inverse of a matrix using a least-squares 

1255 solver. 

1256 

1257 Parameters 

1258 ---------- 

1259 a : (M, N) array_like 

1260 Matrix to be pseudo-inverted. 

1261 cond, rcond : float, optional 

1262 Cutoff factor for 'small' singular values. In `lstsq`, 

1263 singular values less than ``cond*largest_singular_value`` will be 

1264 considered as zero. If both are omitted, the default value 

1265 ``max(M, N) * eps`` is passed to `lstsq` where ``eps`` is the 

1266 corresponding machine precision value of the datatype of ``a``. 

1267 

1268 .. versionchanged:: 1.3.0 

1269 Previously the default cutoff value was just `eps` without the 

1270 factor ``max(M, N)``. 

1271 

1272 return_rank : bool, optional 

1273 if True, return the effective rank of the matrix 

1274 check_finite : bool, optional 

1275 Whether to check that the input matrix contains only finite numbers. 

1276 Disabling may give a performance gain, but may result in problems 

1277 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

1278 

1279 Returns 

1280 ------- 

1281 B : (N, M) ndarray 

1282 The pseudo-inverse of matrix `a`. 

1283 rank : int 

1284 The effective rank of the matrix. Returned if return_rank == True 

1285 

1286 Raises 

1287 ------ 

1288 LinAlgError 

1289 If computation does not converge. 

1290 

1291 Examples 

1292 -------- 

1293 >>> from scipy import linalg 

1294 >>> a = np.random.randn(9, 6) 

1295 >>> B = linalg.pinv(a) 

1296 >>> np.allclose(a, np.dot(a, np.dot(B, a))) 

1297 True 

1298 >>> np.allclose(B, np.dot(B, np.dot(a, B))) 

1299 True 

1300 

1301 """ 

1302 a = _asarray_validated(a, check_finite=check_finite) 

1303 b = np.identity(a.shape[0], dtype=a.dtype) 

1304 

1305 if rcond is not None: 

1306 cond = rcond 

1307 

1308 if cond is None: 

1309 cond = max(a.shape) * np.spacing(a.real.dtype.type(1)) 

1310 

1311 x, resids, rank, s = lstsq(a, b, cond=cond, check_finite=False) 

1312 

1313 if return_rank: 

1314 return x, rank 

1315 else: 

1316 return x 

1317 

1318 

1319def pinv2(a, cond=None, rcond=None, return_rank=False, check_finite=True): 

1320 """ 

1321 Compute the (Moore-Penrose) pseudo-inverse of a matrix. 

1322 

1323 Calculate a generalized inverse of a matrix using its 

1324 singular-value decomposition and including all 'large' singular 

1325 values. 

1326 

1327 Parameters 

1328 ---------- 

1329 a : (M, N) array_like 

1330 Matrix to be pseudo-inverted. 

1331 cond, rcond : float or None 

1332 Cutoff for 'small' singular values; singular values smaller than this 

1333 value are considered as zero. If both are omitted, the default value 

1334 ``max(M,N)*largest_singular_value*eps`` is used where ``eps`` is the 

1335 machine precision value of the datatype of ``a``. 

1336 

1337 .. versionchanged:: 1.3.0 

1338 Previously the default cutoff value was just ``eps*f`` where ``f`` 

1339 was ``1e3`` for single precision and ``1e6`` for double precision. 

1340 

1341 return_rank : bool, optional 

1342 If True, return the effective rank of the matrix. 

1343 check_finite : bool, optional 

1344 Whether to check that the input matrix contains only finite numbers. 

1345 Disabling may give a performance gain, but may result in problems 

1346 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

1347 

1348 Returns 

1349 ------- 

1350 B : (N, M) ndarray 

1351 The pseudo-inverse of matrix `a`. 

1352 rank : int 

1353 The effective rank of the matrix. Returned if `return_rank` is True. 

1354 

1355 Raises 

1356 ------ 

1357 LinAlgError 

1358 If SVD computation does not converge. 

1359 

1360 Examples 

1361 -------- 

1362 >>> from scipy import linalg 

1363 >>> a = np.random.randn(9, 6) 

1364 >>> B = linalg.pinv2(a) 

1365 >>> np.allclose(a, np.dot(a, np.dot(B, a))) 

1366 True 

1367 >>> np.allclose(B, np.dot(B, np.dot(a, B))) 

1368 True 

1369 

1370 """ 

1371 a = _asarray_validated(a, check_finite=check_finite) 

1372 u, s, vh = decomp_svd.svd(a, full_matrices=False, check_finite=False) 

1373 

1374 if rcond is not None: 

1375 cond = rcond 

1376 if cond in [None, -1]: 

1377 t = u.dtype.char.lower() 

1378 cond = np.max(s) * max(a.shape) * np.finfo(t).eps 

1379 

1380 rank = np.sum(s > cond) 

1381 

1382 u = u[:, :rank] 

1383 u /= s[:rank] 

1384 B = np.transpose(np.conjugate(np.dot(u, vh[:rank]))) 

1385 

1386 if return_rank: 

1387 return B, rank 

1388 else: 

1389 return B 

1390 

1391 

1392def pinvh(a, cond=None, rcond=None, lower=True, return_rank=False, 

1393 check_finite=True): 

1394 """ 

1395 Compute the (Moore-Penrose) pseudo-inverse of a Hermitian matrix. 

1396 

1397 Calculate a generalized inverse of a Hermitian or real symmetric matrix 

1398 using its eigenvalue decomposition and including all eigenvalues with 

1399 'large' absolute value. 

1400 

1401 Parameters 

1402 ---------- 

1403 a : (N, N) array_like 

1404 Real symmetric or complex hermetian matrix to be pseudo-inverted 

1405 cond, rcond : float or None 

1406 Cutoff for 'small' singular values; singular values smaller than this 

1407 value are considered as zero. If both are omitted, the default 

1408 ``max(M,N)*largest_eigenvalue*eps`` is used where ``eps`` is the 

1409 machine precision value of the datatype of ``a``. 

1410 

1411 .. versionchanged:: 1.3.0 

1412 Previously the default cutoff value was just ``eps*f`` where ``f`` 

1413 was ``1e3`` for single precision and ``1e6`` for double precision. 

1414 

1415 lower : bool, optional 

1416 Whether the pertinent array data is taken from the lower or upper 

1417 triangle of `a`. (Default: lower) 

1418 return_rank : bool, optional 

1419 If True, return the effective rank of the matrix. 

1420 check_finite : bool, optional 

1421 Whether to check that the input matrix contains only finite numbers. 

1422 Disabling may give a performance gain, but may result in problems 

1423 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

1424 

1425 Returns 

1426 ------- 

1427 B : (N, N) ndarray 

1428 The pseudo-inverse of matrix `a`. 

1429 rank : int 

1430 The effective rank of the matrix. Returned if `return_rank` is True. 

1431 

1432 Raises 

1433 ------ 

1434 LinAlgError 

1435 If eigenvalue does not converge 

1436 

1437 Examples 

1438 -------- 

1439 >>> from scipy.linalg import pinvh 

1440 >>> a = np.random.randn(9, 6) 

1441 >>> a = np.dot(a, a.T) 

1442 >>> B = pinvh(a) 

1443 >>> np.allclose(a, np.dot(a, np.dot(B, a))) 

1444 True 

1445 >>> np.allclose(B, np.dot(B, np.dot(a, B))) 

1446 True 

1447 

1448 """ 

1449 a = _asarray_validated(a, check_finite=check_finite) 

1450 s, u = decomp.eigh(a, lower=lower, check_finite=False) 

1451 

1452 if rcond is not None: 

1453 cond = rcond 

1454 if cond in [None, -1]: 

1455 t = u.dtype.char.lower() 

1456 cond = np.max(np.abs(s)) * max(a.shape) * np.finfo(t).eps 

1457 

1458 # For Hermitian matrices, singular values equal abs(eigenvalues) 

1459 above_cutoff = (abs(s) > cond) 

1460 psigma_diag = 1.0 / s[above_cutoff] 

1461 u = u[:, above_cutoff] 

1462 

1463 B = np.dot(u * psigma_diag, np.conjugate(u).T) 

1464 

1465 if return_rank: 

1466 return B, len(psigma_diag) 

1467 else: 

1468 return B 

1469 

1470 

1471def matrix_balance(A, permute=True, scale=True, separate=False, 

1472 overwrite_a=False): 

1473 """ 

1474 Compute a diagonal similarity transformation for row/column balancing. 

1475 

1476 The balancing tries to equalize the row and column 1-norms by applying 

1477 a similarity transformation such that the magnitude variation of the 

1478 matrix entries is reflected to the scaling matrices. 

1479 

1480 Moreover, if enabled, the matrix is first permuted to isolate the upper 

1481 triangular parts of the matrix and, again if scaling is also enabled, 

1482 only the remaining subblocks are subjected to scaling. 

1483 

1484 The balanced matrix satisfies the following equality 

1485 

1486 .. math:: 

1487 

1488 B = T^{-1} A T 

1489 

1490 The scaling coefficients are approximated to the nearest power of 2 

1491 to avoid round-off errors. 

1492 

1493 Parameters 

1494 ---------- 

1495 A : (n, n) array_like 

1496 Square data matrix for the balancing. 

1497 permute : bool, optional 

1498 The selector to define whether permutation of A is also performed 

1499 prior to scaling. 

1500 scale : bool, optional 

1501 The selector to turn on and off the scaling. If False, the matrix 

1502 will not be scaled. 

1503 separate : bool, optional 

1504 This switches from returning a full matrix of the transformation 

1505 to a tuple of two separate 1-D permutation and scaling arrays. 

1506 overwrite_a : bool, optional 

1507 This is passed to xGEBAL directly. Essentially, overwrites the result 

1508 to the data. It might increase the space efficiency. See LAPACK manual 

1509 for details. This is False by default. 

1510 

1511 Returns 

1512 ------- 

1513 B : (n, n) ndarray 

1514 Balanced matrix 

1515 T : (n, n) ndarray 

1516 A possibly permuted diagonal matrix whose nonzero entries are 

1517 integer powers of 2 to avoid numerical truncation errors. 

1518 scale, perm : (n,) ndarray 

1519 If ``separate`` keyword is set to True then instead of the array 

1520 ``T`` above, the scaling and the permutation vectors are given 

1521 separately as a tuple without allocating the full array ``T``. 

1522 

1523 Notes 

1524 ----- 

1525 

1526 This algorithm is particularly useful for eigenvalue and matrix 

1527 decompositions and in many cases it is already called by various 

1528 LAPACK routines. 

1529 

1530 The algorithm is based on the well-known technique of [1]_ and has 

1531 been modified to account for special cases. See [2]_ for details 

1532 which have been implemented since LAPACK v3.5.0. Before this version 

1533 there are corner cases where balancing can actually worsen the 

1534 conditioning. See [3]_ for such examples. 

1535 

1536 The code is a wrapper around LAPACK's xGEBAL routine family for matrix 

1537 balancing. 

1538 

1539 .. versionadded:: 0.19.0 

1540 

1541 Examples 

1542 -------- 

1543 >>> from scipy import linalg 

1544 >>> x = np.array([[1,2,0], [9,1,0.01], [1,2,10*np.pi]]) 

1545 

1546 >>> y, permscale = linalg.matrix_balance(x) 

1547 >>> np.abs(x).sum(axis=0) / np.abs(x).sum(axis=1) 

1548 array([ 3.66666667, 0.4995005 , 0.91312162]) 

1549 

1550 >>> np.abs(y).sum(axis=0) / np.abs(y).sum(axis=1) 

1551 array([ 1.2 , 1.27041742, 0.92658316]) # may vary 

1552 

1553 >>> permscale # only powers of 2 (0.5 == 2^(-1)) 

1554 array([[ 0.5, 0. , 0. ], # may vary 

1555 [ 0. , 1. , 0. ], 

1556 [ 0. , 0. , 1. ]]) 

1557 

1558 References 

1559 ---------- 

1560 .. [1] : B.N. Parlett and C. Reinsch, "Balancing a Matrix for 

1561 Calculation of Eigenvalues and Eigenvectors", Numerische Mathematik, 

1562 Vol.13(4), 1969, DOI:10.1007/BF02165404 

1563 

1564 .. [2] : R. James, J. Langou, B.R. Lowery, "On matrix balancing and 

1565 eigenvector computation", 2014, Available online: 

1566 https://arxiv.org/abs/1401.5766 

1567 

1568 .. [3] : D.S. Watkins. A case where balancing is harmful. 

1569 Electron. Trans. Numer. Anal, Vol.23, 2006. 

1570 

1571 """ 

1572 

1573 A = np.atleast_2d(_asarray_validated(A, check_finite=True)) 

1574 

1575 if not np.equal(*A.shape): 

1576 raise ValueError('The data matrix for balancing should be square.') 

1577 

1578 gebal = get_lapack_funcs(('gebal'), (A,)) 

1579 B, lo, hi, ps, info = gebal(A, scale=scale, permute=permute, 

1580 overwrite_a=overwrite_a) 

1581 

1582 if info < 0: 

1583 raise ValueError('xGEBAL exited with the internal error ' 

1584 '"illegal value in argument number {}.". See ' 

1585 'LAPACK documentation for the xGEBAL error codes.' 

1586 ''.format(-info)) 

1587 

1588 # Separate the permutations from the scalings and then convert to int 

1589 scaling = np.ones_like(ps, dtype=float) 

1590 scaling[lo:hi+1] = ps[lo:hi+1] 

1591 

1592 # gebal uses 1-indexing 

1593 ps = ps.astype(int, copy=False) - 1 

1594 n = A.shape[0] 

1595 perm = np.arange(n) 

1596 

1597 # LAPACK permutes with the ordering n --> hi, then 0--> lo 

1598 if hi < n: 

1599 for ind, x in enumerate(ps[hi+1:][::-1], 1): 

1600 if n-ind == x: 

1601 continue 

1602 perm[[x, n-ind]] = perm[[n-ind, x]] 

1603 

1604 if lo > 0: 

1605 for ind, x in enumerate(ps[:lo]): 

1606 if ind == x: 

1607 continue 

1608 perm[[x, ind]] = perm[[ind, x]] 

1609 

1610 if separate: 

1611 return B, (scaling, perm) 

1612 

1613 # get the inverse permutation 

1614 iperm = np.empty_like(perm) 

1615 iperm[perm] = np.arange(n) 

1616 

1617 return B, np.diag(scaling)[iperm, :]