Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-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 

8from itertools import product 

9import numpy as np 

10from numpy import atleast_1d, atleast_2d 

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 

16from ._cythonized_array_utils import find_det_from_lu 

17from scipy._lib.deprecation import _NoValue, _deprecate_positional_args 

18 

19# deprecated imports to be removed in SciPy 1.13.0 

20from scipy.linalg._flinalg_py import get_flinalg_funcs # noqa 

21 

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

23 'solve_toeplitz', 'solve_circulant', 'inv', 'det', 'lstsq', 

24 'pinv', 'pinvh', 'matrix_balance', 'matmul_toeplitz'] 

25 

26 

27# The numpy facilities for type-casting checks are too slow for small sized 

28# arrays and eat away the time budget for the checkups. Here we set a 

29# precomputed dict container of the numpy.can_cast() table. 

30 

31# It can be used to determine quickly what a dtype can be cast to LAPACK 

32# compatible types, i.e., 'float32, float64, complex64, complex128'. 

33# Then it can be checked via "casting_dict[arr.dtype.char]" 

34lapack_cast_dict = {x: ''.join([y for y in 'fdFD' if np.can_cast(x, y)]) 

35 for x in np.typecodes['All']} 

36 

37 

38# Linear equations 

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

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

41 if info < 0: 

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

43 '.'.format(-info)) 

44 elif 0 < info: 

45 raise LinAlgError('Matrix is singular.') 

46 

47 if lamch is None: 

48 return 

49 E = lamch('E') 

50 if rcond < E: 

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

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

53 LinAlgWarning, stacklevel=3) 

54 

55 

56def solve(a, b, lower=False, overwrite_a=False, 

57 overwrite_b=False, check_finite=True, assume_a='gen', 

58 transposed=False): 

59 """ 

60 Solves the linear equation set ``a @ x == b`` for the unknown ``x`` 

61 for square `a` matrix. 

62 

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

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

65 The available options are 

66 

67 =================== ======== 

68 generic matrix 'gen' 

69 symmetric 'sym' 

70 hermitian 'her' 

71 positive definite 'pos' 

72 =================== ======== 

73 

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

75 

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

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

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

79 on the data type of the array. 

80 

81 Parameters 

82 ---------- 

83 a : (N, N) array_like 

84 Square input data 

85 b : (N, NRHS) array_like 

86 Input data for the right hand side. 

87 lower : bool, default: False 

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

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

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

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

92 below the diagonal are ignored. 

93 overwrite_a : bool, default: False 

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

95 overwrite_b : bool, default: False 

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

97 check_finite : bool, default: True 

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

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

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

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

102 Valid entries are explained above. 

103 transposed : bool, default: False 

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

105 for complex `a`. 

106 

107 Returns 

108 ------- 

109 x : (N, NRHS) ndarray 

110 The solution array. 

111 

112 Raises 

113 ------ 

114 ValueError 

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

116 LinAlgError 

117 If the matrix is singular. 

118 LinAlgWarning 

119 If an ill-conditioned input a is detected. 

120 NotImplementedError 

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

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 Examples 

135 -------- 

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

137 

138 >>> import numpy as np 

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

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

141 >>> from scipy import linalg 

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

143 >>> x 

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

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

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

147 

148 """ 

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

150 b_is_1D = False 

151 

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

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

154 n = a1.shape[0] 

155 

156 overwrite_a = overwrite_a or _datacopied(a1, a) 

157 overwrite_b = overwrite_b or _datacopied(b1, b) 

158 

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

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

161 

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

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

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

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

166 'input a') 

167 

168 # accommodate empty arrays 

169 if b1.size == 0: 

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

171 

172 # regularize 1-D b arrays to 2D 

173 if b1.ndim == 1: 

174 if n == 1: 

175 b1 = b1[None, :] 

176 else: 

177 b1 = b1[:, None] 

178 b_is_1D = True 

179 

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

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

182 ''.format(assume_a)) 

183 

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

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

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

187 assume_a = 'sym' 

188 

189 # Get the correct lamch function. 

190 # The LAMCH functions only exists for S and D 

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

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

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

194 else: 

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

196 

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

198 # lansy, lanpo, lanhe. 

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

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

201 

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

203 # we can collect them all in this one call 

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

205 # the I-norm should be used 

206 if transposed: 

207 trans = 1 

208 norm = 'I' 

209 if np.iscomplexobj(a1): 

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

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

212 'for complex matrices.') 

213 else: 

214 trans = 0 

215 norm = '1' 

216 

217 anorm = lange(norm, a1) 

218 

219 # Generalized case 'gesv' 

220 if assume_a == 'gen': 

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

222 (a1, b1)) 

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

224 _solve_check(n, info) 

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

226 trans=trans, overwrite_b=overwrite_b) 

227 _solve_check(n, info) 

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

229 # Hermitian case 'hesv' 

230 elif assume_a == 'her': 

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

232 'hesv_lwork'), (a1, b1)) 

233 lwork = _compute_lwork(hesv_lw, n, lower) 

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

235 lower=lower, 

236 overwrite_a=overwrite_a, 

237 overwrite_b=overwrite_b) 

238 _solve_check(n, info) 

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

240 # Symmetric case 'sysv' 

241 elif assume_a == 'sym': 

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

243 'sysv_lwork'), (a1, b1)) 

244 lwork = _compute_lwork(sysv_lw, n, lower) 

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

246 lower=lower, 

247 overwrite_a=overwrite_a, 

248 overwrite_b=overwrite_b) 

249 _solve_check(n, info) 

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

251 # Positive definite case 'posv' 

252 else: 

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

254 (a1, b1)) 

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

256 overwrite_a=overwrite_a, 

257 overwrite_b=overwrite_b) 

258 _solve_check(n, info) 

259 rcond, info = pocon(lu, anorm) 

260 

261 _solve_check(n, info, lamch, rcond) 

262 

263 if b_is_1D: 

264 x = x.ravel() 

265 

266 return x 

267 

268 

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

270 overwrite_b=False, check_finite=True): 

271 """ 

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

273 

274 Parameters 

275 ---------- 

276 a : (M, M) array_like 

277 A triangular matrix 

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

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

280 lower : bool, optional 

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

282 Default is to use upper triangle. 

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

284 Type of system to solve: 

285 

286 ======== ========= 

287 trans system 

288 ======== ========= 

289 0 or 'N' a x = b 

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

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

292 ======== ========= 

293 unit_diagonal : bool, optional 

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

295 will not be referenced. 

296 overwrite_b : bool, optional 

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

298 check_finite : bool, optional 

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

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

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

302 

303 Returns 

304 ------- 

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

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

307 

308 Raises 

309 ------ 

310 LinAlgError 

311 If `a` is singular 

312 

313 Notes 

314 ----- 

315 .. versionadded:: 0.9.0 

316 

317 Examples 

318 -------- 

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

320 

321 [3 0 0 0] [4] 

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

323 [1 0 1 0] [4] 

324 [1 1 1 1] [2] 

325 

326 >>> import numpy as np 

327 >>> from scipy.linalg import solve_triangular 

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

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

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

331 >>> x 

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

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

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

335 

336 """ 

337 

338 a1 = _asarray_validated(a, check_finite=check_finite) 

339 b1 = _asarray_validated(b, check_finite=check_finite) 

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

341 raise ValueError('expected square matrix') 

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

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

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

345 overwrite_b = overwrite_b or _datacopied(b1, b) 

346 

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

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

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

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

351 trans=trans, unitdiag=unit_diagonal) 

352 else: 

353 # transposed system is solved since trtrs expects Fortran ordering 

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

355 trans=not trans, unitdiag=unit_diagonal) 

356 

357 if info == 0: 

358 return x 

359 if info > 0: 

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

361 (info-1)) 

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

363 (-info)) 

364 

365 

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

367 check_finite=True): 

368 """ 

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

370 

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

372 

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

374 

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

376 

377 * a01 a12 a23 a34 a45 

378 a00 a11 a22 a33 a44 a55 

379 a10 a21 a32 a43 a54 * 

380 a20 a31 a42 a53 * * 

381 

382 Parameters 

383 ---------- 

384 (l, u) : (integer, integer) 

385 Number of non-zero lower and upper diagonals 

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

387 Banded matrix 

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

389 Right-hand side 

390 overwrite_ab : bool, optional 

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

392 overwrite_b : bool, optional 

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

394 check_finite : bool, optional 

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

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

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

398 

399 Returns 

400 ------- 

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

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

403 shape of `b`. 

404 

405 Examples 

406 -------- 

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

408 

409 [5 2 -1 0 0] [0] 

410 [1 4 2 -1 0] [1] 

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

412 [0 0 1 2 2] [2] 

413 [0 0 0 1 1] [3] 

414 

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

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

417 

418 [* * -1 -1 -1] 

419 ab = [* 2 2 2 2] 

420 [5 4 3 2 1] 

421 [1 1 1 1 *] 

422 

423 >>> import numpy as np 

424 >>> from scipy.linalg import solve_banded 

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

426 ... [0, 2, 2, 2, 2], 

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

428 ... [1, 1, 1, 1, 0]]) 

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

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

431 >>> x 

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

433 

434 """ 

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 Uses Thomas' Algorithm, which is more efficient than standard LU 

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

481 matrices. 

482 

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

484 diagonal ordered form: 

485 

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

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

488 

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

490 ``u`` =2):: 

491 

492 upper form: 

493 * * a02 a13 a24 a35 

494 * a01 a12 a23 a34 a45 

495 a00 a11 a22 a33 a44 a55 

496 

497 lower form: 

498 a00 a11 a22 a33 a44 a55 

499 a10 a21 a32 a43 a54 * 

500 a20 a31 a42 a53 * * 

501 

502 Cells marked with * are not used. 

503 

504 Parameters 

505 ---------- 

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

507 Banded matrix 

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

509 Right-hand side 

510 overwrite_ab : bool, optional 

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

512 overwrite_b : bool, optional 

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

514 lower : bool, optional 

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

516 check_finite : bool, optional 

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

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

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

520 

521 Returns 

522 ------- 

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

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

525 of `b`. 

526 

527 Notes 

528 ----- 

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

530 `solve_banded` may be used. 

531 

532 Examples 

533 -------- 

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

535 

536 [ 4 2 -1 0 0 0] [1] 

537 [ 2 5 2 -1 0 0] [2] 

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

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

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

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

542 

543 >>> import numpy as np 

544 >>> from scipy.linalg import solveh_banded 

545 

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

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

548 

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

550 ... [ 2, 2, 2, 2, 2, 0], 

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

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

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

554 >>> x 

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

556 0.34733894]) 

557 

558 

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

560 

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

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

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

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

565 

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

567 

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

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

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

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

572 >>> x 

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

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

575 

576 """ 

577 a1 = _asarray_validated(ab, check_finite=check_finite) 

578 b1 = _asarray_validated(b, check_finite=check_finite) 

579 # Validate shapes. 

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

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

582 

583 overwrite_b = overwrite_b or _datacopied(b1, b) 

584 overwrite_ab = overwrite_ab or _datacopied(a1, ab) 

585 

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

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

588 if lower: 

589 d = a1[0, :].real 

590 e = a1[1, :-1] 

591 else: 

592 d = a1[1, :].real 

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

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

595 overwrite_b) 

596 else: 

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

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

599 overwrite_b=overwrite_b) 

600 if info > 0: 

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

602 if info < 0: 

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

604 'pbsv' % -info) 

605 return x 

606 

607 

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

609 """Solve a Toeplitz system using Levinson Recursion 

610 

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

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

613 assumed. 

614 

615 Parameters 

616 ---------- 

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

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

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

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

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

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

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

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

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

626 check_finite : bool, optional 

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

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

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

630 

631 Returns 

632 ------- 

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

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

635 of `b`. 

636 

637 See Also 

638 -------- 

639 toeplitz : Toeplitz matrix 

640 

641 Notes 

642 ----- 

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

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

645 

646 Examples 

647 -------- 

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

649 

650 [ 1 -1 -2 -3] [1] 

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

652 [ 6 3 1 -1] [2] 

653 [10 6 3 1] [5] 

654 

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

656 row are needed. 

657 

658 >>> import numpy as np 

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

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

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

662 

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

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

665 >>> x 

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

667 

668 Check the result by creating the full Toeplitz matrix and 

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

670 

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

672 >>> T.dot(x) 

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

674 

675 """ 

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

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

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

679 

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

681 c_or_cr, b, check_finite, keep_b_shape=True) 

682 

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

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

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

686 if b is None: 

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

688 

689 if b.ndim == 1: 

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

691 else: 

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

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

694 x = x.reshape(*b_shape) 

695 

696 return x 

697 

698 

699def _get_axis_len(aname, a, axis): 

700 ax = axis 

701 if ax < 0: 

702 ax += a.ndim 

703 if 0 <= ax < a.ndim: 

704 return a.shape[ax] 

705 raise ValueError(f"'{aname}axis' entry is out of bounds") 

706 

707 

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

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

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

711 

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

713 

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

715 calculation is:: 

716 

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

718 

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

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

721 solving the system with the full circulant matrix. 

722 

723 Parameters 

724 ---------- 

725 c : array_like 

726 The coefficients of the circulant matrix. 

727 b : array_like 

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

729 singular : str, optional 

730 This argument controls how a near singular circulant matrix is 

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

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

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

734 tol : float, optional 

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

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

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

738 

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

740 

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

742 of the circulant matrix. 

743 caxis : int 

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

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

746 holds the vectors of circulant coefficients. 

747 baxis : int 

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

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

750 right-hand side vectors. 

751 outaxis : int 

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

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

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

755 

756 Returns 

757 ------- 

758 x : ndarray 

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

760 

761 Raises 

762 ------ 

763 LinAlgError 

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

765 

766 See Also 

767 -------- 

768 circulant : circulant matrix 

769 

770 Notes 

771 ----- 

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

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

774 

775 solve_circulant(c, b) 

776 

777 returns the same result as 

778 

779 solve(circulant(c), b) 

780 

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

782 

783 .. versionadded:: 0.16.0 

784 

785 Examples 

786 -------- 

787 >>> import numpy as np 

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

789 

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

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

792 >>> solve_circulant(c, b) 

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

794 

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

796 

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

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

799 

800 A singular example: 

801 

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

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

804 

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

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

807 

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

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

810 

811 Compare to `scipy.linalg.lstsq`: 

812 

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

814 >>> x 

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

816 

817 A broadcasting example: 

818 

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

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

821 (3, 5). For example, 

822 

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

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

825 

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

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

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

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

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

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

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

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

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

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

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

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

838 

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

840 >>> x.shape 

841 (2, 3, 5) 

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

843 >>> x 

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

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

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

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

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

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

850 

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

852 

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

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

855 

856 """ 

857 c = np.atleast_1d(c) 

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

859 b = np.atleast_1d(b) 

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

861 if nc != nb: 

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

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

864 

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

866 abs_fc = np.abs(fc) 

867 if tol is None: 

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

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

870 if tol.shape != (): 

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

872 else: 

873 tol = np.atleast_1d(tol) 

874 

875 near_zeros = abs_fc <= tol 

876 is_near_singular = np.any(near_zeros) 

877 if is_near_singular: 

878 if singular == 'raise': 

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

880 else: 

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

882 # division fb/fc below. 

883 fc[near_zeros] = 1 

884 

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

886 

887 q = fb / fc 

888 

889 if is_near_singular: 

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

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

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

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

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

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

896 q[mask] = 0 

897 

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

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

900 x = x.real 

901 if outaxis != -1: 

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

903 return x 

904 

905 

906# matrix inversion 

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

908 """ 

909 Compute the inverse of a matrix. 

910 

911 Parameters 

912 ---------- 

913 a : array_like 

914 Square matrix to be inverted. 

915 overwrite_a : bool, optional 

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

917 check_finite : bool, optional 

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

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

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

921 

922 Returns 

923 ------- 

924 ainv : ndarray 

925 Inverse of the matrix `a`. 

926 

927 Raises 

928 ------ 

929 LinAlgError 

930 If `a` is singular. 

931 ValueError 

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

933 

934 Examples 

935 -------- 

936 >>> import numpy as np 

937 >>> from scipy import linalg 

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

939 >>> linalg.inv(a) 

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

941 [ 1.5, -0.5]]) 

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

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

944 [ 0., 1.]]) 

945 

946 """ 

947 a1 = _asarray_validated(a, check_finite=check_finite) 

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

949 raise ValueError('expected square matrix') 

950 overwrite_a = overwrite_a or _datacopied(a1, a) 

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

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

953# if finv is not None: 

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

955# if info==0: 

956# return a_inv 

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

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

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

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

961 'getri_lwork'), 

962 (a1,)) 

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

964 if info == 0: 

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

966 

967 # XXX: the following line fixes curious SEGFAULT when 

968 # benchmarking 500x500 matrix inverse. This seems to 

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

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

971 # all tests pass. Further investigation is required if 

972 # more such SEGFAULTs occur. 

973 lwork = int(1.01 * lwork) 

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

975 if info > 0: 

976 raise LinAlgError("singular matrix") 

977 if info < 0: 

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

979 'getrf|getri' % -info) 

980 return inv_a 

981 

982 

983# Determinant 

984 

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

986 """ 

987 Compute the determinant of a matrix 

988 

989 The determinant is a scalar that is a function of the associated square 

990 matrix coefficients. The determinant value is zero for singular matrices. 

991 

992 Parameters 

993 ---------- 

994 a : (..., M, M) array_like 

995 Input array to compute determinants for. 

996 overwrite_a : bool, optional 

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

998 check_finite : bool, optional 

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

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

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

1002 

1003 Returns 

1004 ------- 

1005 det : (...) float or complex