Coverage for /usr/lib/python3/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 

17 

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

19 'solve_toeplitz', 'solve_circulant', 'inv', 'det', 'lstsq', 

20 'pinv', 'pinvh', 'matrix_balance', 'matmul_toeplitz'] 

21 

22 

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

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

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

26 

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

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

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

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

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

32 

33 

34# Linear equations 

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

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

37 if info < 0: 

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

39 '.'.format(-info)) 

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('Ill-conditioned matrix (rcond={:.6g}): ' 

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

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('{} is not a recognized matrix structure' 

178 ''.format(assume_a)) 

179 

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

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

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

183 assume_a = 'sym' 

184 

185 # Get the correct lamch function. 

186 # The LAMCH functions only exists for S and D 

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

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

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

190 else: 

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

192 

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

194 # lansy, lanpo, lanhe. 

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

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

197 

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

199 # we can collect them all in this one call 

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

201 # the I-norm should be used 

202 if transposed: 

203 trans = 1 

204 norm = 'I' 

205 if np.iscomplexobj(a1): 

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

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

208 'for complex matrices.') 

209 else: 

210 trans = 0 

211 norm = '1' 

212 

213 anorm = lange(norm, a1) 

214 

215 # Generalized case 'gesv' 

216 if assume_a == 'gen': 

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

218 (a1, b1)) 

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

220 _solve_check(n, info) 

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

222 trans=trans, overwrite_b=overwrite_b) 

223 _solve_check(n, info) 

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

225 # Hermitian case 'hesv' 

226 elif assume_a == 'her': 

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

228 'hesv_lwork'), (a1, b1)) 

229 lwork = _compute_lwork(hesv_lw, n, lower) 

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

231 lower=lower, 

232 overwrite_a=overwrite_a, 

233 overwrite_b=overwrite_b) 

234 _solve_check(n, info) 

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

236 # Symmetric case 'sysv' 

237 elif assume_a == 'sym': 

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

239 'sysv_lwork'), (a1, b1)) 

240 lwork = _compute_lwork(sysv_lw, n, lower) 

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

242 lower=lower, 

243 overwrite_a=overwrite_a, 

244 overwrite_b=overwrite_b) 

245 _solve_check(n, info) 

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

247 # Positive definite case 'posv' 

248 else: 

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

250 (a1, b1)) 

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

252 overwrite_a=overwrite_a, 

253 overwrite_b=overwrite_b) 

254 _solve_check(n, info) 

255 rcond, info = pocon(lu, anorm) 

256 

257 _solve_check(n, info, lamch, rcond) 

258 

259 if b_is_1D: 

260 x = x.ravel() 

261 

262 return x 

263 

264 

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

266 overwrite_b=False, check_finite=True): 

267 """ 

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

269 

270 Parameters 

271 ---------- 

272 a : (M, M) array_like 

273 A triangular matrix 

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

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

276 lower : bool, optional 

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

278 Default is to use upper triangle. 

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

280 Type of system to solve: 

281 

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

283 trans system 

284 ======== ========= 

285 0 or 'N' a x = b 

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

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

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

289 unit_diagonal : bool, optional 

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

291 will not be referenced. 

292 overwrite_b : bool, optional 

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

294 check_finite : bool, optional 

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

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

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

298 

299 Returns 

300 ------- 

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

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

303 

304 Raises 

305 ------ 

306 LinAlgError 

307 If `a` is singular 

308 

309 Notes 

310 ----- 

311 .. versionadded:: 0.9.0 

312 

313 Examples 

314 -------- 

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

316 

317 [3 0 0 0] [4] 

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

319 [1 0 1 0] [4] 

320 [1 1 1 1] [2] 

321 

322 >>> import numpy as np 

323 >>> from scipy.linalg import solve_triangular 

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

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

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

327 >>> x 

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

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

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

331 

332 """ 

333 

334 a1 = _asarray_validated(a, check_finite=check_finite) 

335 b1 = _asarray_validated(b, check_finite=check_finite) 

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

337 raise ValueError('expected square matrix') 

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

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

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

341 overwrite_b = overwrite_b or _datacopied(b1, b) 

342 

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

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

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

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

347 trans=trans, unitdiag=unit_diagonal) 

348 else: 

349 # transposed system is solved since trtrs expects Fortran ordering 

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

351 trans=not trans, unitdiag=unit_diagonal) 

352 

353 if info == 0: 

354 return x 

355 if info > 0: 

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

357 (info-1)) 

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

359 (-info)) 

360 

361 

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

363 check_finite=True): 

364 """ 

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

366 

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

368 

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

370 

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

372 

373 * a01 a12 a23 a34 a45 

374 a00 a11 a22 a33 a44 a55 

375 a10 a21 a32 a43 a54 * 

376 a20 a31 a42 a53 * * 

377 

378 Parameters 

379 ---------- 

380 (l, u) : (integer, integer) 

381 Number of non-zero lower and upper diagonals 

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

383 Banded matrix 

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

385 Right-hand side 

386 overwrite_ab : bool, optional 

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

388 overwrite_b : bool, optional 

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

390 check_finite : bool, optional 

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

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

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

394 

395 Returns 

396 ------- 

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

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

399 shape of `b`. 

400 

401 Examples 

402 -------- 

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

404 

405 [5 2 -1 0 0] [0] 

406 [1 4 2 -1 0] [1] 

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

408 [0 0 1 2 2] [2] 

409 [0 0 0 1 1] [3] 

410 

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

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

413 

414 [* * -1 -1 -1] 

415 ab = [* 2 2 2 2] 

416 [5 4 3 2 1] 

417 [1 1 1 1 *] 

418 

419 >>> import numpy as np 

420 >>> from scipy.linalg import solve_banded 

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

422 ... [0, 2, 2, 2, 2], 

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

424 ... [1, 1, 1, 1, 0]]) 

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

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

427 >>> x 

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

429 

430 """ 

431 

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

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

434 # Validate shapes. 

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

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

437 (nlower, nupper) = l_and_u 

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

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

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

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

442 

443 overwrite_b = overwrite_b or _datacopied(b1, b) 

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

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

446 b2 /= a1[1, 0] 

447 return b2 

448 if nlower == nupper == 1: 

449 overwrite_ab = overwrite_ab or _datacopied(a1, ab) 

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

451 du = a1[0, 1:] 

452 d = a1[1, :] 

453 dl = a1[2, :-1] 

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

455 overwrite_ab, overwrite_b) 

456 else: 

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

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

459 a2[nlower:, :] = a1 

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

461 overwrite_b=overwrite_b) 

462 if info == 0: 

463 return x 

464 if info > 0: 

465 raise LinAlgError("singular matrix") 

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

467 'gbsv/gtsv' % -info) 

468 

469 

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

471 check_finite=True): 

472 """ 

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

474 

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

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

477 matrices. 

478 

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

480 diagonal ordered form: 

481 

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

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

484 

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

486 ``u`` =2):: 

487 

488 upper form: 

489 * * a02 a13 a24 a35 

490 * a01 a12 a23 a34 a45 

491 a00 a11 a22 a33 a44 a55 

492 

493 lower form: 

494 a00 a11 a22 a33 a44 a55 

495 a10 a21 a32 a43 a54 * 

496 a20 a31 a42 a53 * * 

497 

498 Cells marked with * are not used. 

499 

500 Parameters 

501 ---------- 

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

503 Banded matrix 

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

505 Right-hand side 

506 overwrite_ab : bool, optional 

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

508 overwrite_b : bool, optional 

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

510 lower : bool, optional 

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

512 check_finite : bool, optional 

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

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

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

516 

517 Returns 

518 ------- 

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

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

521 of `b`. 

522 

523 Notes 

524 ----- 

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

526 `solve_banded` may be used. 

527 

528 Examples 

529 -------- 

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

531 

532 [ 4 2 -1 0 0 0] [1] 

533 [ 2 5 2 -1 0 0] [2] 

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

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

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

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

538 

539 >>> import numpy as np 

540 >>> from scipy.linalg import solveh_banded 

541 

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

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

544 

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

546 ... [ 2, 2, 2, 2, 2, 0], 

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

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

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

550 >>> x 

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

552 0.34733894]) 

553 

554 

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

556 

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

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

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

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

561 

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

563 

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

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

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

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

568 >>> x 

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

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

571 

572 """ 

573 a1 = _asarray_validated(ab, check_finite=check_finite) 

574 b1 = _asarray_validated(b, check_finite=check_finite) 

575 # Validate shapes. 

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

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

578 

579 overwrite_b = overwrite_b or _datacopied(b1, b) 

580 overwrite_ab = overwrite_ab or _datacopied(a1, ab) 

581 

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

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

584 if lower: 

585 d = a1[0, :].real 

586 e = a1[1, :-1] 

587 else: 

588 d = a1[1, :].real 

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

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

591 overwrite_b) 

592 else: 

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

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

595 overwrite_b=overwrite_b) 

596 if info > 0: 

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

598 if info < 0: 

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

600 'pbsv' % -info) 

601 return x 

602 

603 

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

605 """Solve a Toeplitz system using Levinson Recursion 

606 

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

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

609 assumed. 

610 

611 Parameters 

612 ---------- 

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

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

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

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

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

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

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

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

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

622 check_finite : bool, optional 

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

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

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

626 

627 Returns 

628 ------- 

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

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

631 of `b`. 

632 

633 See Also 

634 -------- 

635 toeplitz : Toeplitz matrix 

636 

637 Notes 

638 ----- 

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

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

641 

642 Examples 

643 -------- 

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

645 

646 [ 1 -1 -2 -3] [1] 

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

648 [ 6 3 1 -1] [2] 

649 [10 6 3 1] [5] 

650 

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

652 row are needed. 

653 

654 >>> import numpy as np 

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

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

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

658 

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

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

661 >>> x 

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

663 

664 Check the result by creating the full Toeplitz matrix and 

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

666 

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

668 >>> T.dot(x) 

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

670 

671 """ 

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

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

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

675 

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

677 c_or_cr, b, check_finite, keep_b_shape=True) 

678 

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

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

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

682 if b is None: 

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

684 

685 if b.ndim == 1: 

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

687 else: 

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

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

690 x = x.reshape(*b_shape) 

691 

692 return x 

693 

694 

695def _get_axis_len(aname, a, axis): 

696 ax = axis 

697 if ax < 0: 

698 ax += a.ndim 

699 if 0 <= ax < a.ndim: 

700 return a.shape[ax] 

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

702 

703 

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

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

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

707 

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

709 

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

711 calculation is:: 

712 

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

714 

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

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

717 solving the system with the full circulant matrix. 

718 

719 Parameters 

720 ---------- 

721 c : array_like 

722 The coefficients of the circulant matrix. 

723 b : array_like 

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

725 singular : str, optional 

726 This argument controls how a near singular circulant matrix is 

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

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

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

730 tol : float, optional 

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

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

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

734 

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

736 

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

738 of the circulant matrix. 

739 caxis : int 

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

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

742 holds the vectors of circulant coefficients. 

743 baxis : int 

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

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

746 right-hand side vectors. 

747 outaxis : int 

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

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

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

751 

752 Returns 

753 ------- 

754 x : ndarray 

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

756 

757 Raises 

758 ------ 

759 LinAlgError 

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

761 

762 See Also 

763 -------- 

764 circulant : circulant matrix 

765 

766 Notes 

767 ----- 

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

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

770 

771 solve_circulant(c, b) 

772 

773 returns the same result as 

774 

775 solve(circulant(c), b) 

776 

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

778 

779 .. versionadded:: 0.16.0 

780 

781 Examples 

782 -------- 

783 >>> import numpy as np 

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

785 

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

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

788 >>> solve_circulant(c, b) 

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

790 

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

792 

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

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

795 

796 A singular example: 

797 

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

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

800 

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

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

803 

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

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

806 

807 Compare to `scipy.linalg.lstsq`: 

808 

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

810 >>> x 

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

812 

813 A broadcasting example: 

814 

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

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

817 (3, 5). For example, 

818 

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

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

821 

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

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

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

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

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

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

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

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

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

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

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

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

834 

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

836 >>> x.shape 

837 (2, 3, 5) 

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

839 >>> x 

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

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

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

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

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

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

846 

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

848 

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

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

851 

852 """ 

853 c = np.atleast_1d(c) 

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

855 b = np.atleast_1d(b) 

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

857 if nc != nb: 

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

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

860 

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

862 abs_fc = np.abs(fc) 

863 if tol is None: 

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

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

866 if tol.shape != (): 

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

868 else: 

869 tol = np.atleast_1d(tol) 

870 

871 near_zeros = abs_fc <= tol 

872 is_near_singular = np.any(near_zeros) 

873 if is_near_singular: 

874 if singular == 'raise': 

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

876 else: 

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

878 # division fb/fc below. 

879 fc[near_zeros] = 1 

880 

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

882 

883 q = fb / fc 

884 

885 if is_near_singular: 

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

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

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

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

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

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

892 q[mask] = 0 

893 

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

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

896 x = x.real 

897 if outaxis != -1: 

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

899 return x 

900 

901 

902# matrix inversion 

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

904 """ 

905 Compute the inverse of a matrix. 

906 

907 Parameters 

908 ---------- 

909 a : array_like 

910 Square matrix to be inverted. 

911 overwrite_a : bool, optional 

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

913 check_finite : bool, optional 

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

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

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

917 

918 Returns 

919 ------- 

920 ainv : ndarray 

921 Inverse of the matrix `a`. 

922 

923 Raises 

924 ------ 

925 LinAlgError 

926 If `a` is singular. 

927 ValueError 

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

929 

930 Examples 

931 -------- 

932 >>> import numpy as np 

933 >>> from scipy import linalg 

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

935 >>> linalg.inv(a) 

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

937 [ 1.5, -0.5]]) 

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

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

940 [ 0., 1.]]) 

941 

942 """ 

943 a1 = _asarray_validated(a, check_finite=check_finite) 

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

945 raise ValueError('expected square matrix') 

946 overwrite_a = overwrite_a or _datacopied(a1, a) 

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

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

949# if finv is not None: 

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

951# if info==0: 

952# return a_inv 

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

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

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

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

957 'getri_lwork'), 

958 (a1,)) 

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

960 if info == 0: 

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

962 

963 # XXX: the following line fixes curious SEGFAULT when 

964 # benchmarking 500x500 matrix inverse. This seems to 

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

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

967 # all tests pass. Further investigation is required if 

968 # more such SEGFAULTs occur. 

969 lwork = int(1.01 * lwork) 

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

971 if info > 0: 

972 raise LinAlgError("singular matrix") 

973 if info < 0: 

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

975 'getrf|getri' % -info) 

976 return inv_a 

977 

978 

979# Determinant 

980 

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

982 """ 

983 Compute the determinant of a matrix 

984 

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

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

987 

988 Parameters 

989 ---------- 

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

991 Input array to compute determinants for. 

992 overwrite_a : bool, optional 

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

994 check_finite : bool, optional 

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

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

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

998 

999 Returns 

1000 ------- 

1001 det : (...) float or complex 

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

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

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

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

1006 

1007 Notes