Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/numpy/linalg/linalg.py: 16%

Shortcuts on this page

r m x   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

651 statements  

1"""Lite version of scipy.linalg. 

2 

3Notes 

4----- 

5This module is a lite version of the linalg.py module in SciPy which 

6contains high-level Python interface to the LAPACK library. The lite 

7version only accesses the following LAPACK functions: dgesv, zgesv, 

8dgeev, zgeev, dgesdd, zgesdd, dgelsd, zgelsd, dsyevd, zheevd, dgetrf, 

9zgetrf, dpotrf, zpotrf, dgeqrf, zgeqrf, zungqr, dorgqr. 

10""" 

11 

12__all__ = ['matrix_power', 'solve', 'tensorsolve', 'tensorinv', 'inv', 

13 'cholesky', 'eigvals', 'eigvalsh', 'pinv', 'slogdet', 'det', 

14 'svd', 'eig', 'eigh', 'lstsq', 'norm', 'qr', 'cond', 'matrix_rank', 

15 'LinAlgError', 'multi_dot'] 

16 

17import functools 

18import operator 

19import warnings 

20 

21from numpy.core import ( 

22 array, asarray, zeros, empty, empty_like, intc, single, double, 

23 csingle, cdouble, inexact, complexfloating, newaxis, all, Inf, dot, 

24 add, multiply, sqrt, sum, isfinite, 

25 finfo, errstate, geterrobj, moveaxis, amin, amax, product, abs, 

26 atleast_2d, intp, asanyarray, object_, matmul, 

27 swapaxes, divide, count_nonzero, isnan, sign, argsort, sort, 

28 reciprocal 

29) 

30from numpy.core.multiarray import normalize_axis_index 

31from numpy.core.overrides import set_module 

32from numpy.core import overrides 

33from numpy.lib.twodim_base import triu, eye 

34from numpy.linalg import _umath_linalg 

35 

36 

37array_function_dispatch = functools.partial( 

38 overrides.array_function_dispatch, module='numpy.linalg') 

39 

40 

41fortran_int = intc 

42 

43 

44@set_module('numpy.linalg') 

45class LinAlgError(Exception): 

46 """ 

47 Generic Python-exception-derived object raised by linalg functions. 

48 

49 General purpose exception class, derived from Python's exception.Exception 

50 class, programmatically raised in linalg functions when a Linear 

51 Algebra-related condition would prevent further correct execution of the 

52 function. 

53 

54 Parameters 

55 ---------- 

56 None 

57 

58 Examples 

59 -------- 

60 >>> from numpy import linalg as LA 

61 >>> LA.inv(np.zeros((2,2))) 

62 Traceback (most recent call last): 

63 File "<stdin>", line 1, in <module> 

64 File "...linalg.py", line 350, 

65 in inv return wrap(solve(a, identity(a.shape[0], dtype=a.dtype))) 

66 File "...linalg.py", line 249, 

67 in solve 

68 raise LinAlgError('Singular matrix') 

69 numpy.linalg.LinAlgError: Singular matrix 

70 

71 """ 

72 

73 

74def _determine_error_states(): 

75 errobj = geterrobj() 

76 bufsize = errobj[0] 

77 

78 with errstate(invalid='call', over='ignore', 

79 divide='ignore', under='ignore'): 

80 invalid_call_errmask = geterrobj()[1] 

81 

82 return [bufsize, invalid_call_errmask, None] 

83 

84# Dealing with errors in _umath_linalg 

85_linalg_error_extobj = _determine_error_states() 

86del _determine_error_states 

87 

88def _raise_linalgerror_singular(err, flag): 

89 raise LinAlgError("Singular matrix") 

90 

91def _raise_linalgerror_nonposdef(err, flag): 

92 raise LinAlgError("Matrix is not positive definite") 

93 

94def _raise_linalgerror_eigenvalues_nonconvergence(err, flag): 

95 raise LinAlgError("Eigenvalues did not converge") 

96 

97def _raise_linalgerror_svd_nonconvergence(err, flag): 

98 raise LinAlgError("SVD did not converge") 

99 

100def _raise_linalgerror_lstsq(err, flag): 

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

102 

103def _raise_linalgerror_qr(err, flag): 

104 raise LinAlgError("Incorrect argument found while performing " 

105 "QR factorization") 

106 

107def get_linalg_error_extobj(callback): 

108 extobj = list(_linalg_error_extobj) # make a copy 

109 extobj[2] = callback 

110 return extobj 

111 

112def _makearray(a): 

113 new = asarray(a) 

114 wrap = getattr(a, "__array_prepare__", new.__array_wrap__) 

115 return new, wrap 

116 

117def isComplexType(t): 

118 return issubclass(t, complexfloating) 

119 

120_real_types_map = {single : single, 

121 double : double, 

122 csingle : single, 

123 cdouble : double} 

124 

125_complex_types_map = {single : csingle, 

126 double : cdouble, 

127 csingle : csingle, 

128 cdouble : cdouble} 

129 

130def _realType(t, default=double): 

131 return _real_types_map.get(t, default) 

132 

133def _complexType(t, default=cdouble): 

134 return _complex_types_map.get(t, default) 

135 

136def _commonType(*arrays): 

137 # in lite version, use higher precision (always double or cdouble) 

138 result_type = single 

139 is_complex = False 

140 for a in arrays: 

141 if issubclass(a.dtype.type, inexact): 

142 if isComplexType(a.dtype.type): 

143 is_complex = True 

144 rt = _realType(a.dtype.type, default=None) 

145 if rt is None: 

146 # unsupported inexact scalar 

147 raise TypeError("array type %s is unsupported in linalg" % 

148 (a.dtype.name,)) 

149 else: 

150 rt = double 

151 if rt is double: 

152 result_type = double 

153 if is_complex: 

154 t = cdouble 

155 result_type = _complex_types_map[result_type] 

156 else: 

157 t = double 

158 return t, result_type 

159 

160 

161def _to_native_byte_order(*arrays): 

162 ret = [] 

163 for arr in arrays: 

164 if arr.dtype.byteorder not in ('=', '|'): 

165 ret.append(asarray(arr, dtype=arr.dtype.newbyteorder('='))) 

166 else: 

167 ret.append(arr) 

168 if len(ret) == 1: 

169 return ret[0] 

170 else: 

171 return ret 

172 

173 

174def _assert_2d(*arrays): 

175 for a in arrays: 

176 if a.ndim != 2: 

177 raise LinAlgError('%d-dimensional array given. Array must be ' 

178 'two-dimensional' % a.ndim) 

179 

180def _assert_stacked_2d(*arrays): 

181 for a in arrays: 

182 if a.ndim < 2: 

183 raise LinAlgError('%d-dimensional array given. Array must be ' 

184 'at least two-dimensional' % a.ndim) 

185 

186def _assert_stacked_square(*arrays): 

187 for a in arrays: 

188 m, n = a.shape[-2:] 

189 if m != n: 

190 raise LinAlgError('Last 2 dimensions of the array must be square') 

191 

192def _assert_finite(*arrays): 

193 for a in arrays: 

194 if not isfinite(a).all(): 

195 raise LinAlgError("Array must not contain infs or NaNs") 

196 

197def _is_empty_2d(arr): 

198 # check size first for efficiency 

199 return arr.size == 0 and product(arr.shape[-2:]) == 0 

200 

201 

202def transpose(a): 

203 """ 

204 Transpose each matrix in a stack of matrices. 

205 

206 Unlike np.transpose, this only swaps the last two axes, rather than all of 

207 them 

208 

209 Parameters 

210 ---------- 

211 a : (...,M,N) array_like 

212 

213 Returns 

214 ------- 

215 aT : (...,N,M) ndarray 

216 """ 

217 return swapaxes(a, -1, -2) 

218 

219# Linear equations 

220 

221def _tensorsolve_dispatcher(a, b, axes=None): 

222 return (a, b) 

223 

224 

225@array_function_dispatch(_tensorsolve_dispatcher) 

226def tensorsolve(a, b, axes=None): 

227 """ 

228 Solve the tensor equation ``a x = b`` for x. 

229 

230 It is assumed that all indices of `x` are summed over in the product, 

231 together with the rightmost indices of `a`, as is done in, for example, 

232 ``tensordot(a, x, axes=x.ndim)``. 

233 

234 Parameters 

235 ---------- 

236 a : array_like 

237 Coefficient tensor, of shape ``b.shape + Q``. `Q`, a tuple, equals 

238 the shape of that sub-tensor of `a` consisting of the appropriate 

239 number of its rightmost indices, and must be such that 

240 ``prod(Q) == prod(b.shape)`` (in which sense `a` is said to be 

241 'square'). 

242 b : array_like 

243 Right-hand tensor, which can be of any shape. 

244 axes : tuple of ints, optional 

245 Axes in `a` to reorder to the right, before inversion. 

246 If None (default), no reordering is done. 

247 

248 Returns 

249 ------- 

250 x : ndarray, shape Q 

251 

252 Raises 

253 ------ 

254 LinAlgError 

255 If `a` is singular or not 'square' (in the above sense). 

256 

257 See Also 

258 -------- 

259 numpy.tensordot, tensorinv, numpy.einsum 

260 

261 Examples 

262 -------- 

263 >>> a = np.eye(2*3*4) 

264 >>> a.shape = (2*3, 4, 2, 3, 4) 

265 >>> b = np.random.randn(2*3, 4) 

266 >>> x = np.linalg.tensorsolve(a, b) 

267 >>> x.shape 

268 (2, 3, 4) 

269 >>> np.allclose(np.tensordot(a, x, axes=3), b) 

270 True 

271 

272 """ 

273 a, wrap = _makearray(a) 

274 b = asarray(b) 

275 an = a.ndim 

276 

277 if axes is not None: 

278 allaxes = list(range(0, an)) 

279 for k in axes: 

280 allaxes.remove(k) 

281 allaxes.insert(an, k) 

282 a = a.transpose(allaxes) 

283 

284 oldshape = a.shape[-(an-b.ndim):] 

285 prod = 1 

286 for k in oldshape: 

287 prod *= k 

288 

289 if a.size != prod ** 2: 

290 raise LinAlgError( 

291 "Input arrays must satisfy the requirement \ 

292 prod(a.shape[b.ndim:]) == prod(a.shape[:b.ndim])" 

293 ) 

294 

295 a = a.reshape(prod, prod) 

296 b = b.ravel() 

297 res = wrap(solve(a, b)) 

298 res.shape = oldshape 

299 return res 

300 

301 

302def _solve_dispatcher(a, b): 

303 return (a, b) 

304 

305 

306@array_function_dispatch(_solve_dispatcher) 

307def solve(a, b): 

308 """ 

309 Solve a linear matrix equation, or system of linear scalar equations. 

310 

311 Computes the "exact" solution, `x`, of the well-determined, i.e., full 

312 rank, linear matrix equation `ax = b`. 

313 

314 Parameters 

315 ---------- 

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

317 Coefficient matrix. 

318 b : {(..., M,), (..., M, K)}, array_like 

319 Ordinate or "dependent variable" values. 

320 

321 Returns 

322 ------- 

323 x : {(..., M,), (..., M, K)} ndarray 

324 Solution to the system a x = b. Returned shape is identical to `b`. 

325 

326 Raises 

327 ------ 

328 LinAlgError 

329 If `a` is singular or not square. 

330 

331 See Also 

332 -------- 

333 scipy.linalg.solve : Similar function in SciPy. 

334 

335 Notes 

336 ----- 

337 

338 .. versionadded:: 1.8.0 

339 

340 Broadcasting rules apply, see the `numpy.linalg` documentation for 

341 details. 

342 

343 The solutions are computed using LAPACK routine ``_gesv``. 

344 

345 `a` must be square and of full-rank, i.e., all rows (or, equivalently, 

346 columns) must be linearly independent; if either is not true, use 

347 `lstsq` for the least-squares best "solution" of the 

348 system/equation. 

349 

350 References 

351 ---------- 

352 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, 

353 FL, Academic Press, Inc., 1980, pg. 22. 

354 

355 Examples 

356 -------- 

357 Solve the system of equations ``x0 + 2 * x1 = 1`` and ``3 * x0 + 5 * x1 = 2``: 

358 

359 >>> a = np.array([[1, 2], [3, 5]]) 

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

361 >>> x = np.linalg.solve(a, b) 

362 >>> x 

363 array([-1., 1.]) 

364 

365 Check that the solution is correct: 

366 

367 >>> np.allclose(np.dot(a, x), b) 

368 True 

369 

370 """ 

371 a, _ = _makearray(a) 

372 _assert_stacked_2d(a) 

373 _assert_stacked_square(a) 

374 b, wrap = _makearray(b) 

375 t, result_t = _commonType(a, b) 

376 

377 # We use the b = (..., M,) logic, only if the number of extra dimensions 

378 # match exactly 

379 if b.ndim == a.ndim - 1: 

380 gufunc = _umath_linalg.solve1 

381 else: 

382 gufunc = _umath_linalg.solve 

383 

384 signature = 'DD->D' if isComplexType(t) else 'dd->d' 

385 extobj = get_linalg_error_extobj(_raise_linalgerror_singular) 

386 r = gufunc(a, b, signature=signature, extobj=extobj) 

387 

388 return wrap(r.astype(result_t, copy=False)) 

389 

390 

391def _tensorinv_dispatcher(a, ind=None): 

392 return (a,) 

393 

394 

395@array_function_dispatch(_tensorinv_dispatcher) 

396def tensorinv(a, ind=2): 

397 """ 

398 Compute the 'inverse' of an N-dimensional array. 

399 

400 The result is an inverse for `a` relative to the tensordot operation 

401 ``tensordot(a, b, ind)``, i. e., up to floating-point accuracy, 

402 ``tensordot(tensorinv(a), a, ind)`` is the "identity" tensor for the 

403 tensordot operation. 

404 

405 Parameters 

406 ---------- 

407 a : array_like 

408 Tensor to 'invert'. Its shape must be 'square', i. e., 

409 ``prod(a.shape[:ind]) == prod(a.shape[ind:])``. 

410 ind : int, optional 

411 Number of first indices that are involved in the inverse sum. 

412 Must be a positive integer, default is 2. 

413 

414 Returns 

415 ------- 

416 b : ndarray 

417 `a`'s tensordot inverse, shape ``a.shape[ind:] + a.shape[:ind]``. 

418 

419 Raises 

420 ------ 

421 LinAlgError 

422 If `a` is singular or not 'square' (in the above sense). 

423 

424 See Also 

425 -------- 

426 numpy.tensordot, tensorsolve 

427 

428 Examples 

429 -------- 

430 >>> a = np.eye(4*6) 

431 >>> a.shape = (4, 6, 8, 3) 

432 >>> ainv = np.linalg.tensorinv(a, ind=2) 

433 >>> ainv.shape 

434 (8, 3, 4, 6) 

435 >>> b = np.random.randn(4, 6) 

436 >>> np.allclose(np.tensordot(ainv, b), np.linalg.tensorsolve(a, b)) 

437 True 

438 

439 >>> a = np.eye(4*6) 

440 >>> a.shape = (24, 8, 3) 

441 >>> ainv = np.linalg.tensorinv(a, ind=1) 

442 >>> ainv.shape 

443 (8, 3, 24) 

444 >>> b = np.random.randn(24) 

445 >>> np.allclose(np.tensordot(ainv, b, 1), np.linalg.tensorsolve(a, b)) 

446 True 

447 

448 """ 

449 a = asarray(a) 

450 oldshape = a.shape 

451 prod = 1 

452 if ind > 0: 

453 invshape = oldshape[ind:] + oldshape[:ind] 

454 for k in oldshape[ind:]: 

455 prod *= k 

456 else: 

457 raise ValueError("Invalid ind argument.") 

458 a = a.reshape(prod, -1) 

459 ia = inv(a) 

460 return ia.reshape(*invshape) 

461 

462 

463# Matrix inversion 

464 

465def _unary_dispatcher(a): 

466 return (a,) 

467 

468 

469@array_function_dispatch(_unary_dispatcher) 

470def inv(a): 

471 """ 

472 Compute the (multiplicative) inverse of a matrix. 

473 

474 Given a square matrix `a`, return the matrix `ainv` satisfying 

475 ``dot(a, ainv) = dot(ainv, a) = eye(a.shape[0])``. 

476 

477 Parameters 

478 ---------- 

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

480 Matrix to be inverted. 

481 

482 Returns 

483 ------- 

484 ainv : (..., M, M) ndarray or matrix 

485 (Multiplicative) inverse of the matrix `a`. 

486 

487 Raises 

488 ------ 

489 LinAlgError 

490 If `a` is not square or inversion fails. 

491 

492 See Also 

493 -------- 

494 scipy.linalg.inv : Similar function in SciPy. 

495 

496 Notes 

497 ----- 

498 

499 .. versionadded:: 1.8.0 

500 

501 Broadcasting rules apply, see the `numpy.linalg` documentation for 

502 details. 

503 

504 Examples 

505 -------- 

506 >>> from numpy.linalg import inv 

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

508 >>> ainv = inv(a) 

509 >>> np.allclose(np.dot(a, ainv), np.eye(2)) 

510 True 

511 >>> np.allclose(np.dot(ainv, a), np.eye(2)) 

512 True 

513 

514 If a is a matrix object, then the return value is a matrix as well: 

515 

516 >>> ainv = inv(np.matrix(a)) 

517 >>> ainv 

518 matrix([[-2. , 1. ], 

519 [ 1.5, -0.5]]) 

520 

521 Inverses of several matrices can be computed at once: 

522 

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

524 >>> inv(a) 

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

526 [ 1.5 , -0.5 ]], 

527 [[-1.25, 0.75], 

528 [ 0.75, -0.25]]]) 

529 

530 """ 

531 a, wrap = _makearray(a) 

532 _assert_stacked_2d(a) 

533 _assert_stacked_square(a) 

534 t, result_t = _commonType(a) 

535 

536 signature = 'D->D' if isComplexType(t) else 'd->d' 

537 extobj = get_linalg_error_extobj(_raise_linalgerror_singular) 

538 ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj) 

539 return wrap(ainv.astype(result_t, copy=False)) 

540 

541 

542def _matrix_power_dispatcher(a, n): 

543 return (a,) 

544 

545 

546@array_function_dispatch(_matrix_power_dispatcher) 

547def matrix_power(a, n): 

548 """ 

549 Raise a square matrix to the (integer) power `n`. 

550 

551 For positive integers `n`, the power is computed by repeated matrix 

552 squarings and matrix multiplications. If ``n == 0``, the identity matrix 

553 of the same shape as M is returned. If ``n < 0``, the inverse 

554 is computed and then raised to the ``abs(n)``. 

555 

556 .. note:: Stacks of object matrices are not currently supported. 

557 

558 Parameters 

559 ---------- 

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

561 Matrix to be "powered". 

562 n : int 

563 The exponent can be any integer or long integer, positive, 

564 negative, or zero. 

565 

566 Returns 

567 ------- 

568 a**n : (..., M, M) ndarray or matrix object 

569 The return value is the same shape and type as `M`; 

570 if the exponent is positive or zero then the type of the 

571 elements is the same as those of `M`. If the exponent is 

572 negative the elements are floating-point. 

573 

574 Raises 

575 ------ 

576 LinAlgError 

577 For matrices that are not square or that (for negative powers) cannot 

578 be inverted numerically. 

579 

580 Examples 

581 -------- 

582 >>> from numpy.linalg import matrix_power 

583 >>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit 

584 >>> matrix_power(i, 3) # should = -i 

585 array([[ 0, -1], 

586 [ 1, 0]]) 

587 >>> matrix_power(i, 0) 

588 array([[1, 0], 

589 [0, 1]]) 

590 >>> matrix_power(i, -3) # should = 1/(-i) = i, but w/ f.p. elements 

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

592 [-1., 0.]]) 

593 

594 Somewhat more sophisticated example 

595 

596 >>> q = np.zeros((4, 4)) 

597 >>> q[0:2, 0:2] = -i 

598 >>> q[2:4, 2:4] = i 

599 >>> q # one of the three quaternion units not equal to 1 

600 array([[ 0., -1., 0., 0.], 

601 [ 1., 0., 0., 0.], 

602 [ 0., 0., 0., 1.], 

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

604 >>> matrix_power(q, 2) # = -np.eye(4) 

605 array([[-1., 0., 0., 0.], 

606 [ 0., -1., 0., 0.], 

607 [ 0., 0., -1., 0.], 

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

609 

610 """ 

611 a = asanyarray(a) 

612 _assert_stacked_2d(a) 

613 _assert_stacked_square(a) 

614 

615 try: 

616 n = operator.index(n) 

617 except TypeError as e: 

618 raise TypeError("exponent must be an integer") from e 

619 

620 # Fall back on dot for object arrays. Object arrays are not supported by 

621 # the current implementation of matmul using einsum 

622 if a.dtype != object: 

623 fmatmul = matmul 

624 elif a.ndim == 2: 

625 fmatmul = dot 

626 else: 

627 raise NotImplementedError( 

628 "matrix_power not supported for stacks of object arrays") 

629 

630 if n == 0: 

631 a = empty_like(a) 

632 a[...] = eye(a.shape[-2], dtype=a.dtype) 

633 return a 

634 

635 elif n < 0: 

636 a = inv(a) 

637 n = abs(n) 

638 

639 # short-cuts. 

640 if n == 1: 

641 return a 

642 

643 elif n == 2: 

644 return fmatmul(a, a) 

645 

646 elif n == 3: 

647 return fmatmul(fmatmul(a, a), a) 

648 

649 # Use binary decomposition to reduce the number of matrix multiplications. 

650 # Here, we iterate over the bits of n, from LSB to MSB, raise `a` to 

651 # increasing powers of 2, and multiply into the result as needed. 

652 z = result = None 

653 while n > 0: 

654 z = a if z is None else fmatmul(z, z) 

655 n, bit = divmod(n, 2) 

656 if bit: 

657 result = z if result is None else fmatmul(result, z) 

658 

659 return result 

660 

661 

662# Cholesky decomposition 

663 

664 

665@array_function_dispatch(_unary_dispatcher) 

666def cholesky(a): 

667 """ 

668 Cholesky decomposition. 

669 

670 Return the Cholesky decomposition, `L * L.H`, of the square matrix `a`, 

671 where `L` is lower-triangular and .H is the conjugate transpose operator 

672 (which is the ordinary transpose if `a` is real-valued). `a` must be 

673 Hermitian (symmetric if real-valued) and positive-definite. No 

674 checking is performed to verify whether `a` is Hermitian or not. 

675 In addition, only the lower-triangular and diagonal elements of `a` 

676 are used. Only `L` is actually returned. 

677 

678 Parameters 

679 ---------- 

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

681 Hermitian (symmetric if all elements are real), positive-definite 

682 input matrix. 

683 

684 Returns 

685 ------- 

686 L : (..., M, M) array_like 

687 Lower-triangular Cholesky factor of `a`. Returns a matrix object if 

688 `a` is a matrix object. 

689 

690 Raises 

691 ------ 

692 LinAlgError 

693 If the decomposition fails, for example, if `a` is not 

694 positive-definite. 

695 

696 See Also 

697 -------- 

698 scipy.linalg.cholesky : Similar function in SciPy. 

699 scipy.linalg.cholesky_banded : Cholesky decompose a banded Hermitian 

700 positive-definite matrix. 

701 scipy.linalg.cho_factor : Cholesky decomposition of a matrix, to use in 

702 `scipy.linalg.cho_solve`. 

703 

704 Notes 

705 ----- 

706 

707 .. versionadded:: 1.8.0 

708 

709 Broadcasting rules apply, see the `numpy.linalg` documentation for 

710 details. 

711 

712 The Cholesky decomposition is often used as a fast way of solving 

713 

714 .. math:: A \\mathbf{x} = \\mathbf{b} 

715 

716 (when `A` is both Hermitian/symmetric and positive-definite). 

717 

718 First, we solve for :math:`\\mathbf{y}` in 

719 

720 .. math:: L \\mathbf{y} = \\mathbf{b}, 

721 

722 and then for :math:`\\mathbf{x}` in 

723 

724 .. math:: L.H \\mathbf{x} = \\mathbf{y}. 

725 

726 Examples 

727 -------- 

728 >>> A = np.array([[1,-2j],[2j,5]]) 

729 >>> A 

730 array([[ 1.+0.j, -0.-2.j], 

731 [ 0.+2.j, 5.+0.j]]) 

732 >>> L = np.linalg.cholesky(A) 

733 >>> L 

734 array([[1.+0.j, 0.+0.j], 

735 [0.+2.j, 1.+0.j]]) 

736 >>> np.dot(L, L.T.conj()) # verify that L * L.H = A 

737 array([[1.+0.j, 0.-2.j], 

738 [0.+2.j, 5.+0.j]]) 

739 >>> A = [[1,-2j],[2j,5]] # what happens if A is only array_like? 

740 >>> np.linalg.cholesky(A) # an ndarray object is returned 

741 array([[1.+0.j, 0.+0.j], 

742 [0.+2.j, 1.+0.j]]) 

743 >>> # But a matrix object is returned if A is a matrix object 

744 >>> np.linalg.cholesky(np.matrix(A)) 

745 matrix([[ 1.+0.j, 0.+0.j], 

746 [ 0.+2.j, 1.+0.j]]) 

747 

748 """ 

749 extobj = get_linalg_error_extobj(_raise_linalgerror_nonposdef) 

750 gufunc = _umath_linalg.cholesky_lo 

751 a, wrap = _makearray(a) 

752 _assert_stacked_2d(a) 

753 _assert_stacked_square(a) 

754 t, result_t = _commonType(a) 

755 signature = 'D->D' if isComplexType(t) else 'd->d' 

756 r = gufunc(a, signature=signature, extobj=extobj) 

757 return wrap(r.astype(result_t, copy=False)) 

758 

759 

760# QR decomposition 

761 

762def _qr_dispatcher(a, mode=None): 

763 return (a,) 

764 

765 

766@array_function_dispatch(_qr_dispatcher) 

767def qr(a, mode='reduced'): 

768 """ 

769 Compute the qr factorization of a matrix. 

770 

771 Factor the matrix `a` as *qr*, where `q` is orthonormal and `r` is 

772 upper-triangular. 

773 

774 Parameters 

775 ---------- 

776 a : array_like, shape (..., M, N) 

777 An array-like object with the dimensionality of at least 2. 

778 mode : {'reduced', 'complete', 'r', 'raw'}, optional 

779 If K = min(M, N), then 

780 

781 * 'reduced' : returns q, r with dimensions 

782 (..., M, K), (..., K, N) (default) 

783 * 'complete' : returns q, r with dimensions (..., M, M), (..., M, N) 

784 * 'r' : returns r only with dimensions (..., K, N) 

785 * 'raw' : returns h, tau with dimensions (..., N, M), (..., K,) 

786 

787 The options 'reduced', 'complete, and 'raw' are new in numpy 1.8, 

788 see the notes for more information. The default is 'reduced', and to 

789 maintain backward compatibility with earlier versions of numpy both 

790 it and the old default 'full' can be omitted. Note that array h 

791 returned in 'raw' mode is transposed for calling Fortran. The 

792 'economic' mode is deprecated. The modes 'full' and 'economic' may 

793 be passed using only the first letter for backwards compatibility, 

794 but all others must be spelled out. See the Notes for more 

795 explanation. 

796 

797 

798 Returns 

799 ------- 

800 q : ndarray of float or complex, optional 

801 A matrix with orthonormal columns. When mode = 'complete' the 

802 result is an orthogonal/unitary matrix depending on whether or not 

803 a is real/complex. The determinant may be either +/- 1 in that 

804 case. In case the number of dimensions in the input array is 

805 greater than 2 then a stack of the matrices with above properties 

806 is returned. 

807 r : ndarray of float or complex, optional 

808 The upper-triangular matrix or a stack of upper-triangular 

809 matrices if the number of dimensions in the input array is greater 

810 than 2. 

811 (h, tau) : ndarrays of np.double or np.cdouble, optional 

812 The array h contains the Householder reflectors that generate q 

813 along with r. The tau array contains scaling factors for the 

814 reflectors. In the deprecated 'economic' mode only h is returned. 

815 

816 Raises 

817 ------ 

818 LinAlgError 

819 If factoring fails. 

820 

821 See Also 

822 -------- 

823 scipy.linalg.qr : Similar function in SciPy. 

824 scipy.linalg.rq : Compute RQ decomposition of a matrix. 

825 

826 Notes 

827 ----- 

828 This is an interface to the LAPACK routines ``dgeqrf``, ``zgeqrf``, 

829 ``dorgqr``, and ``zungqr``. 

830 

831 For more information on the qr factorization, see for example: 

832 https://en.wikipedia.org/wiki/QR_factorization 

833 

834 Subclasses of `ndarray` are preserved except for the 'raw' mode. So if 

835 `a` is of type `matrix`, all the return values will be matrices too. 

836 

837 New 'reduced', 'complete', and 'raw' options for mode were added in 

838 NumPy 1.8.0 and the old option 'full' was made an alias of 'reduced'. In 

839 addition the options 'full' and 'economic' were deprecated. Because 

840 'full' was the previous default and 'reduced' is the new default, 

841 backward compatibility can be maintained by letting `mode` default. 

842 The 'raw' option was added so that LAPACK routines that can multiply 

843 arrays by q using the Householder reflectors can be used. Note that in 

844 this case the returned arrays are of type np.double or np.cdouble and 

845 the h array is transposed to be FORTRAN compatible. No routines using 

846 the 'raw' return are currently exposed by numpy, but some are available 

847 in lapack_lite and just await the necessary work. 

848 

849 Examples 

850 -------- 

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

852 >>> q, r = np.linalg.qr(a) 

853 >>> np.allclose(a, np.dot(q, r)) # a does equal qr 

854 True 

855 >>> r2 = np.linalg.qr(a, mode='r') 

856 >>> np.allclose(r, r2) # mode='r' returns the same r as mode='full' 

857 True 

858 >>> a = np.random.normal(size=(3, 2, 2)) # Stack of 2 x 2 matrices as input 

859 >>> q, r = np.linalg.qr(a) 

860 >>> q.shape 

861 (3, 2, 2) 

862 >>> r.shape 

863 (3, 2, 2) 

864 >>> np.allclose(a, np.matmul(q, r)) 

865 True 

866 

867 Example illustrating a common use of `qr`: solving of least squares 

868 problems 

869 

870 What are the least-squares-best `m` and `y0` in ``y = y0 + mx`` for 

871 the following data: {(0,1), (1,0), (1,2), (2,1)}. (Graph the points 

872 and you'll see that it should be y0 = 0, m = 1.) The answer is provided 

873 by solving the over-determined matrix equation ``Ax = b``, where:: 

874 

875 A = array([[0, 1], [1, 1], [1, 1], [2, 1]]) 

876 x = array([[y0], [m]]) 

877 b = array([[1], [0], [2], [1]]) 

878 

879 If A = qr such that q is orthonormal (which is always possible via 

880 Gram-Schmidt), then ``x = inv(r) * (q.T) * b``. (In numpy practice, 

881 however, we simply use `lstsq`.) 

882 

883 >>> A = np.array([[0, 1], [1, 1], [1, 1], [2, 1]]) 

884 >>> A 

885 array([[0, 1], 

886 [1, 1], 

887 [1, 1], 

888 [2, 1]]) 

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

890 >>> q, r = np.linalg.qr(A) 

891 >>> p = np.dot(q.T, b) 

892 >>> np.dot(np.linalg.inv(r), p) 

893 array([ 1., 1.]) 

894 

895 """ 

896 if mode not in ('reduced', 'complete', 'r', 'raw'): 

897 if mode in ('f', 'full'): 

898 # 2013-04-01, 1.8 

899 msg = "".join(( 

900 "The 'full' option is deprecated in favor of 'reduced'.\n", 

901 "For backward compatibility let mode default.")) 

902 warnings.warn(msg, DeprecationWarning, stacklevel=3) 

903 mode = 'reduced' 

904 elif mode in ('e', 'economic'): 

905 # 2013-04-01, 1.8 

906 msg = "The 'economic' option is deprecated." 

907 warnings.warn(msg, DeprecationWarning, stacklevel=3) 

908 mode = 'economic' 

909 else: 

910 raise ValueError(f"Unrecognized mode '{mode}'") 

911 

912 a, wrap = _makearray(a) 

913 _assert_stacked_2d(a) 

914 m, n = a.shape[-2:] 

915 t, result_t = _commonType(a) 

916 a = a.astype(t, copy=True) 

917 a = _to_native_byte_order(a) 

918 mn = min(m, n) 

919 

920 if m <= n: 

921 gufunc = _umath_linalg.qr_r_raw_m 

922 else: 

923 gufunc = _umath_linalg.qr_r_raw_n 

924 

925 signature = 'D->D' if isComplexType(t) else 'd->d' 

926 extobj = get_linalg_error_extobj(_raise_linalgerror_qr) 

927 tau = gufunc(a, signature=signature, extobj=extobj) 

928 

929 # handle modes that don't return q 

930 if mode == 'r': 

931 r = triu(a[..., :mn, :]) 

932 r = r.astype(result_t, copy=False) 

933 return wrap(r) 

934 

935 if mode == 'raw': 

936 q = transpose(a) 

937 q = q.astype(result_t, copy=False) 

938 tau = tau.astype(result_t, copy=False) 

939 return wrap(q), tau 

940 

941 if mode == 'economic': 

942 a = a.astype(result_t, copy=False) 

943 return wrap(a) 

944 

945 # mc is the number of columns in the resulting q 

946 # matrix. If the mode is complete then it is  

947 # same as number of rows, and if the mode is reduced, 

948 # then it is the minimum of number of rows and columns. 

949 if mode == 'complete' and m > n: 

950 mc = m 

951 gufunc = _umath_linalg.qr_complete 

952 else: 

953 mc = mn 

954 gufunc = _umath_linalg.qr_reduced 

955 

956 signature = 'DD->D' if isComplexType(t) else 'dd->d' 

957 extobj = get_linalg_error_extobj(_raise_linalgerror_qr) 

958 q = gufunc(a, tau, signature=signature, extobj=extobj) 

959 r = triu(a[..., :mc, :]) 

960 

961 q = q.astype(result_t, copy=False) 

962 r = r.astype(result_t, copy=False) 

963 

964 return wrap(q), wrap(r) 

965 

966# Eigenvalues 

967 

968 

969@array_function_dispatch(_unary_dispatcher) 

970def eigvals(a): 

971 """ 

972 Compute the eigenvalues of a general matrix. 

973 

974 Main difference between `eigvals` and `eig`: the eigenvectors aren't 

975 returned. 

976 

977 Parameters 

978 ---------- 

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

980 A complex- or real-valued matrix whose eigenvalues will be computed. 

981 

982 Returns 

983 ------- 

984 w : (..., M,) ndarray 

985 The eigenvalues, each repeated according to its multiplicity. 

986 They are not necessarily ordered, nor are they necessarily 

987 real for real matrices. 

988 

989 Raises 

990 ------ 

991 LinAlgError 

992 If the eigenvalue computation does not converge. 

993 

994 See Also 

995 -------- 

996 eig : eigenvalues and right eigenvectors of general arrays 

997 eigvalsh : eigenvalues of real symmetric or complex Hermitian 

998 (conjugate symmetric) arrays. 

999 eigh : eigenvalues and eigenvectors of real symmetric or complex 

1000 Hermitian (conjugate symmetric) arrays. 

1001 scipy.linalg.eigvals : Similar function in SciPy. 

1002 

1003 Notes 

1004 ----- 

1005 

1006 .. versionadded:: 1.8.0 

1007 

1008 Broadcasting rules apply, see the `numpy.linalg` documentation for 

1009 details. 

1010 

1011 This is implemented using the ``_geev`` LAPACK routines which compute 

1012 the eigenvalues and eigenvectors of general square arrays. 

1013 

1014 Examples 

1015 -------- 

1016 Illustration, using the fact that the eigenvalues of a diagonal matrix 

1017 are its diagonal elements, that multiplying a matrix on the left 

1018 by an orthogonal matrix, `Q`, and on the right by `Q.T` (the transpose 

1019 of `Q`), preserves the eigenvalues of the "middle" matrix. In other words, 

1020 if `Q` is orthogonal, then ``Q * A * Q.T`` has the same eigenvalues as 

1021 ``A``: 

1022 

1023 >>> from numpy import linalg as LA 

1024 >>> x = np.random.random() 

1025 >>> Q = np.array([[np.cos(x), -np.sin(x)], [np.sin(x), np.cos(x)]]) 

1026 >>> LA.norm(Q[0, :]), LA.norm(Q[1, :]), np.dot(Q[0, :],Q[1, :]) 

1027 (1.0, 1.0, 0.0) 

1028 

1029 Now multiply a diagonal matrix by ``Q`` on one side and by ``Q.T`` on the other: 

1030 

1031 >>> D = np.diag((-1,1)) 

1032 >>> LA.eigvals(D) 

1033 array([-1., 1.]) 

1034 >>> A = np.dot(Q, D) 

1035 >>> A = np.dot(A, Q.T) 

1036 >>> LA.eigvals(A) 

1037 array([ 1., -1.]) # random 

1038 

1039 """ 

1040 a, wrap = _makearray(a) 

1041 _assert_stacked_2d(a) 

1042 _assert_stacked_square(a) 

1043 _assert_finite(a) 

1044 t, result_t = _commonType(a) 

1045 

1046 extobj = get_linalg_error_extobj( 

1047 _raise_linalgerror_eigenvalues_nonconvergence) 

1048 signature = 'D->D' if isComplexType(t) else 'd->D' 

1049 w = _umath_linalg.eigvals(a, signature=signature, extobj=extobj) 

1050 

1051 if not isComplexType(t): 

1052 if all(w.imag == 0): 

1053 w = w.real 

1054 result_t = _realType(result_t) 

1055 else: 

1056 result_t = _complexType(result_t) 

1057 

1058 return w.astype(result_t, copy=False) 

1059 

1060 

1061def _eigvalsh_dispatcher(a, UPLO=None): 

1062 return (a,) 

1063 

1064 

1065@array_function_dispatch(_eigvalsh_dispatcher) 

1066def eigvalsh(a, UPLO='L'): 

1067 """ 

1068 Compute the eigenvalues of a complex Hermitian or real symmetric matrix. 

1069 

1070 Main difference from eigh: the eigenvectors are not computed. 

1071 

1072 Parameters 

1073 ---------- 

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

1075 A complex- or real-valued matrix whose eigenvalues are to be 

1076 computed. 

1077 UPLO : {'L', 'U'}, optional 

1078 Specifies whether the calculation is done with the lower triangular 

1079 part of `a` ('L', default) or the upper triangular part ('U'). 

1080 Irrespective of this value only the real parts of the diagonal will 

1081 be considered in the computation to preserve the notion of a Hermitian 

1082 matrix. It therefore follows that the imaginary part of the diagonal 

1083 will always be treated as zero. 

1084 

1085 Returns 

1086 ------- 

1087 w : (..., M,) ndarray 

1088 The eigenvalues in ascending order, each repeated according to 

1089 its multiplicity. 

1090 

1091 Raises 

1092 ------ 

1093 LinAlgError 

1094 If the eigenvalue computation does not converge. 

1095 

1096 See Also 

1097 -------- 

1098 eigh : eigenvalues and eigenvectors of real symmetric or complex Hermitian 

1099 (conjugate symmetric) arrays. 

1100 eigvals : eigenvalues of general real or complex arrays. 

1101 eig : eigenvalues and right eigenvectors of general real or complex 

1102 arrays. 

1103 scipy.linalg.eigvalsh : Similar function in SciPy. 

1104 

1105 Notes 

1106 ----- 

1107 

1108 .. versionadded:: 1.8.0 

1109 

1110 Broadcasting rules apply, see the `numpy.linalg` documentation for 

1111 details. 

1112 

1113 The eigenvalues are computed using LAPACK routines ``_syevd``, ``_heevd``. 

1114 

1115 Examples 

1116 -------- 

1117 >>> from numpy import linalg as LA 

1118 >>> a = np.array([[1, -2j], [2j, 5]]) 

1119 >>> LA.eigvalsh(a) 

1120 array([ 0.17157288, 5.82842712]) # may vary 

1121 

1122 >>> # demonstrate the treatment of the imaginary part of the diagonal 

1123 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]]) 

1124 >>> a 

1125 array([[5.+2.j, 9.-2.j], 

1126 [0.+2.j, 2.-1.j]]) 

1127 >>> # with UPLO='L' this is numerically equivalent to using LA.eigvals() 

1128 >>> # with: 

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

1130 >>> b 

1131 array([[5.+0.j, 0.-2.j], 

1132 [0.+2.j, 2.+0.j]]) 

1133 >>> wa = LA.eigvalsh(a) 

1134 >>> wb = LA.eigvals(b) 

1135 >>> wa; wb 

1136 array([1., 6.]) 

1137 array([6.+0.j, 1.+0.j]) 

1138 

1139 """ 

1140 UPLO = UPLO.upper() 

1141 if UPLO not in ('L', 'U'): 

1142 raise ValueError("UPLO argument must be 'L' or 'U'") 

1143 

1144 extobj = get_linalg_error_extobj( 

1145 _raise_linalgerror_eigenvalues_nonconvergence) 

1146 if UPLO == 'L': 

1147 gufunc = _umath_linalg.eigvalsh_lo 

1148 else: 

1149 gufunc = _umath_linalg.eigvalsh_up 

1150 

1151 a, wrap = _makearray(a) 

1152 _assert_stacked_2d(a) 

1153 _assert_stacked_square(a) 

1154 t, result_t = _commonType(a) 

1155 signature = 'D->d' if isComplexType(t) else 'd->d' 

1156 w = gufunc(a, signature=signature, extobj=extobj) 

1157 return w.astype(_realType(result_t), copy=False) 

1158 

1159def _convertarray(a): 

1160 t, result_t = _commonType(a) 

1161 a = a.astype(t).T.copy() 

1162 return a, t, result_t 

1163 

1164 

1165# Eigenvectors 

1166 

1167 

1168@array_function_dispatch(_unary_dispatcher) 

1169def eig(a): 

1170 """ 

1171 Compute the eigenvalues and right eigenvectors of a square array. 

1172 

1173 Parameters 

1174 ---------- 

1175 a : (..., M, M) array 

1176 Matrices for which the eigenvalues and right eigenvectors will 

1177 be computed 

1178 

1179 Returns 

1180 ------- 

1181 w : (..., M) array 

1182 The eigenvalues, each repeated according to its multiplicity. 

1183 The eigenvalues are not necessarily ordered. The resulting 

1184 array will be of complex type, unless the imaginary part is 

1185 zero in which case it will be cast to a real type. When `a` 

1186 is real the resulting eigenvalues will be real (0 imaginary 

1187 part) or occur in conjugate pairs 

1188 

1189 v : (..., M, M) array 

1190 The normalized (unit "length") eigenvectors, such that the 

1191 column ``v[:,i]`` is the eigenvector corresponding to the 

1192 eigenvalue ``w[i]``. 

1193 

1194 Raises 

1195 ------ 

1196 LinAlgError 

1197 If the eigenvalue computation does not converge. 

1198 

1199 See Also 

1200 -------- 

1201 eigvals : eigenvalues of a non-symmetric array. 

1202 eigh : eigenvalues and eigenvectors of a real symmetric or complex 

1203 Hermitian (conjugate symmetric) array. 

1204 eigvalsh : eigenvalues of a real symmetric or complex Hermitian 

1205 (conjugate symmetric) array. 

1206 scipy.linalg.eig : Similar function in SciPy that also solves the 

1207 generalized eigenvalue problem. 

1208 scipy.linalg.schur : Best choice for unitary and other non-Hermitian 

1209 normal matrices. 

1210 

1211 Notes 

1212 ----- 

1213 

1214 .. versionadded:: 1.8.0 

1215 

1216 Broadcasting rules apply, see the `numpy.linalg` documentation for 

1217 details. 

1218 

1219 This is implemented using the ``_geev`` LAPACK routines which compute 

1220 the eigenvalues and eigenvectors of general square arrays. 

1221 

1222 The number `w` is an eigenvalue of `a` if there exists a vector 

1223 `v` such that ``a @ v = w * v``. Thus, the arrays `a`, `w`, and 

1224 `v` satisfy the equations ``a @ v[:,i] = w[i] * v[:,i]`` 

1225 for :math:`i \\in \\{0,...,M-1\\}`. 

1226 

1227 The array `v` of eigenvectors may not be of maximum rank, that is, some 

1228 of the columns may be linearly dependent, although round-off error may 

1229 obscure that fact. If the eigenvalues are all different, then theoretically 

1230 the eigenvectors are linearly independent and `a` can be diagonalized by 

1231 a similarity transformation using `v`, i.e, ``inv(v) @ a @ v`` is diagonal. 

1232 

1233 For non-Hermitian normal matrices the SciPy function `scipy.linalg.schur` 

1234 is preferred because the matrix `v` is guaranteed to be unitary, which is 

1235 not the case when using `eig`. The Schur factorization produces an 

1236 upper triangular matrix rather than a diagonal matrix, but for normal 

1237 matrices only the diagonal of the upper triangular matrix is needed, the 

1238 rest is roundoff error. 

1239 

1240 Finally, it is emphasized that `v` consists of the *right* (as in 

1241 right-hand side) eigenvectors of `a`. A vector `y` satisfying 

1242 ``y.T @ a = z * y.T`` for some number `z` is called a *left* 

1243 eigenvector of `a`, and, in general, the left and right eigenvectors 

1244 of a matrix are not necessarily the (perhaps conjugate) transposes 

1245 of each other. 

1246 

1247 References 

1248 ---------- 

1249 G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, FL, 

1250 Academic Press, Inc., 1980, Various pp. 

1251 

1252 Examples 

1253 -------- 

1254 >>> from numpy import linalg as LA 

1255 

1256 (Almost) trivial example with real e-values and e-vectors. 

1257 

1258 >>> w, v = LA.eig(np.diag((1, 2, 3))) 

1259 >>> w; v 

1260 array([1., 2., 3.]) 

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

1262 [0., 1., 0.], 

1263 [0., 0., 1.]]) 

1264 

1265 Real matrix possessing complex e-values and e-vectors; note that the 

1266 e-values are complex conjugates of each other. 

1267 

1268 >>> w, v = LA.eig(np.array([[1, -1], [1, 1]])) 

1269 >>> w; v 

1270 array([1.+1.j, 1.-1.j]) 

1271 array([[0.70710678+0.j , 0.70710678-0.j ], 

1272 [0. -0.70710678j, 0. +0.70710678j]]) 

1273 

1274 Complex-valued matrix with real e-values (but complex-valued e-vectors); 

1275 note that ``a.conj().T == a``, i.e., `a` is Hermitian. 

1276 

1277 >>> a = np.array([[1, 1j], [-1j, 1]]) 

1278 >>> w, v = LA.eig(a) 

1279 >>> w; v 

1280 array([2.+0.j, 0.+0.j]) 

1281 array([[ 0. +0.70710678j, 0.70710678+0.j ], # may vary 

1282 [ 0.70710678+0.j , -0. +0.70710678j]]) 

1283 

1284 Be careful about round-off error! 

1285 

1286 >>> a = np.array([[1 + 1e-9, 0], [0, 1 - 1e-9]]) 

1287 >>> # Theor. e-values are 1 +/- 1e-9 

1288 >>> w, v = LA.eig(a) 

1289 >>> w; v 

1290 array([1., 1.]) 

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

1292 [0., 1.]]) 

1293 

1294 """ 

1295 a, wrap = _makearray(a) 

1296 _assert_stacked_2d(a) 

1297 _assert_stacked_square(a) 

1298 _assert_finite(a) 

1299 t, result_t = _commonType(a) 

1300 

1301 extobj = get_linalg_error_extobj( 

1302 _raise_linalgerror_eigenvalues_nonconvergence) 

1303 signature = 'D->DD' if isComplexType(t) else 'd->DD' 

1304 w, vt = _umath_linalg.eig(a, signature=signature, extobj=extobj) 

1305 

1306 if not isComplexType(t) and all(w.imag == 0.0): 

1307 w = w.real 

1308 vt = vt.real 

1309 result_t = _realType(result_t) 

1310 else: 

1311 result_t = _complexType(result_t) 

1312 

1313 vt = vt.astype(result_t, copy=False) 

1314 return w.astype(result_t, copy=False), wrap(vt) 

1315 

1316 

1317@array_function_dispatch(_eigvalsh_dispatcher) 

1318def eigh(a, UPLO='L'): 

1319 """ 

1320 Return the eigenvalues and eigenvectors of a complex Hermitian 

1321 (conjugate symmetric) or a real symmetric matrix. 

1322 

1323 Returns two objects, a 1-D array containing the eigenvalues of `a`, and 

1324 a 2-D square array or matrix (depending on the input type) of the 

1325 corresponding eigenvectors (in columns). 

1326 

1327 Parameters 

1328 ---------- 

1329 a : (..., M, M) array 

1330 Hermitian or real symmetric matrices whose eigenvalues and 

1331 eigenvectors are to be computed. 

1332 UPLO : {'L', 'U'}, optional 

1333 Specifies whether the calculation is done with the lower triangular 

1334 part of `a` ('L', default) or the upper triangular part ('U'). 

1335 Irrespective of this value only the real parts of the diagonal will 

1336 be considered in the computation to preserve the notion of a Hermitian 

1337 matrix. It therefore follows that the imaginary part of the diagonal 

1338 will always be treated as zero. 

1339 

1340 Returns 

1341 ------- 

1342 w : (..., M) ndarray 

1343 The eigenvalues in ascending order, each repeated according to 

1344 its multiplicity. 

1345 v : {(..., M, M) ndarray, (..., M, M) matrix} 

1346 The column ``v[:, i]`` is the normalized eigenvector corresponding 

1347 to the eigenvalue ``w[i]``. Will return a matrix object if `a` is 

1348 a matrix object. 

1349 

1350 Raises 

1351 ------ 

1352 LinAlgError 

1353 If the eigenvalue computation does not converge. 

1354 

1355 See Also 

1356 -------- 

1357 eigvalsh : eigenvalues of real symmetric or complex Hermitian 

1358 (conjugate symmetric) arrays. 

1359 eig : eigenvalues and right eigenvectors for non-symmetric arrays. 

1360 eigvals : eigenvalues of non-symmetric arrays. 

1361 scipy.linalg.eigh : Similar function in SciPy (but also solves the 

1362 generalized eigenvalue problem). 

1363 

1364 Notes 

1365 ----- 

1366 

1367 .. versionadded:: 1.8.0 

1368 

1369 Broadcasting rules apply, see the `numpy.linalg` documentation for 

1370 details. 

1371 

1372 The eigenvalues/eigenvectors are computed using LAPACK routines ``_syevd``, 

1373 ``_heevd``. 

1374 

1375 The eigenvalues of real symmetric or complex Hermitian matrices are 

1376 always real. [1]_ The array `v` of (column) eigenvectors is unitary 

1377 and `a`, `w`, and `v` satisfy the equations 

1378 ``dot(a, v[:, i]) = w[i] * v[:, i]``. 

1379 

1380 References 

1381 ---------- 

1382 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, 

1383 FL, Academic Press, Inc., 1980, pg. 222. 

1384 

1385 Examples 

1386 -------- 

1387 >>> from numpy import linalg as LA 

1388 >>> a = np.array([[1, -2j], [2j, 5]]) 

1389 >>> a 

1390 array([[ 1.+0.j, -0.-2.j], 

1391 [ 0.+2.j, 5.+0.j]]) 

1392 >>> w, v = LA.eigh(a) 

1393 >>> w; v 

1394 array([0.17157288, 5.82842712]) 

1395 array([[-0.92387953+0.j , -0.38268343+0.j ], # may vary 

1396 [ 0. +0.38268343j, 0. -0.92387953j]]) 

1397 

1398 >>> np.dot(a, v[:, 0]) - w[0] * v[:, 0] # verify 1st e-val/vec pair 

1399 array([5.55111512e-17+0.0000000e+00j, 0.00000000e+00+1.2490009e-16j]) 

1400 >>> np.dot(a, v[:, 1]) - w[1] * v[:, 1] # verify 2nd e-val/vec pair 

1401 array([0.+0.j, 0.+0.j]) 

1402 

1403 >>> A = np.matrix(a) # what happens if input is a matrix object 

1404 >>> A 

1405 matrix([[ 1.+0.j, -0.-2.j], 

1406 [ 0.+2.j, 5.+0.j]]) 

1407 >>> w, v = LA.eigh(A) 

1408 >>> w; v 

1409 array([0.17157288, 5.82842712]) 

1410 matrix([[-0.92387953+0.j , -0.38268343+0.j ], # may vary 

1411 [ 0. +0.38268343j, 0. -0.92387953j]]) 

1412 

1413 >>> # demonstrate the treatment of the imaginary part of the diagonal 

1414 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]]) 

1415 >>> a 

1416 array([[5.+2.j, 9.-2.j], 

1417 [0.+2.j, 2.-1.j]]) 

1418 >>> # with UPLO='L' this is numerically equivalent to using LA.eig() with: 

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

1420 >>> b 

1421 array([[5.+0.j, 0.-2.j], 

1422 [0.+2.j, 2.+0.j]]) 

1423 >>> wa, va = LA.eigh(a) 

1424 >>> wb, vb = LA.eig(b) 

1425 >>> wa; wb 

1426 array([1., 6.]) 

1427 array([6.+0.j, 1.+0.j]) 

1428 >>> va; vb 

1429 array([[-0.4472136 +0.j , -0.89442719+0.j ], # may vary 

1430 [ 0. +0.89442719j, 0. -0.4472136j ]]) 

1431 array([[ 0.89442719+0.j , -0. +0.4472136j], 

1432 [-0. +0.4472136j, 0.89442719+0.j ]]) 

1433 """ 

1434 UPLO = UPLO.upper() 

1435 if UPLO not in ('L', 'U'): 

1436 raise ValueError("UPLO argument must be 'L' or 'U'") 

1437 

1438 a, wrap = _makearray(a) 

1439 _assert_stacked_2d(a) 

1440 _assert_stacked_square(a) 

1441 t, result_t = _commonType(a) 

1442 

1443 extobj = get_linalg_error_extobj( 

1444 _raise_linalgerror_eigenvalues_nonconvergence) 

1445 if UPLO == 'L': 

1446 gufunc = _umath_linalg.eigh_lo 

1447 else: 

1448 gufunc = _umath_linalg.eigh_up 

1449 

1450 signature = 'D->dD' if isComplexType(t) else 'd->dd' 

1451 w, vt = gufunc(a, signature=signature, extobj=extobj) 

1452 w = w.astype(_realType(result_t), copy=False) 

1453 vt = vt.astype(result_t, copy=False) 

1454 return w, wrap(vt) 

1455 

1456 

1457# Singular value decomposition 

1458 

1459def _svd_dispatcher(a, full_matrices=None, compute_uv=None, hermitian=None): 

1460 return (a,) 

1461 

1462 

1463@array_function_dispatch(_svd_dispatcher) 

1464def svd(a, full_matrices=True, compute_uv=True, hermitian=False): 

1465 """ 

1466 Singular Value Decomposition. 

1467 

1468 When `a` is a 2D array, and ``full_matrices=False``, then it is 

1469 factorized as ``u @ np.diag(s) @ vh = (u * s) @ vh``, where 

1470 `u` and the Hermitian transpose of `vh` are 2D arrays with 

1471 orthonormal columns and `s` is a 1D array of `a`'s singular 

1472 values. When `a` is higher-dimensional, SVD is applied in 

1473 stacked mode as explained below. 

1474 

1475 Parameters 

1476 ---------- 

1477 a : (..., M, N) array_like 

1478 A real or complex array with ``a.ndim >= 2``. 

1479 full_matrices : bool, optional 

1480 If True (default), `u` and `vh` have the shapes ``(..., M, M)`` and 

1481 ``(..., N, N)``, respectively. Otherwise, the shapes are 

1482 ``(..., M, K)`` and ``(..., K, N)``, respectively, where 

1483 ``K = min(M, N)``. 

1484 compute_uv : bool, optional 

1485 Whether or not to compute `u` and `vh` in addition to `s`. True 

1486 by default. 

1487 hermitian : bool, optional 

1488 If True, `a` is assumed to be Hermitian (symmetric if real-valued), 

1489 enabling a more efficient method for finding singular values. 

1490 Defaults to False. 

1491 

1492 .. versionadded:: 1.17.0 

1493 

1494 Returns 

1495 ------- 

1496 u : { (..., M, M), (..., M, K) } array 

1497 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same 

1498 size as those of the input `a`. The size of the last two dimensions 

1499 depends on the value of `full_matrices`. Only returned when 

1500 `compute_uv` is True. 

1501 s : (..., K) array 

1502 Vector(s) with the singular values, within each vector sorted in 

1503 descending order. The first ``a.ndim - 2`` dimensions have the same 

1504 size as those of the input `a`. 

1505 vh : { (..., N, N), (..., K, N) } array 

1506 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same 

1507 size as those of the input `a`. The size of the last two dimensions 

1508 depends on the value of `full_matrices`. Only returned when 

1509 `compute_uv` is True. 

1510 

1511 Raises 

1512 ------ 

1513 LinAlgError 

1514 If SVD computation does not converge. 

1515 

1516 See Also 

1517 -------- 

1518 scipy.linalg.svd : Similar function in SciPy. 

1519 scipy.linalg.svdvals : Compute singular values of a matrix. 

1520 

1521 Notes 

1522 ----- 

1523 

1524 .. versionchanged:: 1.8.0 

1525 Broadcasting rules apply, see the `numpy.linalg` documentation for 

1526 details. 

1527 

1528 The decomposition is performed using LAPACK routine ``_gesdd``. 

1529 

1530 SVD is usually described for the factorization of a 2D matrix :math:`A`. 

1531 The higher-dimensional case will be discussed below. In the 2D case, SVD is 

1532 written as :math:`A = U S V^H`, where :math:`A = a`, :math:`U= u`, 

1533 :math:`S= \\mathtt{np.diag}(s)` and :math:`V^H = vh`. The 1D array `s` 

1534 contains the singular values of `a` and `u` and `vh` are unitary. The rows 

1535 of `vh` are the eigenvectors of :math:`A^H A` and the columns of `u` are 

1536 the eigenvectors of :math:`A A^H`. In both cases the corresponding 

1537 (possibly non-zero) eigenvalues are given by ``s**2``. 

1538 

1539 If `a` has more than two dimensions, then broadcasting rules apply, as 

1540 explained in :ref:`routines.linalg-broadcasting`. This means that SVD is 

1541 working in "stacked" mode: it iterates over all indices of the first 

1542 ``a.ndim - 2`` dimensions and for each combination SVD is applied to the 

1543 last two indices. The matrix `a` can be reconstructed from the 

1544 decomposition with either ``(u * s[..., None, :]) @ vh`` or 

1545 ``u @ (s[..., None] * vh)``. (The ``@`` operator can be replaced by the 

1546 function ``np.matmul`` for python versions below 3.5.) 

1547 

1548 If `a` is a ``matrix`` object (as opposed to an ``ndarray``), then so are 

1549 all the return values. 

1550 

1551 Examples 

1552 -------- 

1553 >>> a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6) 

1554 >>> b = np.random.randn(2, 7, 8, 3) + 1j*np.random.randn(2, 7, 8, 3) 

1555 

1556 Reconstruction based on full SVD, 2D case: 

1557 

1558 >>> u, s, vh = np.linalg.svd(a, full_matrices=True) 

1559 >>> u.shape, s.shape, vh.shape 

1560 ((9, 9), (6,), (6, 6)) 

1561 >>> np.allclose(a, np.dot(u[:, :6] * s, vh)) 

1562 True 

1563 >>> smat = np.zeros((9, 6), dtype=complex) 

1564 >>> smat[:6, :6] = np.diag(s) 

1565 >>> np.allclose(a, np.dot(u, np.dot(smat, vh))) 

1566 True 

1567 

1568 Reconstruction based on reduced SVD, 2D case: 

1569 

1570 >>> u, s, vh = np.linalg.svd(a, full_matrices=False) 

1571 >>> u.shape, s.shape, vh.shape 

1572 ((9, 6), (6,), (6, 6)) 

1573 >>> np.allclose(a, np.dot(u * s, vh)) 

1574 True 

1575 >>> smat = np.diag(s) 

1576 >>> np.allclose(a, np.dot(u, np.dot(smat, vh))) 

1577 True 

1578 

1579 Reconstruction based on full SVD, 4D case: 

1580 

1581 >>> u, s, vh = np.linalg.svd(b, full_matrices=True) 

1582 >>> u.shape, s.shape, vh.shape 

1583 ((2, 7, 8, 8), (2, 7, 3), (2, 7, 3, 3)) 

1584 >>> np.allclose(b, np.matmul(u[..., :3] * s[..., None, :], vh)) 

1585 True 

1586 >>> np.allclose(b, np.matmul(u[..., :3], s[..., None] * vh)) 

1587 True 

1588 

1589 Reconstruction based on reduced SVD, 4D case: 

1590 

1591 >>> u, s, vh = np.linalg.svd(b, full_matrices=False) 

1592 >>> u.shape, s.shape, vh.shape 

1593 ((2, 7, 8, 3), (2, 7, 3), (2, 7, 3, 3)) 

1594 >>> np.allclose(b, np.matmul(u * s[..., None, :], vh)) 

1595 True 

1596 >>> np.allclose(b, np.matmul(u, s[..., None] * vh)) 

1597 True 

1598 

1599 """ 

1600 import numpy as _nx 

1601 a, wrap = _makearray(a) 

1602 

1603 if hermitian: 

1604 # note: lapack svd returns eigenvalues with s ** 2 sorted descending, 

1605 # but eig returns s sorted ascending, so we re-order the eigenvalues 

1606 # and related arrays to have the correct order 

1607 if compute_uv: 

1608 s, u = eigh(a) 

1609 sgn = sign(s) 

1610 s = abs(s) 

1611 sidx = argsort(s)[..., ::-1] 

1612 sgn = _nx.take_along_axis(sgn, sidx, axis=-1) 

1613 s = _nx.take_along_axis(s, sidx, axis=-1) 

1614 u = _nx.take_along_axis(u, sidx[..., None, :], axis=-1) 

1615 # singular values are unsigned, move the sign into v 

1616 vt = transpose(u * sgn[..., None, :]).conjugate() 

1617 return wrap(u), s, wrap(vt) 

1618 else: 

1619 s = eigvalsh(a) 

1620 s = abs(s) 

1621 return sort(s)[..., ::-1] 

1622 

1623 _assert_stacked_2d(a) 

1624 t, result_t = _commonType(a) 

1625 

1626 extobj = get_linalg_error_extobj(_raise_linalgerror_svd_nonconvergence) 

1627 

1628 m, n = a.shape[-2:] 

1629 if compute_uv: 

1630 if full_matrices: 

1631 if m < n: 

1632 gufunc = _umath_linalg.svd_m_f 

1633 else: 

1634 gufunc = _umath_linalg.svd_n_f 

1635 else: 

1636 if m < n: 

1637 gufunc = _umath_linalg.svd_m_s 

1638 else: 

1639 gufunc = _umath_linalg.svd_n_s 

1640 

1641 signature = 'D->DdD' if isComplexType(t) else 'd->ddd' 

1642 u, s, vh = gufunc(a, signature=signature, extobj=extobj) 

1643 u = u.astype(result_t, copy=False) 

1644 s = s.astype(_realType(result_t), copy=False) 

1645 vh = vh.astype(result_t, copy=False) 

1646 return wrap(u), s, wrap(vh) 

1647 else: 

1648 if m < n: 

1649 gufunc = _umath_linalg.svd_m 

1650 else: 

1651 gufunc = _umath_linalg.svd_n 

1652 

1653 signature = 'D->d' if isComplexType(t) else 'd->d' 

1654 s = gufunc(a, signature=signature, extobj=extobj) 

1655 s = s.astype(_realType(result_t), copy=False) 

1656 return s 

1657 

1658 

1659def _cond_dispatcher(x, p=None): 

1660 return (x,) 

1661 

1662 

1663@array_function_dispatch(_cond_dispatcher) 

1664def cond(x, p=None): 

1665 """ 

1666 Compute the condition number of a matrix. 

1667 

1668 This function is capable of returning the condition number using 

1669 one of seven different norms, depending on the value of `p` (see 

1670 Parameters below). 

1671 

1672 Parameters 

1673 ---------- 

1674 x : (..., M, N) array_like 

1675 The matrix whose condition number is sought. 

1676 p : {None, 1, -1, 2, -2, inf, -inf, 'fro'}, optional 

1677 Order of the norm used in the condition number computation: 

1678 

1679 ===== ============================ 

1680 p norm for matrices 

1681 ===== ============================ 

1682 None 2-norm, computed directly using the ``SVD`` 

1683 'fro' Frobenius norm 

1684 inf max(sum(abs(x), axis=1)) 

1685 -inf min(sum(abs(x), axis=1)) 

1686 1 max(sum(abs(x), axis=0)) 

1687 -1 min(sum(abs(x), axis=0)) 

1688 2 2-norm (largest sing. value) 

1689 -2 smallest singular value 

1690 ===== ============================ 

1691 

1692 inf means the `numpy.inf` object, and the Frobenius norm is 

1693 the root-of-sum-of-squares norm. 

1694 

1695 Returns 

1696 ------- 

1697 c : {float, inf} 

1698 The condition number of the matrix. May be infinite. 

1699 

1700 See Also 

1701 -------- 

1702 numpy.linalg.norm 

1703 

1704 Notes 

1705 ----- 

1706 The condition number of `x` is defined as the norm of `x` times the 

1707 norm of the inverse of `x` [1]_; the norm can be the usual L2-norm 

1708 (root-of-sum-of-squares) or one of a number of other matrix norms. 

1709 

1710 References 

1711 ---------- 

1712 .. [1] G. Strang, *Linear Algebra and Its Applications*, Orlando, FL, 

1713 Academic Press, Inc., 1980, pg. 285. 

1714 

1715 Examples 

1716 -------- 

1717 >>> from numpy import linalg as LA 

1718 >>> a = np.array([[1, 0, -1], [0, 1, 0], [1, 0, 1]]) 

1719 >>> a 

1720 array([[ 1, 0, -1], 

1721 [ 0, 1, 0], 

1722 [ 1, 0, 1]]) 

1723 >>> LA.cond(a) 

1724 1.4142135623730951 

1725 >>> LA.cond(a, 'fro') 

1726 3.1622776601683795 

1727 >>> LA.cond(a, np.inf) 

1728 2.0 

1729 >>> LA.cond(a, -np.inf) 

1730 1.0 

1731 >>> LA.cond(a, 1) 

1732 2.0 

1733 >>> LA.cond(a, -1) 

1734 1.0 

1735 >>> LA.cond(a, 2) 

1736 1.4142135623730951 

1737 >>> LA.cond(a, -2) 

1738 0.70710678118654746 # may vary 

1739 >>> min(LA.svd(a, compute_uv=False))*min(LA.svd(LA.inv(a), compute_uv=False)) 

1740 0.70710678118654746 # may vary 

1741 

1742 """ 

1743 x = asarray(x) # in case we have a matrix 

1744 if _is_empty_2d(x): 

1745 raise LinAlgError("cond is not defined on empty arrays") 

1746 if p is None or p == 2 or p == -2: 

1747 s = svd(x, compute_uv=False) 

1748 with errstate(all='ignore'): 

1749 if p == -2: 

1750 r = s[..., -1] / s[..., 0] 

1751 else: 

1752 r = s[..., 0] / s[..., -1] 

1753 else: 

1754 # Call inv(x) ignoring errors. The result array will 

1755 # contain nans in the entries where inversion failed. 

1756 _assert_stacked_2d(x) 

1757 _assert_stacked_square(x) 

1758 t, result_t = _commonType(x) 

1759 signature = 'D->D' if isComplexType(t) else 'd->d' 

1760 with errstate(all='ignore'): 

1761 invx = _umath_linalg.inv(x, signature=signature) 

1762 r = norm(x, p, axis=(-2, -1)) * norm(invx, p, axis=(-2, -1)) 

1763 r = r.astype(result_t, copy=False) 

1764 

1765 # Convert nans to infs unless the original array had nan entries 

1766 r = asarray(r) 

1767 nan_mask = isnan(r) 

1768 if nan_mask.any(): 

1769 nan_mask &= ~isnan(x).any(axis=(-2, -1)) 

1770 if r.ndim > 0: 

1771 r[nan_mask] = Inf 

1772 elif nan_mask: 

1773 r[()] = Inf 

1774 

1775 # Convention is to return scalars instead of 0d arrays 

1776 if r.ndim == 0: 

1777 r = r[()] 

1778 

1779 return r 

1780 

1781 

1782def _matrix_rank_dispatcher(A, tol=None, hermitian=None): 

1783 return (A,) 

1784 

1785 

1786@array_function_dispatch(_matrix_rank_dispatcher) 

1787def matrix_rank(A, tol=None, hermitian=False): 

1788 """ 

1789 Return matrix rank of array using SVD method 

1790 

1791 Rank of the array is the number of singular values of the array that are 

1792 greater than `tol`. 

1793 

1794 .. versionchanged:: 1.14 

1795 Can now operate on stacks of matrices 

1796 

1797 Parameters 

1798 ---------- 

1799 A : {(M,), (..., M, N)} array_like 

1800 Input vector or stack of matrices. 

1801 tol : (...) array_like, float, optional 

1802 Threshold below which SVD values are considered zero. If `tol` is 

1803 None, and ``S`` is an array with singular values for `M`, and 

1804 ``eps`` is the epsilon value for datatype of ``S``, then `tol` is 

1805 set to ``S.max() * max(M, N) * eps``. 

1806 

1807 .. versionchanged:: 1.14 

1808 Broadcasted against the stack of matrices 

1809 hermitian : bool, optional 

1810 If True, `A` is assumed to be Hermitian (symmetric if real-valued), 

1811 enabling a more efficient method for finding singular values. 

1812 Defaults to False. 

1813 

1814 .. versionadded:: 1.14 

1815 

1816 Returns 

1817 ------- 

1818 rank : (...) array_like 

1819 Rank of A. 

1820 

1821 Notes 

1822 ----- 

1823 The default threshold to detect rank deficiency is a test on the magnitude 

1824 of the singular values of `A`. By default, we identify singular values less 

1825 than ``S.max() * max(M, N) * eps`` as indicating rank deficiency (with 

1826 the symbols defined above). This is the algorithm MATLAB uses [1]. It also 

1827 appears in *Numerical recipes* in the discussion of SVD solutions for linear 

1828 least squares [2]. 

1829 

1830 This default threshold is designed to detect rank deficiency accounting for 

1831 the numerical errors of the SVD computation. Imagine that there is a column 

1832 in `A` that is an exact (in floating point) linear combination of other 

1833 columns in `A`. Computing the SVD on `A` will not produce a singular value 

1834 exactly equal to 0 in general: any difference of the smallest SVD value from 

1835 0 will be caused by numerical imprecision in the calculation of the SVD. 

1836 Our threshold for small SVD values takes this numerical imprecision into 

1837 account, and the default threshold will detect such numerical rank 

1838 deficiency. The threshold may declare a matrix `A` rank deficient even if 

1839 the linear combination of some columns of `A` is not exactly equal to 

1840 another column of `A` but only numerically very close to another column of 

1841 `A`. 

1842 

1843 We chose our default threshold because it is in wide use. Other thresholds 

1844 are possible. For example, elsewhere in the 2007 edition of *Numerical 

1845 recipes* there is an alternative threshold of ``S.max() * 

1846 np.finfo(A.dtype).eps / 2. * np.sqrt(m + n + 1.)``. The authors describe 

1847 this threshold as being based on "expected roundoff error" (p 71). 

1848 

1849 The thresholds above deal with floating point roundoff error in the 

1850 calculation of the SVD. However, you may have more information about the 

1851 sources of error in `A` that would make you consider other tolerance values 

1852 to detect *effective* rank deficiency. The most useful measure of the 

1853 tolerance depends on the operations you intend to use on your matrix. For 

1854 example, if your data come from uncertain measurements with uncertainties 

1855 greater than floating point epsilon, choosing a tolerance near that 

1856 uncertainty may be preferable. The tolerance may be absolute if the 

1857 uncertainties are absolute rather than relative. 

1858 

1859 References 

1860 ---------- 

1861 .. [1] MATLAB reference documentation, "Rank" 

1862 https://www.mathworks.com/help/techdoc/ref/rank.html 

1863 .. [2] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, 

1864 "Numerical Recipes (3rd edition)", Cambridge University Press, 2007, 

1865 page 795. 

1866 

1867 Examples 

1868 -------- 

1869 >>> from numpy.linalg import matrix_rank 

1870 >>> matrix_rank(np.eye(4)) # Full rank matrix 

1871 4 

1872 >>> I=np.eye(4); I[-1,-1] = 0. # rank deficient matrix 

1873 >>> matrix_rank(I) 

1874 3 

1875 >>> matrix_rank(np.ones((4,))) # 1 dimension - rank 1 unless all 0 

1876 1 

1877 >>> matrix_rank(np.zeros((4,))) 

1878 0 

1879 """ 

1880 A = asarray(A) 

1881 if A.ndim < 2: 

1882 return int(not all(A==0)) 

1883 S = svd(A, compute_uv=False, hermitian=hermitian) 

1884 if tol is None: 

1885 tol = S.max(axis=-1, keepdims=True) * max(A.shape[-2:]) * finfo(S.dtype).eps 

1886 else: 

1887 tol = asarray(tol)[..., newaxis] 

1888 return count_nonzero(S > tol, axis=-1) 

1889 

1890 

1891# Generalized inverse 

1892 

1893def _pinv_dispatcher(a, rcond=None, hermitian=None): 

1894 return (a,) 

1895 

1896 

1897@array_function_dispatch(_pinv_dispatcher) 

1898def pinv(a, rcond=1e-15, hermitian=False): 

1899 """ 

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

1901 

1902 Calculate the generalized inverse of a matrix using its 

1903 singular-value decomposition (SVD) and including all 

1904 *large* singular values. 

1905 

1906 .. versionchanged:: 1.14 

1907 Can now operate on stacks of matrices 

1908 

1909 Parameters 

1910 ---------- 

1911 a : (..., M, N) array_like 

1912 Matrix or stack of matrices to be pseudo-inverted. 

1913 rcond : (...) array_like of float 

1914 Cutoff for small singular values. 

1915 Singular values less than or equal to 

1916 ``rcond * largest_singular_value`` are set to zero. 

1917 Broadcasts against the stack of matrices. 

1918 hermitian : bool, optional 

1919 If True, `a` is assumed to be Hermitian (symmetric if real-valued), 

1920 enabling a more efficient method for finding singular values. 

1921 Defaults to False. 

1922 

1923 .. versionadded:: 1.17.0 

1924 

1925 Returns 

1926 ------- 

1927 B : (..., N, M) ndarray 

1928 The pseudo-inverse of `a`. If `a` is a `matrix` instance, then so 

1929 is `B`. 

1930 

1931 Raises 

1932 ------ 

1933 LinAlgError 

1934 If the SVD computation does not converge. 

1935 

1936 See Also 

1937 -------- 

1938 scipy.linalg.pinv : Similar function in SciPy. 

1939 scipy.linalg.pinvh : Compute the (Moore-Penrose) pseudo-inverse of a 

1940 Hermitian matrix. 

1941 

1942 Notes 

1943 ----- 

1944 The pseudo-inverse of a matrix A, denoted :math:`A^+`, is 

1945 defined as: "the matrix that 'solves' [the least-squares problem] 

1946 :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then 

1947 :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`. 

1948 

1949 It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular 

1950 value decomposition of A, then 

1951 :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are 

1952 orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting 

1953 of A's so-called singular values, (followed, typically, by 

1954 zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix 

1955 consisting of the reciprocals of A's singular values 

1956 (again, followed by zeros). [1]_ 

1957 

1958 References 

1959 ---------- 

1960 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, 

1961 FL, Academic Press, Inc., 1980, pp. 139-142. 

1962 

1963 Examples 

1964 -------- 

1965 The following example checks that ``a * a+ * a == a`` and 

1966 ``a+ * a * a+ == a+``: 

1967 

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

1969 >>> B = np.linalg.pinv(a) 

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

1971 True 

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

1973 True 

1974 

1975 """ 

1976 a, wrap = _makearray(a) 

1977 rcond = asarray(rcond) 

1978 if _is_empty_2d(a): 

1979 m, n = a.shape[-2:] 

1980 res = empty(a.shape[:-2] + (n, m), dtype=a.dtype) 

1981 return wrap(res) 

1982 a = a.conjugate() 

1983 u, s, vt = svd(a, full_matrices=False, hermitian=hermitian) 

1984 

1985 # discard small singular values 

1986 cutoff = rcond[..., newaxis] * amax(s, axis=-1, keepdims=True) 

1987 large = s > cutoff 

1988 s = divide(1, s, where=large, out=s) 

1989 s[~large] = 0 

1990 

1991 res = matmul(transpose(vt), multiply(s[..., newaxis], transpose(u))) 

1992 return wrap(res) 

1993 

1994 

1995# Determinant 

1996 

1997 

1998@array_function_dispatch(_unary_dispatcher) 

1999def slogdet(a): 

2000 """ 

2001 Compute the sign and (natural) logarithm of the determinant of an array. 

2002 

2003 If an array has a very small or very large determinant, then a call to 

2004 `det` may overflow or underflow. This routine is more robust against such 

2005 issues, because it computes the logarithm of the determinant rather than 

2006 the determinant itself. 

2007 

2008 Parameters 

2009 ---------- 

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

2011 Input array, has to be a square 2-D array. 

2012 

2013 Returns 

2014 ------- 

2015 sign : (...) array_like 

2016 A number representing the sign of the determinant. For a real matrix, 

2017 this is 1, 0, or -1. For a complex matrix, this is a complex number 

2018 with absolute value 1 (i.e., it is on the unit circle), or else 0. 

2019 logdet : (...) array_like 

2020 The natural log of the absolute value of the determinant. 

2021 

2022 If the determinant is zero, then `sign` will be 0 and `logdet` will be 

2023 -Inf. In all cases, the determinant is equal to ``sign * np.exp(logdet)``. 

2024 

2025 See Also 

2026 -------- 

2027 det 

2028 

2029 Notes 

2030 ----- 

2031 

2032 .. versionadded:: 1.8.0 

2033 

2034 Broadcasting rules apply, see the `numpy.linalg` documentation for 

2035 details. 

2036 

2037 .. versionadded:: 1.6.0 

2038 

2039 The determinant is computed via LU factorization using the LAPACK 

2040 routine ``z/dgetrf``. 

2041 

2042 

2043 Examples 

2044 -------- 

2045 The determinant of a 2-D array ``[[a, b], [c, d]]`` is ``ad - bc``: 

2046 

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

2048 >>> (sign, logdet) = np.linalg.slogdet(a) 

2049 >>> (sign, logdet) 

2050 (-1, 0.69314718055994529) # may vary 

2051 >>> sign * np.exp(logdet) 

2052 -2.0 

2053 

2054 Computing log-determinants for a stack of matrices: 

2055 

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

2057 >>> a.shape 

2058 (3, 2, 2) 

2059 >>> sign, logdet = np.linalg.slogdet(a) 

2060 >>> (sign, logdet) 

2061 (array([-1., -1., -1.]), array([ 0.69314718, 1.09861229, 2.07944154])) 

2062 >>> sign * np.exp(logdet) 

2063 array([-2., -3., -8.]) 

2064 

2065 This routine succeeds where ordinary `det` does not: 

2066 

2067 >>> np.linalg.det(np.eye(500) * 0.1) 

2068 0.0 

2069 >>> np.linalg.slogdet(np.eye(500) * 0.1) 

2070 (1, -1151.2925464970228) 

2071 

2072 """ 

2073 a = asarray(a) 

2074 _assert_stacked_2d(a) 

2075 _assert_stacked_square(a) 

2076 t, result_t = _commonType(a) 

2077 real_t = _realType(result_t) 

2078 signature = 'D->Dd' if isComplexType(t) else 'd->dd' 

2079 sign, logdet = _umath_linalg.slogdet(a, signature=signature) 

2080 sign = sign.astype(result_t, copy=False) 

2081 logdet = logdet.astype(real_t, copy=False) 

2082 return sign, logdet 

2083 

2084 

2085@array_function_dispatch(_unary_dispatcher) 

2086def det(a): 

2087 """ 

2088 Compute the determinant of an array. 

2089 

2090 Parameters 

2091 ---------- 

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

2093 Input array to compute determinants for. 

2094 

2095 Returns 

2096 ------- 

2097 det : (...) array_like 

2098 Determinant of `a`. 

2099 

2100 See Also 

2101 -------- 

2102 slogdet : Another way to represent the determinant, more suitable 

2103 for large matrices where underflow/overflow may occur. 

2104 scipy.linalg.det : Similar function in SciPy. 

2105 

2106 Notes 

2107 ----- 

2108 

2109 .. versionadded:: 1.8.0 

2110 

2111 Broadcasting rules apply, see the `numpy.linalg` documentation for 

2112 details. 

2113 

2114 The determinant is computed via LU factorization using the LAPACK 

2115 routine ``z/dgetrf``. 

2116 

2117 Examples 

2118 -------- 

2119 The determinant of a 2-D array [[a, b], [c, d]] is ad - bc: 

2120 

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

2122 >>> np.linalg.det(a) 

2123 -2.0 # may vary 

2124 

2125 Computing determinants for a stack of matrices: 

2126 

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

2128 >>> a.shape 

2129 (3, 2, 2) 

2130 >>> np.linalg.det(a) 

2131 array([-2., -3., -8.]) 

2132 

2133 """ 

2134 a = asarray(a) 

2135 _assert_stacked_2d(a) 

2136 _assert_stacked_square(a) 

2137 t, result_t = _commonType(a) 

2138 signature = 'D->D' if isComplexType(t) else 'd->d' 

2139 r = _umath_linalg.det(a, signature=signature) 

2140 r = r.astype(result_t, copy=False) 

2141 return r 

2142 

2143 

2144# Linear Least Squares 

2145 

2146def _lstsq_dispatcher(a, b, rcond=None): 

2147 return (a, b) 

2148 

2149 

2150@array_function_dispatch(_lstsq_dispatcher) 

2151def lstsq(a, b, rcond="warn"): 

2152 r""" 

2153 Return the least-squares solution to a linear matrix equation. 

2154 

2155 Computes the vector `x` that approximately solves the equation 

2156 ``a @ x = b``. The equation may be under-, well-, or over-determined 

2157 (i.e., the number of linearly independent rows of `a` can be less than, 

2158 equal to, or greater than its number of linearly independent columns). 

2159 If `a` is square and of full rank, then `x` (but for round-off error) 

2160 is the "exact" solution of the equation. Else, `x` minimizes the 

2161 Euclidean 2-norm :math:`||b - ax||`. If there are multiple minimizing 

2162 solutions, the one with the smallest 2-norm :math:`||x||` is returned. 

2163 

2164 Parameters 

2165 ---------- 

2166 a : (M, N) array_like 

2167 "Coefficient" matrix. 

2168 b : {(M,), (M, K)} array_like 

2169 Ordinate or "dependent variable" values. If `b` is two-dimensional, 

2170 the least-squares solution is calculated for each of the `K` columns 

2171 of `b`. 

2172 rcond : float, optional 

2173 Cut-off ratio for small singular values of `a`. 

2174 For the purposes of rank determination, singular values are treated 

2175 as zero if they are smaller than `rcond` times the largest singular 

2176 value of `a`. 

2177 

2178 .. versionchanged:: 1.14.0 

2179 If not set, a FutureWarning is given. The previous default 

2180 of ``-1`` will use the machine precision as `rcond` parameter, 

2181 the new default will use the machine precision times `max(M, N)`. 

2182 To silence the warning and use the new default, use ``rcond=None``, 

2183 to keep using the old behavior, use ``rcond=-1``. 

2184 

2185 Returns 

2186 ------- 

2187 x : {(N,), (N, K)} ndarray 

2188 Least-squares solution. If `b` is two-dimensional, 

2189 the solutions are in the `K` columns of `x`. 

2190 residuals : {(1,), (K,), (0,)} ndarray 

2191 Sums of squared residuals: Squared Euclidean 2-norm for each column in 

2192 ``b - a @ x``. 

2193 If the rank of `a` is < N or M <= N, this is an empty array. 

2194 If `b` is 1-dimensional, this is a (1,) shape array. 

2195 Otherwise the shape is (K,). 

2196 rank : int 

2197 Rank of matrix `a`. 

2198 s : (min(M, N),) ndarray 

2199 Singular values of `a`. 

2200 

2201 Raises 

2202 ------ 

2203 LinAlgError 

2204 If computation does not converge. 

2205 

2206 See Also 

2207 -------- 

2208 scipy.linalg.lstsq : Similar function in SciPy. 

2209 

2210 Notes 

2211 ----- 

2212 If `b` is a matrix, then all array results are returned as matrices. 

2213 

2214 Examples 

2215 -------- 

2216 Fit a line, ``y = mx + c``, through some noisy data-points: 

2217 

2218 >>> x = np.array([0, 1, 2, 3]) 

2219 >>> y = np.array([-1, 0.2, 0.9, 2.1]) 

2220 

2221 By examining the coefficients, we see that the line should have a 

2222 gradient of roughly 1 and cut the y-axis at, more or less, -1. 

2223 

2224 We can rewrite the line equation as ``y = Ap``, where ``A = [[x 1]]`` 

2225 and ``p = [[m], [c]]``. Now use `lstsq` to solve for `p`: 

2226 

2227 >>> A = np.vstack([x, np.ones(len(x))]).T 

2228 >>> A 

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

2230 [ 1., 1.], 

2231 [ 2., 1.], 

2232 [ 3., 1.]]) 

2233 

2234 >>> m, c = np.linalg.lstsq(A, y, rcond=None)[0] 

2235 >>> m, c 

2236 (1.0 -0.95) # may vary 

2237 

2238 Plot the data along with the fitted line: 

2239 

2240 >>> import matplotlib.pyplot as plt 

2241 >>> _ = plt.plot(x, y, 'o', label='Original data', markersize=10) 

2242 >>> _ = plt.plot(x, m*x + c, 'r', label='Fitted line') 

2243 >>> _ = plt.legend() 

2244 >>> plt.show() 

2245 

2246 """ 

2247 a, _ = _makearray(a) 

2248 b, wrap = _makearray(b) 

2249 is_1d = b.ndim == 1 

2250 if is_1d: 

2251 b = b[:, newaxis] 

2252 _assert_2d(a, b) 

2253 m, n = a.shape[-2:] 

2254 m2, n_rhs = b.shape[-2:] 

2255 if m != m2: 

2256 raise LinAlgError('Incompatible dimensions') 

2257 

2258 t, result_t = _commonType(a, b) 

2259 result_real_t = _realType(result_t) 

2260 

2261 # Determine default rcond value 

2262 if rcond == "warn": 

2263 # 2017-08-19, 1.14.0 

2264 warnings.warn("`rcond` parameter will change to the default of " 

2265 "machine precision times ``max(M, N)`` where M and N " 

2266 "are the input matrix dimensions.\n" 

2267 "To use the future default and silence this warning " 

2268 "we advise to pass `rcond=None`, to keep using the old, " 

2269 "explicitly pass `rcond=-1`.", 

2270 FutureWarning, stacklevel=3) 

2271 rcond = -1 

2272 if rcond is None: 

2273 rcond = finfo(t).eps * max(n, m) 

2274 

2275 if m <= n: 

2276 gufunc = _umath_linalg.lstsq_m 

2277 else: 

2278 gufunc = _umath_linalg.lstsq_n 

2279 

2280 signature = 'DDd->Ddid' if isComplexType(t) else 'ddd->ddid' 

2281 extobj = get_linalg_error_extobj(_raise_linalgerror_lstsq) 

2282 if n_rhs == 0: 

2283 # lapack can't handle n_rhs = 0 - so allocate the array one larger in that axis 

2284 b = zeros(b.shape[:-2] + (m, n_rhs + 1), dtype=b.dtype) 

2285 x, resids, rank, s = gufunc(a, b, rcond, signature=signature, extobj=extobj) 

2286 if m == 0: 

2287 x[...] = 0 

2288 if n_rhs == 0: 

2289 # remove the item we added 

2290 x = x[..., :n_rhs] 

2291 resids = resids[..., :n_rhs] 

2292 

2293 # remove the axis we added 

2294 if is_1d: 

2295 x = x.squeeze(axis=-1) 

2296 # we probably should squeeze resids too, but we can't 

2297 # without breaking compatibility. 

2298 

2299 # as documented 

2300 if rank != n or m <= n: 

2301 resids = array([], result_real_t) 

2302 

2303 # coerce output arrays 

2304 s = s.astype(result_real_t, copy=False) 

2305 resids = resids.astype(result_real_t, copy=False) 

2306 x = x.astype(result_t, copy=True) # Copying lets the memory in r_parts be freed 

2307 return wrap(x), wrap(resids), rank, s 

2308 

2309 

2310def _multi_svd_norm(x, row_axis, col_axis, op): 

2311 """Compute a function of the singular values of the 2-D matrices in `x`. 

2312 

2313 This is a private utility function used by `numpy.linalg.norm()`. 

2314 

2315 Parameters 

2316 ---------- 

2317 x : ndarray 

2318 row_axis, col_axis : int 

2319 The axes of `x` that hold the 2-D matrices. 

2320 op : callable 

2321 This should be either numpy.amin or `numpy.amax` or `numpy.sum`. 

2322 

2323 Returns 

2324 ------- 

2325 result : float or ndarray 

2326 If `x` is 2-D, the return values is a float. 

2327 Otherwise, it is an array with ``x.ndim - 2`` dimensions. 

2328 The return values are either the minimum or maximum or sum of the 

2329 singular values of the matrices, depending on whether `op` 

2330 is `numpy.amin` or `numpy.amax` or `numpy.sum`. 

2331 

2332 """ 

2333 y = moveaxis(x, (row_axis, col_axis), (-2, -1)) 

2334 result = op(svd(y, compute_uv=False), axis=-1) 

2335 return result 

2336 

2337 

2338def _norm_dispatcher(x, ord=None, axis=None, keepdims=None): 

2339 return (x,) 

2340 

2341 

2342@array_function_dispatch(_norm_dispatcher) 

2343def norm(x, ord=None, axis=None, keepdims=False): 

2344 """ 

2345 Matrix or vector norm. 

2346 

2347 This function is able to return one of eight different matrix norms, 

2348 or one of an infinite number of vector norms (described below), depending 

2349 on the value of the ``ord`` parameter. 

2350 

2351 Parameters 

2352 ---------- 

2353 x : array_like 

2354 Input array. If `axis` is None, `x` must be 1-D or 2-D, unless `ord` 

2355 is None. If both `axis` and `ord` are None, the 2-norm of 

2356 ``x.ravel`` will be returned. 

2357 ord : {non-zero int, inf, -inf, 'fro', 'nuc'}, optional 

2358 Order of the norm (see table under ``Notes``). inf means numpy's 

2359 `inf` object. The default is None. 

2360 axis : {None, int, 2-tuple of ints}, optional. 

2361 If `axis` is an integer, it specifies the axis of `x` along which to 

2362 compute the vector norms. If `axis` is a 2-tuple, it specifies the 

2363 axes that hold 2-D matrices, and the matrix norms of these matrices 

2364 are computed. If `axis` is None then either a vector norm (when `x` 

2365 is 1-D) or a matrix norm (when `x` is 2-D) is returned. The default 

2366 is None. 

2367 

2368 .. versionadded:: 1.8.0 

2369 

2370 keepdims : bool, optional 

2371 If this is set to True, the axes which are normed over are left in the 

2372 result as dimensions with size one. With this option the result will 

2373 broadcast correctly against the original `x`. 

2374 

2375 .. versionadded:: 1.10.0 

2376 

2377 Returns 

2378 ------- 

2379 n : float or ndarray 

2380 Norm of the matrix or vector(s). 

2381 

2382 See Also 

2383 -------- 

2384 scipy.linalg.norm : Similar function in SciPy. 

2385 

2386 Notes 

2387 ----- 

2388 For values of ``ord < 1``, the result is, strictly speaking, not a 

2389 mathematical 'norm', but it may still be useful for various numerical 

2390 purposes. 

2391 

2392 The following norms can be calculated: 

2393 

2394 ===== ============================ ========================== 

2395 ord norm for matrices norm for vectors 

2396 ===== ============================ ========================== 

2397 None Frobenius norm 2-norm 

2398 'fro' Frobenius norm -- 

2399 'nuc' nuclear norm -- 

2400 inf max(sum(abs(x), axis=1)) max(abs(x)) 

2401 -inf min(sum(abs(x), axis=1)) min(abs(x)) 

2402 0 -- sum(x != 0) 

2403 1 max(sum(abs(x), axis=0)) as below 

2404 -1 min(sum(abs(x), axis=0)) as below 

2405 2 2-norm (largest sing. value) as below 

2406 -2 smallest singular value as below 

2407 other -- sum(abs(x)**ord)**(1./ord) 

2408 ===== ============================ ========================== 

2409 

2410 The Frobenius norm is given by [1]_: 

2411 

2412 :math:`||A||_F = [\\sum_{i,j} abs(a_{i,j})^2]^{1/2}` 

2413 

2414 The nuclear norm is the sum of the singular values. 

2415 

2416 Both the Frobenius and nuclear norm orders are only defined for 

2417 matrices and raise a ValueError when ``x.ndim != 2``. 

2418 

2419 References 

2420 ---------- 

2421 .. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*, 

2422 Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15 

2423 

2424 Examples 

2425 -------- 

2426 >>> from numpy import linalg as LA 

2427 >>> a = np.arange(9) - 4 

2428 >>> a 

2429 array([-4, -3, -2, ..., 2, 3, 4]) 

2430 >>> b = a.reshape((3, 3)) 

2431 >>> b 

2432 array([[-4, -3, -2], 

2433 [-1, 0, 1], 

2434 [ 2, 3, 4]]) 

2435 

2436 >>> LA.norm(a) 

2437 7.745966692414834 

2438 >>> LA.norm(b) 

2439 7.745966692414834 

2440 >>> LA.norm(b, 'fro') 

2441 7.745966692414834 

2442 >>> LA.norm(a, np.inf) 

2443 4.0 

2444 >>> LA.norm(b, np.inf) 

2445 9.0 

2446 >>> LA.norm(a, -np.inf) 

2447 0.0 

2448 >>> LA.norm(b, -np.inf) 

2449 2.0 

2450 

2451 >>> LA.norm(a, 1) 

2452 20.0 

2453 >>> LA.norm(b, 1) 

2454 7.0 

2455 >>> LA.norm(a, -1) 

2456 -4.6566128774142013e-010 

2457 >>> LA.norm(b, -1) 

2458 6.0 

2459 >>> LA.norm(a, 2) 

2460 7.745966692414834 

2461 >>> LA.norm(b, 2) 

2462 7.3484692283495345 

2463 

2464 >>> LA.norm(a, -2) 

2465 0.0 

2466 >>> LA.norm(b, -2) 

2467 1.8570331885190563e-016 # may vary 

2468 >>> LA.norm(a, 3) 

2469 5.8480354764257312 # may vary 

2470 >>> LA.norm(a, -3) 

2471 0.0 

2472 

2473 Using the `axis` argument to compute vector norms: 

2474 

2475 >>> c = np.array([[ 1, 2, 3], 

2476 ... [-1, 1, 4]]) 

2477 >>> LA.norm(c, axis=0) 

2478 array([ 1.41421356, 2.23606798, 5. ]) 

2479 >>> LA.norm(c, axis=1) 

2480 array([ 3.74165739, 4.24264069]) 

2481 >>> LA.norm(c, ord=1, axis=1) 

2482 array([ 6., 6.]) 

2483 

2484 Using the `axis` argument to compute matrix norms: 

2485 

2486 >>> m = np.arange(8).reshape(2,2,2) 

2487 >>> LA.norm(m, axis=(1,2)) 

2488 array([ 3.74165739, 11.22497216]) 

2489 >>> LA.norm(m[0, :, :]), LA.norm(m[1, :, :]) 

2490 (3.7416573867739413, 11.224972160321824) 

2491 

2492 """ 

2493 x = asarray(x) 

2494 

2495 if not issubclass(x.dtype.type, (inexact, object_)): 

2496 x = x.astype(float) 

2497 

2498 # Immediately handle some default, simple, fast, and common cases. 

2499 if axis is None: 

2500 ndim = x.ndim 

2501 if ((ord is None) or 

2502 (ord in ('f', 'fro') and ndim == 2) or 

2503 (ord == 2 and ndim == 1)): 

2504 

2505 x = x.ravel(order='K') 

2506 if isComplexType(x.dtype.type): 

2507 x_real = x.real 

2508 x_imag = x.imag 

2509 sqnorm = x_real.dot(x_real) + x_imag.dot(x_imag) 

2510 else: 

2511 sqnorm = x.dot(x) 

2512 ret = sqrt(sqnorm) 

2513 if keepdims: 

2514 ret = ret.reshape(ndim*[1]) 

2515 return ret 

2516 

2517 # Normalize the `axis` argument to a tuple. 

2518 nd = x.ndim 

2519 if axis is None: 

2520 axis = tuple(range(nd)) 

2521 elif not isinstance(axis, tuple): 

2522 try: 

2523 axis = int(axis) 

2524 except Exception as e: 

2525 raise TypeError("'axis' must be None, an integer or a tuple of integers") from e 

2526 axis = (axis,) 

2527 

2528 if len(axis) == 1: 

2529 if ord == Inf: 

2530 return abs(x).max(axis=axis, keepdims=keepdims) 

2531 elif ord == -Inf: 

2532 return abs(x).min(axis=axis, keepdims=keepdims) 

2533 elif ord == 0: 

2534 # Zero norm 

2535 return (x != 0).astype(x.real.dtype).sum(axis=axis, keepdims=keepdims) 

2536 elif ord == 1: 

2537 # special case for speedup 

2538 return add.reduce(abs(x), axis=axis, keepdims=keepdims) 

2539 elif ord is None or ord == 2: 

2540 # special case for speedup 

2541 s = (x.conj() * x).real 

2542 return sqrt(add.reduce(s, axis=axis, keepdims=keepdims)) 

2543 # None of the str-type keywords for ord ('fro', 'nuc') 

2544 # are valid for vectors 

2545 elif isinstance(ord, str): 

2546 raise ValueError(f"Invalid norm order '{ord}' for vectors") 

2547 else: 

2548 absx = abs(x) 

2549 absx **= ord 

2550 ret = add.reduce(absx, axis=axis, keepdims=keepdims) 

2551 ret **= reciprocal(ord, dtype=ret.dtype) 

2552 return ret 

2553 elif len(axis) == 2: 

2554 row_axis, col_axis = axis 

2555 row_axis = normalize_axis_index(row_axis, nd) 

2556 col_axis = normalize_axis_index(col_axis, nd) 

2557 if row_axis == col_axis: 

2558 raise ValueError('Duplicate axes given.') 

2559 if ord == 2: 

2560 ret = _multi_svd_norm(x, row_axis, col_axis, amax) 

2561 elif ord == -2: 

2562 ret = _multi_svd_norm(x, row_axis, col_axis, amin) 

2563 elif ord == 1: 

2564 if col_axis > row_axis: 

2565 col_axis -= 1 

2566 ret = add.reduce(abs(x), axis=row_axis).max(axis=col_axis) 

2567 elif ord == Inf: 

2568 if row_axis > col_axis: 

2569 row_axis -= 1 

2570 ret = add.reduce(abs(x), axis=col_axis).max(axis=row_axis) 

2571 elif ord == -1: 

2572 if col_axis > row_axis: 

2573 col_axis -= 1 

2574 ret = add.reduce(abs(x), axis=row_axis).min(axis=col_axis) 

2575 elif ord == -Inf: 

2576 if row_axis > col_axis: 

2577 row_axis -= 1 

2578 ret = add.reduce(abs(x), axis=col_axis).min(axis=row_axis) 

2579 elif ord in [None, 'fro', 'f']: 

2580 ret = sqrt(add.reduce((x.conj() * x).real, axis=axis)) 

2581 elif ord == 'nuc': 

2582 ret = _multi_svd_norm(x, row_axis, col_axis, sum) 

2583 else: 

2584 raise ValueError("Invalid norm order for matrices.") 

2585 if keepdims: 

2586 ret_shape = list(x.shape) 

2587 ret_shape[axis[0]] = 1 

2588 ret_shape[axis[1]] = 1 

2589 ret = ret.reshape(ret_shape) 

2590 return ret 

2591 else: 

2592 raise ValueError("Improper number of dimensions to norm.") 

2593 

2594 

2595# multi_dot 

2596 

2597def _multidot_dispatcher(arrays, *, out=None): 

2598 yield from arrays 

2599 yield out 

2600 

2601 

2602@array_function_dispatch(_multidot_dispatcher) 

2603def multi_dot(arrays, *, out=None): 

2604 """ 

2605 Compute the dot product of two or more arrays in a single function call, 

2606 while automatically selecting the fastest evaluation order. 

2607 

2608 `multi_dot` chains `numpy.dot` and uses optimal parenthesization 

2609 of the matrices [1]_ [2]_. Depending on the shapes of the matrices, 

2610 this can speed up the multiplication a lot. 

2611 

2612 If the first argument is 1-D it is treated as a row vector. 

2613 If the last argument is 1-D it is treated as a column vector. 

2614 The other arguments must be 2-D. 

2615 

2616 Think of `multi_dot` as:: 

2617 

2618 def multi_dot(arrays): return functools.reduce(np.dot, arrays) 

2619 

2620 

2621 Parameters 

2622 ---------- 

2623 arrays : sequence of array_like 

2624 If the first argument is 1-D it is treated as row vector. 

2625 If the last argument is 1-D it is treated as column vector. 

2626 The other arguments must be 2-D. 

2627 out : ndarray, optional 

2628 Output argument. This must have the exact kind that would be returned 

2629 if it was not used. In particular, it must have the right type, must be 

2630 C-contiguous, and its dtype must be the dtype that would be returned 

2631 for `dot(a, b)`. This is a performance feature. Therefore, if these 

2632 conditions are not met, an exception is raised, instead of attempting 

2633 to be flexible. 

2634 

2635 .. versionadded:: 1.19.0 

2636 

2637 Returns 

2638 ------- 

2639 output : ndarray 

2640 Returns the dot product of the supplied arrays. 

2641 

2642 See Also 

2643 -------- 

2644 numpy.dot : dot multiplication with two arguments. 

2645 

2646 References 

2647 ---------- 

2648 

2649 .. [1] Cormen, "Introduction to Algorithms", Chapter 15.2, p. 370-378 

2650 .. [2] https://en.wikipedia.org/wiki/Matrix_chain_multiplication 

2651 

2652 Examples 

2653 -------- 

2654 `multi_dot` allows you to write:: 

2655 

2656 >>> from numpy.linalg import multi_dot 

2657 >>> # Prepare some data 

2658 >>> A = np.random.random((10000, 100)) 

2659 >>> B = np.random.random((100, 1000)) 

2660 >>> C = np.random.random((1000, 5)) 

2661 >>> D = np.random.random((5, 333)) 

2662 >>> # the actual dot multiplication 

2663 >>> _ = multi_dot([A, B, C, D]) 

2664 

2665 instead of:: 

2666 

2667 >>> _ = np.dot(np.dot(np.dot(A, B), C), D) 

2668 >>> # or 

2669 >>> _ = A.dot(B).dot(C).dot(D) 

2670 

2671 Notes 

2672 ----- 

2673 The cost for a matrix multiplication can be calculated with the 

2674 following function:: 

2675 

2676 def cost(A, B): 

2677 return A.shape[0] * A.shape[1] * B.shape[1] 

2678 

2679 Assume we have three matrices 

2680 :math:`A_{10x100}, B_{100x5}, C_{5x50}`. 

2681 

2682 The costs for the two different parenthesizations are as follows:: 

2683 

2684 cost((AB)C) = 10*100*5 + 10*5*50 = 5000 + 2500 = 7500 

2685 cost(A(BC)) = 10*100*50 + 100*5*50 = 50000 + 25000 = 75000 

2686 

2687 """ 

2688 n = len(arrays) 

2689 # optimization only makes sense for len(arrays) > 2 

2690 if n < 2: 

2691 raise ValueError("Expecting at least two arrays.") 

2692 elif n == 2: 

2693 return dot(arrays[0], arrays[1], out=out) 

2694 

2695 arrays = [asanyarray(a) for a in arrays] 

2696 

2697 # save original ndim to reshape the result array into the proper form later 

2698 ndim_first, ndim_last = arrays[0].ndim, arrays[-1].ndim 

2699 # Explicitly convert vectors to 2D arrays to keep the logic of the internal 

2700 # _multi_dot_* functions as simple as possible. 

2701 if arrays[0].ndim == 1: 

2702 arrays[0] = atleast_2d(arrays[0]) 

2703 if arrays[-1].ndim == 1: 

2704 arrays[-1] = atleast_2d(arrays[-1]).T 

2705 _assert_2d(*arrays) 

2706 

2707 # _multi_dot_three is much faster than _multi_dot_matrix_chain_order 

2708 if n == 3: 

2709 result = _multi_dot_three(arrays[0], arrays[1], arrays[2], out=out) 

2710 else: 

2711 order = _multi_dot_matrix_chain_order(arrays) 

2712 result = _multi_dot(arrays, order, 0, n - 1, out=out) 

2713 

2714 # return proper shape 

2715 if ndim_first == 1 and ndim_last == 1: 

2716 return result[0, 0] # scalar 

2717 elif ndim_first == 1 or ndim_last == 1: 

2718 return result.ravel() # 1-D 

2719 else: 

2720 return result 

2721 

2722 

2723def _multi_dot_three(A, B, C, out=None): 

2724 """ 

2725 Find the best order for three arrays and do the multiplication. 

2726 

2727 For three arguments `_multi_dot_three` is approximately 15 times faster 

2728 than `_multi_dot_matrix_chain_order` 

2729 

2730 """ 

2731 a0, a1b0 = A.shape 

2732 b1c0, c1 = C.shape 

2733 # cost1 = cost((AB)C) = a0*a1b0*b1c0 + a0*b1c0*c1 

2734 cost1 = a0 * b1c0 * (a1b0 + c1) 

2735 # cost2 = cost(A(BC)) = a1b0*b1c0*c1 + a0*a1b0*c1 

2736 cost2 = a1b0 * c1 * (a0 + b1c0) 

2737 

2738 if cost1 < cost2: 

2739 return dot(dot(A, B), C, out=out) 

2740 else: 

2741 return dot(A, dot(B, C), out=out) 

2742 

2743 

2744def _multi_dot_matrix_chain_order(arrays, return_costs=False): 

2745 """ 

2746 Return a np.array that encodes the optimal order of mutiplications. 

2747 

2748 The optimal order array is then used by `_multi_dot()` to do the 

2749 multiplication. 

2750 

2751 Also return the cost matrix if `return_costs` is `True` 

2752 

2753 The implementation CLOSELY follows Cormen, "Introduction to Algorithms", 

2754 Chapter 15.2, p. 370-378. Note that Cormen uses 1-based indices. 

2755 

2756 cost[i, j] = min([ 

2757 cost[prefix] + cost[suffix] + cost_mult(prefix, suffix) 

2758 for k in range(i, j)]) 

2759 

2760 """ 

2761 n = len(arrays) 

2762 # p stores the dimensions of the matrices 

2763 # Example for p: A_{10x100}, B_{100x5}, C_{5x50} --> p = [10, 100, 5, 50] 

2764 p = [a.shape[0] for a in arrays] + [arrays[-1].shape[1]] 

2765 # m is a matrix of costs of the subproblems 

2766 # m[i,j]: min number of scalar multiplications needed to compute A_{i..j} 

2767 m = zeros((n, n), dtype=double) 

2768 # s is the actual ordering 

2769 # s[i, j] is the value of k at which we split the product A_i..A_j 

2770 s = empty((n, n), dtype=intp) 

2771 

2772 for l in range(1, n): 

2773 for i in range(n - l): 

2774 j = i + l 

2775 m[i, j] = Inf 

2776 for k in range(i, j): 

2777 q = m[i, k] + m[k+1, j] + p[i]*p[k+1]*p[j+1] 

2778 if q < m[i, j]: 

2779 m[i, j] = q 

2780 s[i, j] = k # Note that Cormen uses 1-based index 

2781 

2782 return (s, m) if return_costs else s 

2783 

2784 

2785def _multi_dot(arrays, order, i, j, out=None): 

2786 """Actually do the multiplication with the given order.""" 

2787 if i == j: 

2788 # the initial call with non-None out should never get here 

2789 assert out is None 

2790 

2791 return arrays[i] 

2792 else: 

2793 return dot(_multi_dot(arrays, order, i, order[i, j]), 

2794 _multi_dot(arrays, order, order[i, j] + 1, j), 

2795 out=out)