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__all__ = ['solve', 'solve_triangular', 'solveh_banded', 'solve_banded', 

20 'solve_toeplitz', 'solve_circulant', 'inv', 'det', 'lstsq', 

21 'pinv', 'pinvh', 'matrix_balance', 'matmul_toeplitz'] 

22 

23 

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

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

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

27 

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

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

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

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

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

33 

34 

35# Linear equations 

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

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

38 if info < 0: 

39 raise ValueError(f'LAPACK reported an illegal value in {-info}-th argument.') 

40 elif 0 < info: 

41 raise LinAlgError('Matrix is singular.') 

42 

43 if lamch is None: 

44 return 

45 E = lamch('E') 

46 if rcond < E: 

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

48 'result may not be accurate.', 

49 LinAlgWarning, stacklevel=3) 

50 

51 

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

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

54 transposed=False): 

55 """ 

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

57 for square `a` matrix. 

58 

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

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

61 The available options are 

62 

63 =================== ======== 

64 generic matrix 'gen' 

65 symmetric 'sym' 

66 hermitian 'her' 

67 positive definite 'pos' 

68 =================== ======== 

69 

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

71 

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

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

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

75 on the data type of the array. 

76 

77 Parameters 

78 ---------- 

79 a : (N, N) array_like 

80 Square input data 

81 b : (N, NRHS) array_like 

82 Input data for the right hand side. 

83 lower : bool, default: False 

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

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

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

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

88 below the diagonal are ignored. 

89 overwrite_a : bool, default: False 

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

91 overwrite_b : bool, default: False 

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

93 check_finite : bool, default: True 

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

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

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

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

98 Valid entries are explained above. 

99 transposed : bool, default: False 

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

101 for complex `a`. 

102 

103 Returns 

104 ------- 

105 x : (N, NRHS) ndarray 

106 The solution array. 

107 

108 Raises 

109 ------ 

110 ValueError 

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

112 LinAlgError 

113 If the matrix is singular. 

114 LinAlgWarning 

115 If an ill-conditioned input a is detected. 

116 NotImplementedError 

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

118 

119 Notes 

120 ----- 

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

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

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

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

125 

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

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

128 LAPACK respectively. 

129 

130 Examples 

131 -------- 

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

133 

134 >>> import numpy as np 

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

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

137 >>> from scipy import linalg 

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

139 >>> x 

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

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

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

143 

144 """ 

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

146 b_is_1D = False 

147 

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

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

150 n = a1.shape[0] 

151 

152 overwrite_a = overwrite_a or _datacopied(a1, a) 

153 overwrite_b = overwrite_b or _datacopied(b1, b) 

154 

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

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

157 

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

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

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

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

162 'input a') 

163 

164 # accommodate empty arrays 

165 if b1.size == 0: 

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

167 

168 # regularize 1-D b arrays to 2D 

169 if b1.ndim == 1: 

170 if n == 1: 

171 b1 = b1[None, :] 

172 else: 

173 b1 = b1[:, None] 

174 b_is_1D = True 

175 

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

177 raise ValueError(f'{assume_a} is not a recognized matrix structure') 

178 

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

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

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

182 assume_a = 'sym' 

183 

184 # Get the correct lamch function. 

185 # The LAMCH functions only exists for S and D 

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

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

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

189 else: 

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

191 

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

193 # lansy, lanpo, lanhe. 

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

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

196 

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

198 # we can collect them all in this one call 

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

200 # the I-norm should be used 

201 if transposed: 

202 trans = 1 

203 norm = 'I' 

204 if np.iscomplexobj(a1): 

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

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

207 'for complex matrices.') 

208 else: 

209 trans = 0 

210 norm = '1' 

211 

212 anorm = lange(norm, a1) 

213 

214 # Generalized case 'gesv' 

215 if assume_a == 'gen': 

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

217 (a1, b1)) 

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

219 _solve_check(n, info) 

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

221 trans=trans, overwrite_b=overwrite_b) 

222 _solve_check(n, info) 

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

224 # Hermitian case 'hesv' 

225 elif assume_a == 'her': 

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

227 'hesv_lwork'), (a1, b1)) 

228 lwork = _compute_lwork(hesv_lw, n, lower) 

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

230 lower=lower, 

231 overwrite_a=overwrite_a, 

232 overwrite_b=overwrite_b) 

233 _solve_check(n, info) 

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

235 # Symmetric case 'sysv' 

236 elif assume_a == 'sym': 

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

238 'sysv_lwork'), (a1, b1)) 

239 lwork = _compute_lwork(sysv_lw, n, lower) 

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

241 lower=lower, 

242 overwrite_a=overwrite_a, 

243 overwrite_b=overwrite_b) 

244 _solve_check(n, info) 

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

246 # Positive definite case 'posv' 

247 else: 

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

249 (a1, b1)) 

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

251 overwrite_a=overwrite_a, 

252 overwrite_b=overwrite_b) 

253 _solve_check(n, info) 

254 rcond, info = pocon(lu, anorm) 

255 

256 _solve_check(n, info, lamch, rcond) 

257 

258 if b_is_1D: 

259 x = x.ravel() 

260 

261 return x 

262 

263 

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

265 overwrite_b=False, check_finite=True): 

266 """ 

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

268 

269 Parameters 

270 ---------- 

271 a : (M, M) array_like 

272 A triangular matrix 

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

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

275 lower : bool, optional 

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

277 Default is to use upper triangle. 

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

279 Type of system to solve: 

280 

281 ======== ========= 

282 trans system 

283 ======== ========= 

284 0 or 'N' a x = b 

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

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

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

288 unit_diagonal : bool, optional 

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

290 will not be referenced. 

291 overwrite_b : bool, optional 

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

293 check_finite : bool, optional 

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

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

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

297 

298 Returns 

299 ------- 

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

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

302 

303 Raises 

304 ------ 

305 LinAlgError 

306 If `a` is singular 

307 

308 Notes 

309 ----- 

310 .. versionadded:: 0.9.0 

311 

312 Examples 

313 -------- 

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

315 

316 [3 0 0 0] [4] 

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

318 [1 0 1 0] [4] 

319 [1 1 1 1] [2] 

320 

321 >>> import numpy as np 

322 >>> from scipy.linalg import solve_triangular 

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

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

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

326 >>> x 

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

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

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

330 

331 """ 

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(f'shapes of a {a1.shape} and b {b1.shape} are incompatible') 

339 overwrite_b = overwrite_b or _datacopied(b1, b) 

340 

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

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

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

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

345 trans=trans, unitdiag=unit_diagonal) 

346 else: 

347 # transposed system is solved since trtrs expects Fortran ordering 

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

349 trans=not trans, unitdiag=unit_diagonal) 

350 

351 if info == 0: 

352 return x 

353 if info > 0: 

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

355 (info-1)) 

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

357 (-info)) 

358 

359 

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

361 check_finite=True): 

362 """ 

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

364 

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

366 

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

368 

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

370 

371 * a01 a12 a23 a34 a45 

372 a00 a11 a22 a33 a44 a55 

373 a10 a21 a32 a43 a54 * 

374 a20 a31 a42 a53 * * 

375 

376 Parameters 

377 ---------- 

378 (l, u) : (integer, integer) 

379 Number of non-zero lower and upper diagonals 

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

381 Banded matrix 

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

383 Right-hand side 

384 overwrite_ab : bool, optional 

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

386 overwrite_b : bool, optional 

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

388 check_finite : bool, optional 

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

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

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

392 

393 Returns 

394 ------- 

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

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

397 shape of `b`. 

398 

399 Examples 

400 -------- 

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

402 

403 [5 2 -1 0 0] [0] 

404 [1 4 2 -1 0] [1] 

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

406 [0 0 1 2 2] [2] 

407 [0 0 0 1 1] [3] 

408 

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

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

411 

412 [* * -1 -1 -1] 

413 ab = [* 2 2 2 2] 

414 [5 4 3 2 1] 

415 [1 1 1 1 *] 

416 

417 >>> import numpy as np 

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 a1 = _asarray_validated(ab, check_finite=check_finite, as_inexact=True) 

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

432 # Validate shapes. 

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

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

435 (nlower, nupper) = l_and_u 

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

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

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

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

440 

441 overwrite_b = overwrite_b or _datacopied(b1, b) 

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

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

444 b2 /= a1[1, 0] 

445 return b2 

446 if nlower == nupper == 1: 

447 overwrite_ab = overwrite_ab or _datacopied(a1, ab) 

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

449 du = a1[0, 1:] 

450 d = a1[1, :] 

451 dl = a1[2, :-1] 

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

453 overwrite_ab, overwrite_b) 

454 else: 

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

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

457 a2[nlower:, :] = a1 

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

459 overwrite_b=overwrite_b) 

460 if info == 0: 

461 return x 

462 if info > 0: 

463 raise LinAlgError("singular matrix") 

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

465 'gbsv/gtsv' % -info) 

466 

467 

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

469 check_finite=True): 

470 """ 

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

472 

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

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

475 matrices. 

476 

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

478 diagonal ordered form: 

479 

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

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

482 

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

484 ``u`` =2):: 

485 

486 upper form: 

487 * * a02 a13 a24 a35 

488 * a01 a12 a23 a34 a45 

489 a00 a11 a22 a33 a44 a55 

490 

491 lower form: 

492 a00 a11 a22 a33 a44 a55 

493 a10 a21 a32 a43 a54 * 

494 a20 a31 a42 a53 * * 

495 

496 Cells marked with * are not used. 

497 

498 Parameters 

499 ---------- 

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

501 Banded matrix 

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

503 Right-hand side 

504 overwrite_ab : bool, optional 

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

506 overwrite_b : bool, optional 

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

508 lower : bool, optional 

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

510 check_finite : bool, optional 

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

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

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

514 

515 Returns 

516 ------- 

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

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

519 of `b`. 

520 

521 Notes 

522 ----- 

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

524 `solve_banded` may be used. 

525 

526 Examples 

527 -------- 

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

529 

530 [ 4 2 -1 0 0 0] [1] 

531 [ 2 5 2 -1 0 0] [2] 

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

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

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

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

536 

537 >>> import numpy as np 

538 >>> from scipy.linalg import solveh_banded 

539 

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

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

542 

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

544 ... [ 2, 2, 2, 2, 2, 0], 

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

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

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

548 >>> x 

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

550 0.34733894]) 

551 

552 

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

554 

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

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

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

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

559 

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

561 

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

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

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

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

566 >>> x 

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

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

569 

570 """ 

571 a1 = _asarray_validated(ab, check_finite=check_finite) 

572 b1 = _asarray_validated(b, check_finite=check_finite) 

573 # Validate shapes. 

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

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

576 

577 overwrite_b = overwrite_b or _datacopied(b1, b) 

578 overwrite_ab = overwrite_ab or _datacopied(a1, ab) 

579 

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

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

582 if lower: 

583 d = a1[0, :].real 

584 e = a1[1, :-1] 

585 else: 

586 d = a1[1, :].real 

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

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

589 overwrite_b) 

590 else: 

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

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

593 overwrite_b=overwrite_b) 

594 if info > 0: 

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

596 if info < 0: 

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

598 'pbsv' % -info) 

599 return x 

600 

601 

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

603 """Solve a Toeplitz system using Levinson Recursion 

604 

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

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

607 assumed. 

608 

609 Parameters 

610 ---------- 

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

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

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

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

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

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

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

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

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

620 check_finite : bool, optional 

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

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

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

624 

625 Returns 

626 ------- 

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

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

629 of `b`. 

630 

631 See Also 

632 -------- 

633 toeplitz : Toeplitz matrix 

634 

635 Notes 

636 ----- 

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

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

639 

640 Examples 

641 -------- 

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

643 

644 [ 1 -1 -2 -3] [1] 

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

646 [ 6 3 1 -1] [2] 

647 [10 6 3 1] [5] 

648 

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

650 row are needed. 

651 

652 >>> import numpy as np 

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

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

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

656 

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

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

659 >>> x 

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

661 

662 Check the result by creating the full Toeplitz matrix and 

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

664 

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

666 >>> T.dot(x) 

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

668 

669 """ 

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

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

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

673 

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

675 c_or_cr, b, check_finite, keep_b_shape=True) 

676 

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

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

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

680 if b is None: 

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

682 

683 if b.ndim == 1: 

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

685 else: 

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

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

688 x = x.reshape(*b_shape) 

689 

690 return x 

691 

692 

693def _get_axis_len(aname, a, axis): 

694 ax = axis 

695 if ax < 0: 

696 ax += a.ndim 

697 if 0 <= ax < a.ndim: 

698 return a.shape[ax] 

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

700 

701 

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

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

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

705 

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

707 

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

709 calculation is:: 

710 

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

712 

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

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

715 solving the system with the full circulant matrix. 

716 

717 Parameters 

718 ---------- 

719 c : array_like 

720 The coefficients of the circulant matrix. 

721 b : array_like 

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

723 singular : str, optional 

724 This argument controls how a near singular circulant matrix is 

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

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

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

728 tol : float, optional 

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

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

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

732 

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

734 

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

736 of the circulant matrix. 

737 caxis : int 

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

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

740 holds the vectors of circulant coefficients. 

741 baxis : int 

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

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

744 right-hand side vectors. 

745 outaxis : int 

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

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

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

749 

750 Returns 

751 ------- 

752 x : ndarray 

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

754 

755 Raises 

756 ------ 

757 LinAlgError 

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

759 

760 See Also 

761 -------- 

762 circulant : circulant matrix 

763 

764 Notes 

765 ----- 

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

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

768 

769 solve_circulant(c, b) 

770 

771 returns the same result as 

772 

773 solve(circulant(c), b) 

774 

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

776 

777 .. versionadded:: 0.16.0 

778 

779 Examples 

780 -------- 

781 >>> import numpy as np 

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

783 

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

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

786 >>> solve_circulant(c, b) 

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

788 

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

790 

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

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

793 

794 A singular example: 

795 

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

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

798 

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

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

801 

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

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

804 

805 Compare to `scipy.linalg.lstsq`: 

806 

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

808 >>> x 

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

810 

811 A broadcasting example: 

812 

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

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

815 (3, 5). For example, 

816 

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

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

819 

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

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

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

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

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

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

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

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

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

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

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

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

832 

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

834 >>> x.shape 

835 (2, 3, 5) 

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

837 >>> x 

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

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

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

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

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

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

844 

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

846 

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

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

849 

850 """ 

851 c = np.atleast_1d(c) 

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

853 b = np.atleast_1d(b) 

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

855 if nc != nb: 

856 raise ValueError(f'Shapes of c {c.shape} and b {b.shape} are incompatible') 

857 

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

859 abs_fc = np.abs(fc) 

860 if tol is None: 

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

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

863 if tol.shape != (): 

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

865 else: 

866 tol = np.atleast_1d(tol) 

867 

868 near_zeros = abs_fc <= tol 

869 is_near_singular = np.any(near_zeros) 

870 if is_near_singular: 

871 if singular == 'raise': 

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

873 else: 

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

875 # division fb/fc below. 

876 fc[near_zeros] = 1 

877 

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

879 

880 q = fb / fc 

881 

882 if is_near_singular: 

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

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

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

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

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

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

889 q[mask] = 0 

890 

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

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

893 x = x.real 

894 if outaxis != -1: 

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

896 return x 

897 

898 

899# matrix inversion 

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

901 """ 

902 Compute the inverse of a matrix. 

903 

904 Parameters 

905 ---------- 

906 a : array_like 

907 Square matrix to be inverted. 

908 overwrite_a : bool, optional 

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

910 check_finite : bool, optional 

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

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

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

914 

915 Returns 

916 ------- 

917 ainv : ndarray 

918 Inverse of the matrix `a`. 

919 

920 Raises 

921 ------ 

922 LinAlgError 

923 If `a` is singular. 

924 ValueError 

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

926 

927 Examples 

928 -------- 

929 >>> import numpy as np 

930 >>> from scipy import linalg 

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

932 >>> linalg.inv(a) 

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

934 [ 1.5, -0.5]]) 

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

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

937 [ 0., 1.]]) 

938 

939 """ 

940 a1 = _asarray_validated(a, check_finite=check_finite) 

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

942 raise ValueError('expected square matrix') 

943 overwrite_a = overwrite_a or _datacopied(a1, a) 

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

945 'getri_lwork'), 

946 (a1,)) 

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

948 if info == 0: 

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

950 

951 # XXX: the following line fixes curious SEGFAULT when 

952 # benchmarking 500x500 matrix inverse. This seems to 

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

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

955 # all tests pass. Further investigation is required if 

956 # more such SEGFAULTs occur. 

957 lwork = int(1.01 * lwork) 

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

959 if info > 0: 

960 raise LinAlgError("singular matrix") 

961 if info < 0: 

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

963 'getrf|getri' % -info) 

964 return inv_a 

965 

966 

967# Determinant 

968 

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

970 """ 

971 Compute the determinant of a matrix 

972 

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

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

975 

976 Parameters 

977 ---------- 

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

979 Input array to compute determinants for. 

980 overwrite_a : bool, optional 

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

982 check_finite : bool, optional 

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

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

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

986 

987 Returns 

988 ------- 

989 det : (...) float or complex 

990 Determinant of `a`. For stacked arrays, a scalar is returned for each 

991 (m, m) slice in the last two dimensions of the input. For example, an 

992 input of shape (p, q, m, m) will produce a result of shape (p, q). If 

993 all dimensions are 1 a scalar is returned regardless of ndim. 

994 

995 Notes 

996 ----- 

997 The determinant is computed by performing an LU factorization of the 

998 input with LAPACK routine 'getrf', and then calculating the product of 

999 diagonal entries of the U factor. 

1000 

1001 Even the input array is single precision (float32 or complex64), the result 

1002 will be returned in double precision (float64 or complex128) to prevent 

1003 overflows. 

1004 

1005 Examples 

1006 -------- 

1007 >>> import numpy as np 

1008 >>> from scipy import linalg 

1009 >>> a = np.array([[1,2,3], [4,5,6], [7,8,9]]) # A singular matrix 

1010 >>> linalg.det(a) 

1011 0.0 

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