Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/linalg/_basic.py: 7%

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_py 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', 'pinvh', 'matrix_balance', 'matmul_toeplitz'] 

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, 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, default: False, deprecated 

72 Assume `a` is symmetric and positive definite. 

73 

74 .. deprecated:: 0.19.0 

75 This keyword is deprecated and should be replaced by using 

76 ``assume_a = 'pos'``. `sym_pos` will be removed in SciPy 1.11.0. 

77 

78 lower : bool, default: False 

79 Ignored if ``assume_a == 'gen'`` (the default). If True, the 

80 calculation uses only the data in the lower triangle of `a`; 

81 entries above the diagonal are ignored. If False (default), the 

82 calculation uses only the data in the upper triangle of `a`; entries 

83 below the diagonal are ignored. 

84 overwrite_a : bool, default: False 

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

86 overwrite_b : bool, default: False 

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

88 check_finite : bool, default: True 

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

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

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

92 assume_a : str, {'gen', 'sym', 'her', 'pos'} 

93 Valid entries are explained above. 

94 transposed : bool, default: False 

95 If True, solve ``a.T @ x == b``. Raises `NotImplementedError` 

96 for complex `a`. 

97 

98 Returns 

99 ------- 

100 x : (N, NRHS) ndarray 

101 The solution array. 

102 

103 Raises 

104 ------ 

105 ValueError 

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

107 LinAlgError 

108 If the matrix is singular. 

109 LinAlgWarning 

110 If an ill-conditioned input a is detected. 

111 NotImplementedError 

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

113 

114 Notes 

115 ----- 

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

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

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

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

120 

121 The generic, symmetric, Hermitian and positive definite solutions are 

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

123 LAPACK respectively. 

124 

125 Examples 

126 -------- 

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

128 

129 >>> import numpy as np 

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

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

132 >>> from scipy import linalg 

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

134 >>> x 

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

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

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

138 

139 """ 

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

141 b_is_1D = False 

142 

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

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

145 n = a1.shape[0] 

146 

147 overwrite_a = overwrite_a or _datacopied(a1, a) 

148 overwrite_b = overwrite_b or _datacopied(b1, b) 

149 

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

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

152 

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

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

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

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

157 'input a') 

158 

159 # accommodate empty arrays 

160 if b1.size == 0: 

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

162 

163 # regularize 1-D b arrays to 2D 

164 if b1.ndim == 1: 

165 if n == 1: 

166 b1 = b1[None, :] 

167 else: 

168 b1 = b1[:, None] 

169 b_is_1D = True 

170 

171 # Backwards compatibility - old keyword. 

172 if sym_pos: 

173 message = ("The 'sym_pos' keyword is deprecated and should be " 

174 "replaced by using 'assume_a = \"pos\"'. 'sym_pos' will be" 

175 " removed in SciPy 1.11.0.") 

176 warn(message, DeprecationWarning, stacklevel=2) 

177 assume_a = 'pos' 

178 

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

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

181 ''.format(assume_a)) 

182 

183 # for a real matrix, describe it as "symmetric", not "hermitian" 

184 # (lapack doesn't know what to do with real hermitian matrices) 

185 if assume_a == 'her' and not np.iscomplexobj(a1): 

186 assume_a = 'sym' 

187 

188 # Get the correct lamch function. 

189 # The LAMCH functions only exists for S and D 

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

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

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

193 else: 

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

195 

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

197 # lansy, lanpo, lanhe. 

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

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

200 

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

202 # we can collect them all in this one call 

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

204 # the I-norm should be used 

205 if transposed: 

206 trans = 1 

207 norm = 'I' 

208 if np.iscomplexobj(a1): 

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

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

211 'for complex matrices.') 

212 else: 

213 trans = 0 

214 norm = '1' 

215 

216 anorm = lange(norm, a1) 

217 

218 # Generalized case 'gesv' 

219 if assume_a == 'gen': 

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

221 (a1, b1)) 

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

223 _solve_check(n, info) 

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

225 trans=trans, overwrite_b=overwrite_b) 

226 _solve_check(n, info) 

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

228 # Hermitian case 'hesv' 

229 elif assume_a == 'her': 

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

231 'hesv_lwork'), (a1, b1)) 

232 lwork = _compute_lwork(hesv_lw, n, lower) 

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

234 lower=lower, 

235 overwrite_a=overwrite_a, 

236 overwrite_b=overwrite_b) 

237 _solve_check(n, info) 

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

239 # Symmetric case 'sysv' 

240 elif assume_a == 'sym': 

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

242 'sysv_lwork'), (a1, b1)) 

243 lwork = _compute_lwork(sysv_lw, n, lower) 

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

245 lower=lower, 

246 overwrite_a=overwrite_a, 

247 overwrite_b=overwrite_b) 

248 _solve_check(n, info) 

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

250 # Positive definite case 'posv' 

251 else: 

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

253 (a1, b1)) 

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

255 overwrite_a=overwrite_a, 

256 overwrite_b=overwrite_b) 

257 _solve_check(n, info) 

258 rcond, info = pocon(lu, anorm) 

259 

260 _solve_check(n, info, lamch, rcond) 

261 

262 if b_is_1D: 

263 x = x.ravel() 

264 

265 return x 

266 

267 

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

269 overwrite_b=False, check_finite=True): 

270 """ 

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

272 

273 Parameters 

274 ---------- 

275 a : (M, M) array_like 

276 A triangular matrix 

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

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

279 lower : bool, optional 

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

281 Default is to use upper triangle. 

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

283 Type of system to solve: 

284 

285 ======== ========= 

286 trans system 

287 ======== ========= 

288 0 or 'N' a x = b 

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

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

291 ======== ========= 

292 unit_diagonal : bool, optional 

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

294 will not be referenced. 

295 overwrite_b : bool, optional 

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

297 check_finite : bool, optional 

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

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

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

301 

302 Returns 

303 ------- 

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

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

306 

307 Raises 

308 ------ 

309 LinAlgError 

310 If `a` is singular 

311 

312 Notes 

313 ----- 

314 .. versionadded:: 0.9.0 

315 

316 Examples 

317 -------- 

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

319 

320 [3 0 0 0] [4] 

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

322 [1 0 1 0] [4] 

323 [1 1 1 1] [2] 

324 

325 >>> import numpy as np 

326 >>> from scipy.linalg import solve_triangular 

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

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

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

330 >>> x 

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

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

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

334 

335 """ 

336 

337 a1 = _asarray_validated(a, check_finite=check_finite) 

338 b1 = _asarray_validated(b, check_finite=check_finite) 

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

340 raise ValueError('expected square matrix') 

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

342 raise ValueError('shapes of a {} and b {} are incompatible' 

343 .format(a1.shape, b1.shape)) 

344 overwrite_b = overwrite_b or _datacopied(b1, b) 

345 

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

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

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

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

350 trans=trans, unitdiag=unit_diagonal) 

351 else: 

352 # transposed system is solved since trtrs expects Fortran ordering 

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

354 trans=not trans, unitdiag=unit_diagonal) 

355 

356 if info == 0: 

357 return x 

358 if info > 0: 

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

360 (info-1)) 

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

362 (-info)) 

363 

364 

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

366 check_finite=True): 

367 """ 

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

369 

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

371 

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

373 

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

375 

376 * a01 a12 a23 a34 a45 

377 a00 a11 a22 a33 a44 a55 

378 a10 a21 a32 a43 a54 * 

379 a20 a31 a42 a53 * * 

380 

381 Parameters 

382 ---------- 

383 (l, u) : (integer, integer) 

384 Number of non-zero lower and upper diagonals 

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

386 Banded matrix 

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

388 Right-hand side 

389 overwrite_ab : bool, optional 

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

391 overwrite_b : bool, optional 

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

393 check_finite : bool, optional 

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

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

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

397 

398 Returns 

399 ------- 

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

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

402 shape of `b`. 

403 

404 Examples 

405 -------- 

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

407 

408 [5 2 -1 0 0] [0] 

409 [1 4 2 -1 0] [1] 

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

411 [0 0 1 2 2] [2] 

412 [0 0 0 1 1] [3] 

413 

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

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

416 

417 [* * -1 -1 -1] 

418 ab = [* 2 2 2 2] 

419 [5 4 3 2 1] 

420 [1 1 1 1 *] 

421 

422 >>> import numpy as np 

423 >>> from scipy.linalg import solve_banded 

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

425 ... [0, 2, 2, 2, 2], 

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

427 ... [1, 1, 1, 1, 0]]) 

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

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

430 >>> x 

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

432 

433 """ 

434 

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

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

437 # Validate shapes. 

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

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

440 (nlower, nupper) = l_and_u 

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

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

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

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

445 

446 overwrite_b = overwrite_b or _datacopied(b1, b) 

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

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

449 b2 /= a1[1, 0] 

450 return b2 

451 if nlower == nupper == 1: 

452 overwrite_ab = overwrite_ab or _datacopied(a1, ab) 

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

454 du = a1[0, 1:] 

455 d = a1[1, :] 

456 dl = a1[2, :-1] 

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

458 overwrite_ab, overwrite_b) 

459 else: 

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

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

462 a2[nlower:, :] = a1 

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

464 overwrite_b=overwrite_b) 

465 if info == 0: 

466 return x 

467 if info > 0: 

468 raise LinAlgError("singular matrix") 

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

470 'gbsv/gtsv' % -info) 

471 

472 

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

474 check_finite=True): 

475 """ 

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

477 

478 Uses Thomas' Algorithm, which is more efficient than standard LU 

479 factorization, but should only be used for Hermitian positive-definite 

480 matrices. 

481 

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

483 diagonal ordered form: 

484 

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

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

487 

488 Example of `ab` (shape of ``a`` is (6, 6), number of upper diagonals, 

489 ``u`` =2):: 

490 

491 upper form: 

492 * * a02 a13 a24 a35 

493 * a01 a12 a23 a34 a45 

494 a00 a11 a22 a33 a44 a55 

495 

496 lower form: 

497 a00 a11 a22 a33 a44 a55 

498 a10 a21 a32 a43 a54 * 

499 a20 a31 a42 a53 * * 

500 

501 Cells marked with * are not used. 

502 

503 Parameters 

504 ---------- 

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

506 Banded matrix 

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

508 Right-hand side 

509 overwrite_ab : bool, optional 

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

511 overwrite_b : bool, optional 

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

513 lower : bool, optional 

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

515 check_finite : bool, optional 

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

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

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

519 

520 Returns 

521 ------- 

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

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

524 of `b`. 

525 

526 Notes 

527 ----- 

528 In the case of a non-positive definite matrix ``a``, the solver 

529 `solve_banded` may be used. 

530 

531 Examples 

532 -------- 

533 Solve the banded system ``A x = b``, where:: 

534 

535 [ 4 2 -1 0 0 0] [1] 

536 [ 2 5 2 -1 0 0] [2] 

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

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

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

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

541 

542 >>> import numpy as np 

543 >>> from scipy.linalg import solveh_banded 

544 

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

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

547 

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

549 ... [ 2, 2, 2, 2, 2, 0], 

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

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

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

553 >>> x 

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

555 0.34733894]) 

556 

557 

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

559 

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

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

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

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

564 

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

566 

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

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

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

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

571 >>> x 

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

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

574 

575 """ 

576 a1 = _asarray_validated(ab, check_finite=check_finite) 

577 b1 = _asarray_validated(b, check_finite=check_finite) 

578 # Validate shapes. 

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

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

581 

582 overwrite_b = overwrite_b or _datacopied(b1, b) 

583 overwrite_ab = overwrite_ab or _datacopied(a1, ab) 

584 

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

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

587 if lower: 

588 d = a1[0, :].real 

589 e = a1[1, :-1] 

590 else: 

591 d = a1[1, :].real 

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

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

594 overwrite_b) 

595 else: 

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

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

598 overwrite_b=overwrite_b) 

599 if info > 0: 

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

601 if info < 0: 

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

603 'pbsv' % -info) 

604 return x 

605 

606 

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

608 """Solve a Toeplitz system using Levinson Recursion 

609 

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

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

612 assumed. 

613 

614 Parameters 

615 ---------- 

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

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

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

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

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

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

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

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

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

625 check_finite : bool, optional 

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

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

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

629 

630 Returns 

631 ------- 

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

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

634 of `b`. 

635 

636 See Also 

637 -------- 

638 toeplitz : Toeplitz matrix 

639 

640 Notes 

641 ----- 

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

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

644 

645 Examples 

646 -------- 

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

648 

649 [ 1 -1 -2 -3] [1] 

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

651 [ 6 3 1 -1] [2] 

652 [10 6 3 1] [5] 

653 

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

655 row are needed. 

656 

657 >>> import numpy as np 

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

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

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

661 

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

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

664 >>> x 

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

666 

667 Check the result by creating the full Toeplitz matrix and 

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

669 

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

671 >>> T.dot(x) 

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

673 

674 """ 

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

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

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

678 

679 r, c, b, dtype, b_shape = _validate_args_for_toeplitz_ops( 

680 c_or_cr, b, check_finite, keep_b_shape=True) 

681 

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

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

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

685 if b is None: 

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

687 

688 if b.ndim == 1: 

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

690 else: 

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

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

693 x = x.reshape(*b_shape) 

694 

695 return x 

696 

697 

698def _get_axis_len(aname, a, axis): 

699 ax = axis 

700 if ax < 0: 

701 ax += a.ndim 

702 if 0 <= ax < a.ndim: 

703 return a.shape[ax] 

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

705 

706 

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

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

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

710 

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

712 

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

714 calculation is:: 

715 

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

717 

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

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

720 solving the system with the full circulant matrix. 

721 

722 Parameters 

723 ---------- 

724 c : array_like 

725 The coefficients of the circulant matrix. 

726 b : array_like 

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

728 singular : str, optional 

729 This argument controls how a near singular circulant matrix is 

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

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

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

733 tol : float, optional 

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

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

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

737 

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

739 

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

741 of the circulant matrix. 

742 caxis : int 

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

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

745 holds the vectors of circulant coefficients. 

746 baxis : int 

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

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

749 right-hand side vectors. 

750 outaxis : int 

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

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

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

754 

755 Returns 

756 ------- 

757 x : ndarray 

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

759 

760 Raises 

761 ------ 

762 LinAlgError 

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

764 

765 See Also 

766 -------- 

767 circulant : circulant matrix 

768 

769 Notes 

770 ----- 

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

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

773 

774 solve_circulant(c, b) 

775 

776 returns the same result as 

777 

778 solve(circulant(c), b) 

779 

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

781 

782 .. versionadded:: 0.16.0 

783 

784 Examples 

785 -------- 

786 >>> import numpy as np 

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

788 

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

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

791 >>> solve_circulant(c, b) 

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

793 

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

795 

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

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

798 

799 A singular example: 

800 

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

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

803 

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

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

806 

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

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

809 

810 Compare to `scipy.linalg.lstsq`: 

811 

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

813 >>> x 

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

815 

816 A broadcasting example: 

817 

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

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

820 (3, 5). For example, 

821 

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

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

824 

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

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

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

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

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

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

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

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

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

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

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

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

837 

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

839 >>> x.shape 

840 (2, 3, 5) 

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

842 >>> x 

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

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

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

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

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

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

849 

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

851 

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

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

854 

855 """ 

856 c = np.atleast_1d(c) 

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

858 b = np.atleast_1d(b) 

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

860 if nc != nb: 

861 raise ValueError('Shapes of c {} and b {} are incompatible' 

862 .format(c.shape, b.shape)) 

863 

864 fc = np.fft.fft(np.moveaxis(c, caxis, -1), axis=-1) 

865 abs_fc = np.abs(fc) 

866 if tol is None: 

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

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

869 if tol.shape != (): 

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

871 else: 

872 tol = np.atleast_1d(tol) 

873 

874 near_zeros = abs_fc <= tol 

875 is_near_singular = np.any(near_zeros) 

876 if is_near_singular: 

877 if singular == 'raise': 

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

879 else: 

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

881 # division fb/fc below. 

882 fc[near_zeros] = 1 

883 

884 fb = np.fft.fft(np.moveaxis(b, baxis, -1), axis=-1) 

885 

886 q = fb / fc 

887 

888 if is_near_singular: 

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

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

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

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

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

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

895 q[mask] = 0 

896 

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

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

899 x = x.real 

900 if outaxis != -1: 

901 x = np.moveaxis(x, -1, outaxis) 

902 return x 

903 

904 

905# matrix inversion 

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

907 """ 

908 Compute the inverse of a matrix. 

909 

910 Parameters 

911 ---------- 

912 a : array_like 

913 Square matrix to be inverted. 

914 overwrite_a : bool, optional 

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

916 check_finite : bool, optional 

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

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

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

920 

921 Returns 

922 ------- 

923 ainv : ndarray 

924 Inverse of the matrix `a`. 

925 

926 Raises 

927 ------ 

928 LinAlgError 

929 If `a` is singular. 

930 ValueError 

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

932 

933 Examples 

934 -------- 

935 >>> import numpy as np 

936 >>> from scipy import linalg 

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

938 >>> linalg.inv(a) 

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

940 [ 1.5, -0.5]]) 

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

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

943 [ 0., 1.]]) 

944 

945 """ 

946 a1 = _asarray_validated(a, check_finite=check_finite) 

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

948 raise ValueError('expected square matrix') 

949 overwrite_a = overwrite_a or _datacopied(a1, a) 

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

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

952# if finv is not None: 

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

954# if info==0: 

955# return a_inv 

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

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

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

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

960 'getri_lwork'), 

961 (a1,)) 

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

963 if info == 0: 

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

965 

966 # XXX: the following line fixes curious SEGFAULT when 

967 # benchmarking 500x500 matrix inverse. This seems to 

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

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

970 # all tests pass. Further investigation is required if 

971 # more such SEGFAULTs occur. 

972 lwork = int(1.01 * lwork) 

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

974 if info > 0: 

975 raise LinAlgError("singular matrix") 

976 if info < 0: 

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

978 'getrf|getri' % -info) 

979 return inv_a 

980 

981 

982# Determinant 

983 

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

985 """ 

986 Compute the determinant of a matrix 

987 

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

989 from the coefficients of the matrix. 

990 

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

992 

993 a b c 

994 d e f = A 

995 g h i 

996 

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

998 

999 Parameters 

1000 ---------- 

1001 a : (M, M) array_like 

1002 A square matrix. 

1003 overwrite_a : bool, optional 

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

1005 check_finite : bool, optional 

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

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

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

1009 

1010 Returns