Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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, fastCopyAndTranspose, 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) 

29from numpy.core.multiarray import normalize_axis_index 

30from numpy.core.overrides import set_module 

31from numpy.core import overrides 

32from numpy.lib.twodim_base import triu, eye 

33from numpy.linalg import lapack_lite, _umath_linalg 

34 

35 

36array_function_dispatch = functools.partial( 

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

38 

39 

40fortran_int = intc 

41 

42 

43@set_module('numpy.linalg') 

44class LinAlgError(Exception): 

45 """ 

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

47 

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

49 class, programmatically raised in linalg functions when a Linear 

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

51 function. 

52 

53 Parameters 

54 ---------- 

55 None 

56 

57 Examples 

58 -------- 

59 >>> from numpy import linalg as LA 

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

61 Traceback (most recent call last): 

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

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

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

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

66 in solve 

67 raise LinAlgError('Singular matrix') 

68 numpy.linalg.LinAlgError: Singular matrix 

69 

70 """ 

71 

72 

73def _determine_error_states(): 

74 errobj = geterrobj() 

75 bufsize = errobj[0] 

76 

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

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

79 invalid_call_errmask = geterrobj()[1] 

80 

81 return [bufsize, invalid_call_errmask, None] 

82 

83# Dealing with errors in _umath_linalg 

84_linalg_error_extobj = _determine_error_states() 

85del _determine_error_states 

86 

87def _raise_linalgerror_singular(err, flag): 

88 raise LinAlgError("Singular matrix") 

89 

90def _raise_linalgerror_nonposdef(err, flag): 

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

92 

93def _raise_linalgerror_eigenvalues_nonconvergence(err, flag): 

94 raise LinAlgError("Eigenvalues did not converge") 

95 

96def _raise_linalgerror_svd_nonconvergence(err, flag): 

97 raise LinAlgError("SVD did not converge") 

98 

99def _raise_linalgerror_lstsq(err, flag): 

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

101 

102def get_linalg_error_extobj(callback): 

103 extobj = list(_linalg_error_extobj) # make a copy 

104 extobj[2] = callback 

105 return extobj 

106 

107def _makearray(a): 

108 new = asarray(a) 

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

110 return new, wrap 

111 

112def isComplexType(t): 

113 return issubclass(t, complexfloating) 

114 

115_real_types_map = {single : single, 

116 double : double, 

117 csingle : single, 

118 cdouble : double} 

119 

120_complex_types_map = {single : csingle, 

121 double : cdouble, 

122 csingle : csingle, 

123 cdouble : cdouble} 

124 

125def _realType(t, default=double): 

126 return _real_types_map.get(t, default) 

127 

128def _complexType(t, default=cdouble): 

129 return _complex_types_map.get(t, default) 

130 

131def _linalgRealType(t): 

132 """Cast the type t to either double or cdouble.""" 

133 return double 

134 

135def _commonType(*arrays): 

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

137 result_type = single 

138 is_complex = False 

139 for a in arrays: 

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

141 if isComplexType(a.dtype.type): 

142 is_complex = True 

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

144 if rt is None: 

145 # unsupported inexact scalar 

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

147 (a.dtype.name,)) 

148 else: 

149 rt = double 

150 if rt is double: 

151 result_type = double 

152 if is_complex: 

153 t = cdouble 

154 result_type = _complex_types_map[result_type] 

155 else: 

156 t = double 

157 return t, result_type 

158 

159 

160# _fastCopyAndTranpose assumes the input is 2D (as all the calls in here are). 

161 

162_fastCT = fastCopyAndTranspose 

163 

164def _to_native_byte_order(*arrays): 

165 ret = [] 

166 for arr in arrays: 

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

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

169 else: 

170 ret.append(arr) 

171 if len(ret) == 1: 

172 return ret[0] 

173 else: 

174 return ret 

175 

176def _fastCopyAndTranspose(type, *arrays): 

177 cast_arrays = () 

178 for a in arrays: 

179 if a.dtype.type is type: 

180 cast_arrays = cast_arrays + (_fastCT(a),) 

181 else: 

182 cast_arrays = cast_arrays + (_fastCT(a.astype(type)),) 

183 if len(cast_arrays) == 1: 

184 return cast_arrays[0] 

185 else: 

186 return cast_arrays 

187 

188def _assert_2d(*arrays): 

189 for a in arrays: 

190 if a.ndim != 2: 

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

192 'two-dimensional' % a.ndim) 

193 

194def _assert_stacked_2d(*arrays): 

195 for a in arrays: 

196 if a.ndim < 2: 

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

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

199 

200def _assert_stacked_square(*arrays): 

201 for a in arrays: 

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

203 if m != n: 

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

205 

206def _assert_finite(*arrays): 

207 for a in arrays: 

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

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

210 

211def _is_empty_2d(arr): 

212 # check size first for efficiency 

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

214 

215 

216def transpose(a): 

217 """ 

218 Transpose each matrix in a stack of matrices. 

219 

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

221 them 

222 

223 Parameters 

224 ---------- 

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

226 

227 Returns 

228 ------- 

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

230 """ 

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

232 

233# Linear equations 

234 

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

236 return (a, b) 

237 

238 

239@array_function_dispatch(_tensorsolve_dispatcher) 

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

241 """ 

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

243 

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

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

246 ``tensordot(a, x, axes=b.ndim)``. 

247 

248 Parameters 

249 ---------- 

250 a : array_like 

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

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

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

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

255 'square'). 

256 b : array_like 

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

258 axes : tuple of ints, optional 

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

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

261 

262 Returns 

263 ------- 

264 x : ndarray, shape Q 

265 

266 Raises 

267 ------ 

268 LinAlgError 

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

270 

271 See Also 

272 -------- 

273 numpy.tensordot, tensorinv, numpy.einsum 

274 

275 Examples 

276 -------- 

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

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

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

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

281 >>> x.shape 

282 (2, 3, 4) 

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

284 True 

285 

286 """ 

287 a, wrap = _makearray(a) 

288 b = asarray(b) 

289 an = a.ndim 

290 

291 if axes is not None: 

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

293 for k in axes: 

294 allaxes.remove(k) 

295 allaxes.insert(an, k) 

296 a = a.transpose(allaxes) 

297 

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

299 prod = 1 

300 for k in oldshape: 

301 prod *= k 

302 

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

304 b = b.ravel() 

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

306 res.shape = oldshape 

307 return res 

308 

309 

310def _solve_dispatcher(a, b): 

311 return (a, b) 

312 

313 

314@array_function_dispatch(_solve_dispatcher) 

315def solve(a, b): 

316 """ 

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

318 

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

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

321 

322 Parameters 

323 ---------- 

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

325 Coefficient matrix. 

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

327 Ordinate or "dependent variable" values. 

328 

329 Returns 

330 ------- 

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

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

333 

334 Raises 

335 ------ 

336 LinAlgError 

337 If `a` is singular or not square. 

338 

339 See Also 

340 -------- 

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

342 

343 Notes 

344 ----- 

345 

346 .. versionadded:: 1.8.0 

347 

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

349 details. 

350 

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

352 

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

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

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

356 system/equation. 

357 

358 References 

359 ---------- 

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

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

362 

363 Examples 

364 -------- 

365 Solve the system of equations ``3 * x0 + x1 = 9`` and ``x0 + 2 * x1 = 8``: 

366 

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

368 >>> b = np.array([9,8]) 

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

370 >>> x 

371 array([2., 3.]) 

372 

373 Check that the solution is correct: 

374 

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

376 True 

377 

378 """ 

379 a, _ = _makearray(a) 

380 _assert_stacked_2d(a) 

381 _assert_stacked_square(a) 

382 b, wrap = _makearray(b) 

383 t, result_t = _commonType(a, b) 

384 

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

386 # match exactly 

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

388 gufunc = _umath_linalg.solve1 

389 else: 

390 gufunc = _umath_linalg.solve 

391 

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

393 extobj = get_linalg_error_extobj(_raise_linalgerror_singular) 

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

395 

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

397 

398 

399def _tensorinv_dispatcher(a, ind=None): 

400 return (a,) 

401 

402 

403@array_function_dispatch(_tensorinv_dispatcher) 

404def tensorinv(a, ind=2): 

405 """ 

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

407 

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

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

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

411 tensordot operation. 

412 

413 Parameters 

414 ---------- 

415 a : array_like 

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

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

418 ind : int, optional 

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

420 Must be a positive integer, default is 2. 

421 

422 Returns 

423 ------- 

424 b : ndarray 

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

426 

427 Raises 

428 ------ 

429 LinAlgError 

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

431 

432 See Also 

433 -------- 

434 numpy.tensordot, tensorsolve 

435 

436 Examples 

437 -------- 

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

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

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

441 >>> ainv.shape 

442 (8, 3, 4, 6) 

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

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

445 True 

446 

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

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

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

450 >>> ainv.shape 

451 (8, 3, 24) 

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

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

454 True 

455 

456 """ 

457 a = asarray(a) 

458 oldshape = a.shape 

459 prod = 1 

460 if ind > 0: 

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

462 for k in oldshape[ind:]: 

463 prod *= k 

464 else: 

465 raise ValueError("Invalid ind argument.") 

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

467 ia = inv(a) 

468 return ia.reshape(*invshape) 

469 

470 

471# Matrix inversion 

472 

473def _unary_dispatcher(a): 

474 return (a,) 

475 

476 

477@array_function_dispatch(_unary_dispatcher) 

478def inv(a): 

479 """ 

480 Compute the (multiplicative) inverse of a matrix. 

481 

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

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

484 

485 Parameters 

486 ---------- 

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

488 Matrix to be inverted. 

489 

490 Returns 

491 ------- 

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

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

494 

495 Raises 

496 ------ 

497 LinAlgError 

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

499 

500 See Also 

501 -------- 

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

503 

504 Notes 

505 ----- 

506 

507 .. versionadded:: 1.8.0 

508 

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

510 details. 

511 

512 Examples 

513 -------- 

514 >>> from numpy.linalg import inv 

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

516 >>> ainv = inv(a) 

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

518 True 

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

520 True 

521 

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

523 

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

525 >>> ainv 

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

527 [ 1.5, -0.5]]) 

528 

529 Inverses of several matrices can be computed at once: 

530 

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

532 >>> inv(a) 

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

534 [ 1.5 , -0.5 ]], 

535 [[-1.25, 0.75], 

536 [ 0.75, -0.25]]]) 

537 

538 """ 

539 a, wrap = _makearray(a) 

540 _assert_stacked_2d(a) 

541 _assert_stacked_square(a) 

542 t, result_t = _commonType(a) 

543 

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

545 extobj = get_linalg_error_extobj(_raise_linalgerror_singular) 

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

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

548 

549 

550def _matrix_power_dispatcher(a, n): 

551 return (a,) 

552 

553 

554@array_function_dispatch(_matrix_power_dispatcher) 

555def matrix_power(a, n): 

556 """ 

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

558 

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

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

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

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

563 

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

565 

566 Parameters 

567 ---------- 

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

569 Matrix to be "powered". 

570 n : int 

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

572 negative, or zero. 

573 

574 Returns 

575 ------- 

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

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

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

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

580 negative the elements are floating-point. 

581 

582 Raises 

583 ------ 

584 LinAlgError 

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

586 be inverted numerically. 

587 

588 Examples 

589 -------- 

590 >>> from numpy.linalg import matrix_power 

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

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

593 array([[ 0, -1], 

594 [ 1, 0]]) 

595 >>> matrix_power(i, 0) 

596 array([[1, 0], 

597 [0, 1]]) 

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

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

600 [-1., 0.]]) 

601 

602 Somewhat more sophisticated example 

603 

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

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

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

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

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

609 [ 1., 0., 0., 0.], 

610 [ 0., 0., 0., 1.], 

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

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

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

614 [ 0., -1., 0., 0.], 

615 [ 0., 0., -1., 0.], 

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

617 

618 """ 

619 a = asanyarray(a) 

620 _assert_stacked_2d(a) 

621 _assert_stacked_square(a) 

622 

623 try: 

624 n = operator.index(n) 

625 except TypeError as e: 

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

627 

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

629 # the current implementation of matmul using einsum 

630 if a.dtype != object: 

631 fmatmul = matmul 

632 elif a.ndim == 2: 

633 fmatmul = dot 

634 else: 

635 raise NotImplementedError( 

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

637 

638 if n == 0: 

639 a = empty_like(a) 

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

641 return a 

642 

643 elif n < 0: 

644 a = inv(a) 

645 n = abs(n) 

646 

647 # short-cuts. 

648 if n == 1: 

649 return a 

650 

651 elif n == 2: 

652 return fmatmul(a, a) 

653 

654 elif n == 3: 

655 return fmatmul(fmatmul(a, a), a) 

656 

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

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

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

660 z = result = None 

661 while n > 0: 

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

663 n, bit = divmod(n, 2) 

664 if bit: 

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

666 

667 return result 

668 

669 

670# Cholesky decomposition 

671 

672 

673@array_function_dispatch(_unary_dispatcher) 

674def cholesky(a): 

675 """ 

676 Cholesky decomposition. 

677 

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

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

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

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

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

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

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

685 

686 Parameters 

687 ---------- 

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

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

690 input matrix. 

691 

692 Returns 

693 ------- 

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

695 Upper or lower-triangular Cholesky factor of `a`. Returns a 

696 matrix object if `a` is a matrix object. 

697 

698 Raises 

699 ------ 

700 LinAlgError 

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

702 positive-definite. 

703 

704 See Also 

705 -------- 

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

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

708 positive-definite matrix. 

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

710 `scipy.linalg.cho_solve`. 

711 

712 Notes 

713 ----- 

714 

715 .. versionadded:: 1.8.0 

716 

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

718 details. 

719 

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

721 

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

723 

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

725 

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

727 

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

729 

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

731 

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

733 

734 Examples 

735 -------- 

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

737 >>> A 

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

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

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

741 >>> L 

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

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

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

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

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

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

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

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

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

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

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

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

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

755 

756 """ 

757 extobj = get_linalg_error_extobj(_raise_linalgerror_nonposdef) 

758 gufunc = _umath_linalg.cholesky_lo 

759 a, wrap = _makearray(a) 

760 _assert_stacked_2d(a) 

761 _assert_stacked_square(a) 

762 t, result_t = _commonType(a) 

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

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

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

766 

767 

768# QR decomposition 

769 

770def _qr_dispatcher(a, mode=None): 

771 return (a,) 

772 

773 

774@array_function_dispatch(_qr_dispatcher) 

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

776 """ 

777 Compute the qr factorization of a matrix. 

778 

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

780 upper-triangular. 

781 

782 Parameters 

783 ---------- 

784 a : array_like, shape (M, N) 

785 Matrix to be factored. 

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

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

788 

789 * 'reduced' : returns q, r with dimensions (M, K), (K, N) (default) 

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

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

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

793 

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

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

796 maintain backward compatibility with earlier versions of numpy both 

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

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

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

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

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

802 explanation. 

803 

804 

805 Returns 

806 ------- 

807 q : ndarray of float or complex, optional 

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

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

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

811 case. 

812 r : ndarray of float or complex, optional 

813 The upper-triangular matrix. 

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

815 The array h contains the Householder reflectors that generate q 

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

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

818 

819 Raises 

820 ------ 

821 LinAlgError 

822 If factoring fails. 

823 

824 See Also 

825 -------- 

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

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

828 

829 Notes 

830 ----- 

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

832 ``dorgqr``, and ``zungqr``. 

833 

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

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

836 

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

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

839 

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

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

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

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

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

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

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

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

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

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

850 in lapack_lite and just await the necessary work. 

851 

852 Examples 

853 -------- 

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

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

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

857 True 

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

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

860 True 

861 

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

863 problems 

864 

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

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

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

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

869 

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

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

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

873 

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

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

876 however, we simply use `lstsq`.) 

877 

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

879 >>> A 

880 array([[0, 1], 

881 [1, 1], 

882 [1, 1], 

883 [2, 1]]) 

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

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

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

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

888 array([ 1.1e-16, 1.0e+00]) 

889 

890 """ 

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

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

893 # 2013-04-01, 1.8 

894 msg = "".join(( 

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

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

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

898 mode = 'reduced' 

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

900 # 2013-04-01, 1.8 

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

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

903 mode = 'economic' 

904 else: 

905 raise ValueError("Unrecognized mode '%s'" % mode) 

906 

907 a, wrap = _makearray(a) 

908 _assert_2d(a) 

909 m, n = a.shape 

910 t, result_t = _commonType(a) 

911 a = _fastCopyAndTranspose(t, a) 

912 a = _to_native_byte_order(a) 

913 mn = min(m, n) 

914 tau = zeros((mn,), t) 

915 

916 if isComplexType(t): 

917 lapack_routine = lapack_lite.zgeqrf 

918 routine_name = 'zgeqrf' 

919 else: 

920 lapack_routine = lapack_lite.dgeqrf 

921 routine_name = 'dgeqrf' 

922 

923 # calculate optimal size of work data 'work' 

924 lwork = 1 

925 work = zeros((lwork,), t) 

926 results = lapack_routine(m, n, a, max(1, m), tau, work, -1, 0) 

927 if results['info'] != 0: 

928 raise LinAlgError('%s returns %d' % (routine_name, results['info'])) 

929 

930 # do qr decomposition 

931 lwork = max(1, n, int(abs(work[0]))) 

932 work = zeros((lwork,), t) 

933 results = lapack_routine(m, n, a, max(1, m), tau, work, lwork, 0) 

934 if results['info'] != 0: 

935 raise LinAlgError('%s returns %d' % (routine_name, results['info'])) 

936 

937 # handle modes that don't return q 

938 if mode == 'r': 

939 r = _fastCopyAndTranspose(result_t, a[:, :mn]) 

940 return wrap(triu(r)) 

941 

942 if mode == 'raw': 

943 return a, tau 

944 

945 if mode == 'economic': 

946 if t != result_t : 

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

948 return wrap(a.T) 

949 

950 # generate q from a 

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

952 mc = m 

953 q = empty((m, m), t) 

954 else: 

955 mc = mn 

956 q = empty((n, m), t) 

957 q[:n] = a 

958 

959 if isComplexType(t): 

960 lapack_routine = lapack_lite.zungqr 

961 routine_name = 'zungqr' 

962 else: 

963 lapack_routine = lapack_lite.dorgqr 

964 routine_name = 'dorgqr' 

965 

966 # determine optimal lwork 

967 lwork = 1 

968 work = zeros((lwork,), t) 

969 results = lapack_routine(m, mc, mn, q, max(1, m), tau, work, -1, 0) 

970 if results['info'] != 0: 

971 raise LinAlgError('%s returns %d' % (routine_name, results['info'])) 

972 

973 # compute q 

974 lwork = max(1, n, int(abs(work[0]))) 

975 work = zeros((lwork,), t) 

976 results = lapack_routine(m, mc, mn, q, max(1, m), tau, work, lwork, 0) 

977 if results['info'] != 0: 

978 raise LinAlgError('%s returns %d' % (routine_name, results['info'])) 

979 

980 q = _fastCopyAndTranspose(result_t, q[:mc]) 

981 r = _fastCopyAndTranspose(result_t, a[:, :mc]) 

982 

983 return wrap(q), wrap(triu(r)) 

984 

985 

986# Eigenvalues 

987 

988 

989@array_function_dispatch(_unary_dispatcher) 

990def eigvals(a): 

991 """ 

992 Compute the eigenvalues of a general matrix. 

993 

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

995 returned. 

996 

997 Parameters 

998 ---------- 

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

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

1001 

1002 Returns 

1003 ------- 

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

1005 The eigenvalues, each repeated according to its multiplicity. 

1006 They are not necessarily ordered, nor are they necessarily 

1007 real for real matrices. 

1008 

1009 Raises 

1010 ------ 

1011 LinAlgError 

1012 If the eigenvalue computation does not converge. 

1013 

1014 See Also 

1015 -------- 

1016 eig : eigenvalues and right eigenvectors of general arrays 

1017 eigvalsh : eigenvalues of real symmetric or complex Hermitian 

1018 (conjugate symmetric) arrays. 

1019 eigh : eigenvalues and eigenvectors of real symmetric or complex 

1020 Hermitian (conjugate symmetric) arrays. 

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

1022 

1023 Notes 

1024 ----- 

1025 

1026 .. versionadded:: 1.8.0 

1027 

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

1029 details. 

1030 

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

1032 the eigenvalues and eigenvectors of general square arrays. 

1033 

1034 Examples 

1035 -------- 

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

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

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

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

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

1041 ``A``: 

1042 

1043 >>> from numpy import linalg as LA 

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

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

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

1047 (1.0, 1.0, 0.0) 

1048 

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

1050 

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

1052 >>> LA.eigvals(D) 

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

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

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

1056 >>> LA.eigvals(A) 

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

1058 

1059 """ 

1060 a, wrap = _makearray(a) 

1061 _assert_stacked_2d(a) 

1062 _assert_stacked_square(a) 

1063 _assert_finite(a) 

1064 t, result_t = _commonType(a) 

1065 

1066 extobj = get_linalg_error_extobj( 

1067 _raise_linalgerror_eigenvalues_nonconvergence) 

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

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

1070 

1071 if not isComplexType(t): 

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

1073 w = w.real 

1074 result_t = _realType(result_t) 

1075 else: 

1076 result_t = _complexType(result_t) 

1077 

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

1079 

1080 

1081def _eigvalsh_dispatcher(a, UPLO=None): 

1082 return (a,) 

1083 

1084 

1085@array_function_dispatch(_eigvalsh_dispatcher) 

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

1087 """ 

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

1089 

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

1091 

1092 Parameters 

1093 ---------- 

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

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

1096 computed. 

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

1098 Specifies whether the calculation is done with the lower triangular 

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

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

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

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

1103 will always be treated as zero. 

1104 

1105 Returns 

1106 ------- 

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

1108 The eigenvalues in ascending order, each repeated according to 

1109 its multiplicity. 

1110 

1111 Raises 

1112 ------ 

1113 LinAlgError 

1114 If the eigenvalue computation does not converge. 

1115 

1116 See Also 

1117 -------- 

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

1119 (conjugate symmetric) arrays. 

1120 eigvals : eigenvalues of general real or complex arrays. 

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

1122 arrays. 

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

1124 

1125 Notes 

1126 ----- 

1127 

1128 .. versionadded:: 1.8.0 

1129 

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

1131 details. 

1132 

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

1134 

1135 Examples 

1136 -------- 

1137 >>> from numpy import linalg as LA 

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

1139 >>> LA.eigvalsh(a) 

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

1141 

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

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

1144 >>> a 

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

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

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

1148 >>> # with: 

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

1150 >>> b 

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

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

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

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

1155 >>> wa; wb 

1156 array([1., 6.]) 

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

1158 

1159 """ 

1160 UPLO = UPLO.upper() 

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

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

1163 

1164 extobj = get_linalg_error_extobj( 

1165 _raise_linalgerror_eigenvalues_nonconvergence) 

1166 if UPLO == 'L': 

1167 gufunc = _umath_linalg.eigvalsh_lo 

1168 else: 

1169 gufunc = _umath_linalg.eigvalsh_up 

1170 

1171 a, wrap = _makearray(a) 

1172 _assert_stacked_2d(a) 

1173 _assert_stacked_square(a) 

1174 t, result_t = _commonType(a) 

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

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

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

1178 

1179def _convertarray(a): 

1180 t, result_t = _commonType(a) 

1181 a = _fastCT(a.astype(t)) 

1182 return a, t, result_t 

1183 

1184 

1185# Eigenvectors 

1186 

1187 

1188@array_function_dispatch(_unary_dispatcher) 

1189def eig(a): 

1190 """ 

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

1192 

1193 Parameters 

1194 ---------- 

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

1196 Matrices for which the eigenvalues and right eigenvectors will 

1197 be computed 

1198 

1199 Returns 

1200 ------- 

1201 w : (..., M) array 

1202 The eigenvalues, each repeated according to its multiplicity. 

1203 The eigenvalues are not necessarily ordered. The resulting 

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

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

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

1207 part) or occur in conjugate pairs 

1208 

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

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

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

1212 eigenvalue ``w[i]``. 

1213 

1214 Raises 

1215 ------ 

1216 LinAlgError 

1217 If the eigenvalue computation does not converge. 

1218 

1219 See Also 

1220 -------- 

1221 eigvals : eigenvalues of a non-symmetric array. 

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

1223 Hermitian (conjugate symmetric) array. 

1224 eigvalsh : eigenvalues of a real symmetric or complex Hermitian 

1225 (conjugate symmetric) array. 

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

1227 generalized eigenvalue problem. 

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

1229 normal matrices. 

1230 

1231 Notes 

1232 ----- 

1233 

1234 .. versionadded:: 1.8.0 

1235 

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

1237 details. 

1238 

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

1240 the eigenvalues and eigenvectors of general square arrays. 

1241 

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

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

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

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

1246 

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

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

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

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

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

1252 

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

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

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

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

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

1258 rest is roundoff error. 

1259 

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

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

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

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

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

1265 of each other. 

1266 

1267 References 

1268 ---------- 

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

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

1271 

1272 Examples 

1273 -------- 

1274 >>> from numpy import linalg as LA 

1275 

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

1277 

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

1279 >>> w; v 

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

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

1282 [0., 1., 0.], 

1283 [0., 0., 1.]]) 

1284 

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

1286 e-values are complex conjugates of each other. 

1287 

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

1289 >>> w; v 

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

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

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

1293 

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

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

1296 

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

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

1299 >>> w; v 

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

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

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

1303 

1304 Be careful about round-off error! 

1305 

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

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

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

1309 >>> w; v 

1310 array([1., 1.]) 

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

1312 [0., 1.]]) 

1313 

1314 """ 

1315 a, wrap = _makearray(a) 

1316 _assert_stacked_2d(a) 

1317 _assert_stacked_square(a) 

1318 _assert_finite(a) 

1319 t, result_t = _commonType(a) 

1320 

1321 extobj = get_linalg_error_extobj( 

1322 _raise_linalgerror_eigenvalues_nonconvergence) 

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

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

1325 

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

1327 w = w.real 

1328 vt = vt.real 

1329 result_t = _realType(result_t) 

1330 else: 

1331 result_t = _complexType(result_t) 

1332 

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

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

1335 

1336 

1337@array_function_dispatch(_eigvalsh_dispatcher) 

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

1339 """ 

1340 Return the eigenvalues and eigenvectors of a complex Hermitian 

1341 (conjugate symmetric) or a real symmetric matrix. 

1342 

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

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

1345 corresponding eigenvectors (in columns). 

1346 

1347 Parameters 

1348 ---------- 

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

1350 Hermitian or real symmetric matrices whose eigenvalues and 

1351 eigenvectors are to be computed. 

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

1353 Specifies whether the calculation is done with the lower triangular 

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

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

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

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

1358 will always be treated as zero. 

1359 

1360 Returns 

1361 ------- 

1362 w : (..., M) ndarray 

1363 The eigenvalues in ascending order, each repeated according to 

1364 its multiplicity. 

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

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

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

1368 a matrix object. 

1369 

1370 Raises 

1371 ------ 

1372 LinAlgError 

1373 If the eigenvalue computation does not converge. 

1374 

1375 See Also 

1376 -------- 

1377 eigvalsh : eigenvalues of real symmetric or complex Hermitian 

1378 (conjugate symmetric) arrays. 

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

1380 eigvals : eigenvalues of non-symmetric arrays. 

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

1382 generalized eigenvalue problem). 

1383 

1384 Notes 

1385 ----- 

1386 

1387 .. versionadded:: 1.8.0 

1388 

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

1390 details. 

1391 

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

1393 ``_heevd``. 

1394 

1395 The eigenvalues of real symmetric or complex Hermitian matrices are 

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

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

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

1399 

1400 References 

1401 ---------- 

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

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

1404 

1405 Examples 

1406 -------- 

1407 >>> from numpy import linalg as LA 

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

1409 >>> a 

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

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

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

1413 >>> w; v 

1414 array([0.17157288, 5.82842712]) 

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

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

1417 

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

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

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

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

1422 

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

1424 >>> A 

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

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

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

1428 >>> w; v 

1429 array([0.17157288, 5.82842712]) 

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

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

1432 

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

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

1435 >>> a 

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

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

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

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

1440 >>> b 

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

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

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

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

1445 >>> wa; wb 

1446 array([1., 6.]) 

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

1448 >>> va; vb 

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

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

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

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

1453 """ 

1454 UPLO = UPLO.upper() 

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

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

1457 

1458 a, wrap = _makearray(a) 

1459 _assert_stacked_2d(a) 

1460 _assert_stacked_square(a) 

1461 t, result_t = _commonType(a) 

1462 

1463 extobj = get_linalg_error_extobj( 

1464 _raise_linalgerror_eigenvalues_nonconvergence) 

1465 if UPLO == 'L': 

1466 gufunc = _umath_linalg.eigh_lo 

1467 else: 

1468 gufunc = _umath_linalg.eigh_up 

1469 

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

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

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

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

1474 return w, wrap(vt) 

1475 

1476 

1477# Singular value decomposition 

1478 

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

1480 return (a,) 

1481 

1482 

1483@array_function_dispatch(_svd_dispatcher) 

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

1485 """ 

1486 Singular Value Decomposition. 

1487 

1488 When `a` is a 2D array, it is factorized as ``u @ np.diag(s) @ vh 

1489 = (u * s) @ vh``, where `u` and `vh` are 2D unitary arrays and `s` is a 1D 

1490 array of `a`'s singular values. When `a` is higher-dimensional, SVD is 

1491 applied in stacked mode as explained below. 

1492 

1493 Parameters 

1494 ---------- 

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

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

1497 full_matrices : bool, optional 

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

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

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

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

1502 compute_uv : bool, optional 

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

1504 by default. 

1505 hermitian : bool, optional 

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

1507 enabling a more efficient method for finding singular values. 

1508 Defaults to False. 

1509 

1510 .. versionadded:: 1.17.0 

1511 

1512 Returns 

1513 ------- 

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

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

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

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

1518 `compute_uv` is True. 

1519 s : (..., K) array 

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

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

1522 size as those of the input `a`. 

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

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

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

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

1527 `compute_uv` is True. 

1528 

1529 Raises 

1530 ------ 

1531 LinAlgError 

1532 If SVD computation does not converge. 

1533 

1534 See Also 

1535 -------- 

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

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

1538 

1539 Notes 

1540 ----- 

1541 

1542 .. versionchanged:: 1.8.0 

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

1544 details. 

1545 

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

1547 

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

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

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

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

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

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

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

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

1556 

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

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

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

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

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

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

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

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

1565 

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

1567 all the return values. 

1568 

1569 Examples 

1570 -------- 

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

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

1573 

1574 Reconstruction based on full SVD, 2D case: 

1575 

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

1577 >>> u.shape, s.shape, vh.shape 

1578 ((9, 9), (6,), (6, 6)) 

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

1580 True 

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

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

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

1584 True 

1585 

1586 Reconstruction based on reduced SVD, 2D case: 

1587 

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

1589 >>> u.shape, s.shape, vh.shape 

1590 ((9, 6), (6,), (6, 6)) 

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

1592 True 

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

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

1595 True 

1596 

1597 Reconstruction based on full SVD, 4D case: 

1598 

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

1600 >>> u.shape, s.shape, vh.shape 

1601 ((2, 7, 8, 8), (2, 7, 3), (2, 7, 3, 3)) 

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

1603 True 

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

1605 True 

1606 

1607 Reconstruction based on reduced SVD, 4D case: 

1608 

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

1610 >>> u.shape, s.shape, vh.shape 

1611 ((2, 7, 8, 3), (2, 7, 3), (2, 7, 3, 3)) 

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

1613 True 

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

1615 True 

1616 

1617 """ 

1618 import numpy as _nx 

1619 a, wrap = _makearray(a) 

1620 

1621 if hermitian: 

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

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

1624 # and related arrays to have the correct order 

1625 if compute_uv: 

1626 s, u = eigh(a) 

1627 sgn = sign(s) 

1628 s = abs(s) 

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

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

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

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

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

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

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

1636 else: 

1637 s = eigvalsh(a) 

1638 s = s[..., ::-1] 

1639 s = abs(s) 

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

1641 

1642 _assert_stacked_2d(a) 

1643 t, result_t = _commonType(a) 

1644 

1645 extobj = get_linalg_error_extobj(_raise_linalgerror_svd_nonconvergence) 

1646 

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

1648 if compute_uv: 

1649 if full_matrices: 

1650 if m < n: 

1651 gufunc = _umath_linalg.svd_m_f 

1652 else: 

1653 gufunc = _umath_linalg.svd_n_f 

1654 else: 

1655 if m < n: 

1656 gufunc = _umath_linalg.svd_m_s 

1657 else: 

1658 gufunc = _umath_linalg.svd_n_s 

1659 

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

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

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

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

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

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

1666 else: 

1667 if m < n: 

1668 gufunc = _umath_linalg.svd_m 

1669 else: 

1670 gufunc = _umath_linalg.svd_n 

1671 

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

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

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

1675 return s 

1676 

1677 

1678def _cond_dispatcher(x, p=None): 

1679 return (x,) 

1680 

1681 

1682@array_function_dispatch(_cond_dispatcher) 

1683def cond(x, p=None): 

1684 """ 

1685 Compute the condition number of a matrix. 

1686 

1687 This function is capable of returning the condition number using 

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

1689 Parameters below). 

1690 

1691 Parameters 

1692 ---------- 

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

1694 The matrix whose condition number is sought. 

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

1696 Order of the norm: 

1697 

1698 ===== ============================ 

1699 p norm for matrices 

1700 ===== ============================ 

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

1702 'fro' Frobenius norm 

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

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

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

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

1707 2 2-norm (largest sing. value) 

1708 -2 smallest singular value 

1709 ===== ============================ 

1710 

1711 inf means the numpy.inf object, and the Frobenius norm is 

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

1713 

1714 Returns 

1715 ------- 

1716 c : {float, inf} 

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

1718 

1719 See Also 

1720 -------- 

1721 numpy.linalg.norm 

1722 

1723 Notes 

1724 ----- 

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

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

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

1728 

1729 References 

1730 ---------- 

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

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

1733 

1734 Examples 

1735 -------- 

1736 >>> from numpy import linalg as LA 

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

1738 >>> a 

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

1740 [ 0, 1, 0], 

1741 [ 1, 0, 1]]) 

1742 >>> LA.cond(a) 

1743 1.4142135623730951 

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

1745 3.1622776601683795 

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

1747 2.0 

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

1749 1.0 

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

1751 2.0 

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

1753 1.0 

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

1755 1.4142135623730951 

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

1757 0.70710678118654746 # may vary 

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

1759 0.70710678118654746 # may vary 

1760 

1761 """ 

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

1763 if _is_empty_2d(x): 

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

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

1766 s = svd(x, compute_uv=False) 

1767 with errstate(all='ignore'): 

1768 if p == -2: 

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

1770 else: 

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

1772 else: 

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

1774 # contain nans in the entries where inversion failed. 

1775 _assert_stacked_2d(x) 

1776 _assert_stacked_square(x) 

1777 t, result_t = _commonType(x) 

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

1779 with errstate(all='ignore'): 

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

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

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

1783 

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

1785 r = asarray(r) 

1786 nan_mask = isnan(r) 

1787 if nan_mask.any(): 

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

1789 if r.ndim > 0: 

1790 r[nan_mask] = Inf 

1791 elif nan_mask: 

1792 r[()] = Inf 

1793 

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

1795 if r.ndim == 0: 

1796 r = r[()] 

1797 

1798 return r 

1799 

1800 

1801def _matrix_rank_dispatcher(M, tol=None, hermitian=None): 

1802 return (M,) 

1803 

1804 

1805@array_function_dispatch(_matrix_rank_dispatcher) 

1806def matrix_rank(M, tol=None, hermitian=False): 

1807 """ 

1808 Return matrix rank of array using SVD method 

1809 

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

1811 greater than `tol`. 

1812 

1813 .. versionchanged:: 1.14 

1814 Can now operate on stacks of matrices 

1815 

1816 Parameters 

1817 ---------- 

1818 M : {(M,), (..., M, N)} array_like 

1819 Input vector or stack of matrices. 

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

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

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

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

1824 set to ``S.max() * max(M.shape) * eps``. 

1825 

1826 .. versionchanged:: 1.14 

1827 Broadcasted against the stack of matrices 

1828 hermitian : bool, optional 

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

1830 enabling a more efficient method for finding singular values. 

1831 Defaults to False. 

1832 

1833 .. versionadded:: 1.14 

1834 

1835 Returns 

1836 ------- 

1837 rank : (...) array_like 

1838 Rank of M. 

1839 

1840 Notes 

1841 ----- 

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

1843 of the singular values of `M`. By default, we identify singular values less 

1844 than ``S.max() * max(M.shape) * eps`` as indicating rank deficiency (with 

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

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

1847 least squares [2]. 

1848 

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

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

1851 in `M` that is an exact (in floating point) linear combination of other 

1852 columns in `M`. Computing the SVD on `M` will not produce a singular value 

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

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

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

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

1857 deficiency. The threshold may declare a matrix `M` rank deficient even if 

1858 the linear combination of some columns of `M` is not exactly equal to 

1859 another column of `M` but only numerically very close to another column of 

1860 `M`. 

1861 

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

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

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

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

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

1867 

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

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

1870 sources of error in `M` that would make you consider other tolerance values 

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

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

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

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

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

1876 uncertainties are absolute rather than relative. 

1877 

1878 References 

1879 ---------- 

1880 .. [1] MATLAB reference documention, "Rank" 

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

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

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

1884 page 795. 

1885 

1886 Examples 

1887 -------- 

1888 >>> from numpy.linalg import matrix_rank 

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

1890 4 

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

1892 >>> matrix_rank(I) 

1893 3 

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

1895 1 

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

1897 0 

1898 """ 

1899 M = asarray(M) 

1900 if M.ndim < 2: 

1901 return int(not all(M==0)) 

1902 S = svd(M, compute_uv=False, hermitian=hermitian) 

1903 if tol is None: 

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

1905 else: 

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

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

1908 

1909 

1910# Generalized inverse 

1911 

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

1913 return (a,) 

1914 

1915 

1916@array_function_dispatch(_pinv_dispatcher) 

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

1918 """ 

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

1920 

1921 Calculate the generalized inverse of a matrix using its 

1922 singular-value decomposition (SVD) and including all 

1923 *large* singular values. 

1924 

1925 .. versionchanged:: 1.14 

1926 Can now operate on stacks of matrices 

1927 

1928 Parameters 

1929 ---------- 

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

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

1932 rcond : (...) array_like of float 

1933 Cutoff for small singular values. 

1934 Singular values less than or equal to 

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

1936 Broadcasts against the stack of matrices. 

1937 hermitian : bool, optional 

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

1939 enabling a more efficient method for finding singular values. 

1940 Defaults to False. 

1941 

1942 .. versionadded:: 1.17.0 

1943 

1944 Returns 

1945 ------- 

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

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

1948 is `B`. 

1949 

1950 Raises 

1951 ------ 

1952 LinAlgError 

1953 If the SVD computation does not converge. 

1954 

1955 See Also 

1956 -------- 

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

1958 scipy.linalg.pinv2 : Similar function in SciPy (SVD-based). 

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

1960 Hermitian matrix. 

1961 

1962 Notes 

1963 ----- 

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

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

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

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

1968 

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

1970 value decomposition of A, then 

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

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

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

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

1975 consisting of the reciprocals of A's singular values 

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

1977 

1978 References 

1979 ---------- 

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

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

1982 

1983 Examples 

1984 -------- 

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

1986 ``a+ * a * a+ == a+``: 

1987 

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

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

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

1991 True 

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

1993 True 

1994 

1995 """ 

1996 a, wrap = _makearray(a) 

1997 rcond = asarray(rcond) 

1998 if _is_empty_2d(a): 

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

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

2001 return wrap(res) 

2002 a = a.conjugate() 

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

2004 

2005 # discard small singular values 

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

2007 large = s > cutoff 

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

2009 s[~large] = 0 

2010 

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

2012 return wrap(res) 

2013 

2014 

2015# Determinant 

2016 

2017 

2018@array_function_dispatch(_unary_dispatcher) 

2019def slogdet(a): 

2020 """ 

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

2022 

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

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

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

2026 the determinant itself. 

2027 

2028 Parameters 

2029 ---------- 

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

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

2032 

2033 Returns 

2034 ------- 

2035 sign : (...) array_like 

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

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

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

2039 logdet : (...) array_like 

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

2041 

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

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

2044 

2045 See Also 

2046 -------- 

2047 det 

2048 

2049 Notes 

2050 ----- 

2051 

2052 .. versionadded:: 1.8.0 

2053 

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

2055 details. 

2056 

2057 .. versionadded:: 1.6.0 

2058 

2059 The determinant is computed via LU factorization using the LAPACK 

2060 routine ``z/dgetrf``. 

2061 

2062 

2063 Examples 

2064 -------- 

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

2066 

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

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

2069 >>> (sign, logdet) 

2070 (-1, 0.69314718055994529) # may vary 

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

2072 -2.0 

2073 

2074 Computing log-determinants for a stack of matrices: 

2075 

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

2077 >>> a.shape 

2078 (3, 2, 2) 

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

2080 >>> (sign, logdet) 

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

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

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

2084 

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

2086 

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

2088 0.0 

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

2090 (1, -1151.2925464970228) 

2091 

2092 """ 

2093 a = asarray(a) 

2094 _assert_stacked_2d(a) 

2095 _assert_stacked_square(a) 

2096 t, result_t = _commonType(a) 

2097 real_t = _realType(result_t) 

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

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

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

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

2102 return sign, logdet 

2103 

2104 

2105@array_function_dispatch(_unary_dispatcher) 

2106def det(a): 

2107 """ 

2108 Compute the determinant of an array. 

2109 

2110 Parameters 

2111 ---------- 

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

2113 Input array to compute determinants for. 

2114 

2115 Returns 

2116 ------- 

2117 det : (...) array_like 

2118 Determinant of `a`. 

2119 

2120 See Also 

2121 -------- 

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

2123 for large matrices where underflow/overflow may occur. 

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

2125 

2126 Notes 

2127 ----- 

2128 

2129 .. versionadded:: 1.8.0 

2130 

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

2132 details. 

2133 

2134 The determinant is computed via LU factorization using the LAPACK 

2135 routine ``z/dgetrf``. 

2136 

2137 Examples 

2138 -------- 

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

2140 

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

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

2143 -2.0 # may vary 

2144 

2145 Computing determinants for a stack of matrices: 

2146 

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

2148 >>> a.shape 

2149 (3, 2, 2) 

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

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

2152 

2153 """ 

2154 a = asarray(a) 

2155 _assert_stacked_2d(a) 

2156 _assert_stacked_square(a) 

2157 t, result_t = _commonType(a) 

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

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

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

2161 return r 

2162 

2163 

2164# Linear Least Squares 

2165 

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

2167 return (a, b) 

2168 

2169 

2170@array_function_dispatch(_lstsq_dispatcher) 

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

2172 r""" 

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

2174 

2175 Computes the vector x that approximatively solves the equation 

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

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

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

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

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

2181 Euclidean 2-norm :math:`|| b - a x ||`. 

2182 

2183 Parameters 

2184 ---------- 

2185 a : (M, N) array_like 

2186 "Coefficient" matrix. 

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

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

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

2190 of `b`. 

2191 rcond : float, optional 

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

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

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

2195 value of `a`. 

2196 

2197 .. versionchanged:: 1.14.0 

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

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

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

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

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

2203 

2204 Returns 

2205 ------- 

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

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

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

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

2210 Sums of residuals; squared Euclidean 2-norm for each column in 

2211 ``b - a*x``. 

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

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

2214 Otherwise the shape is (K,). 

2215 rank : int 

2216 Rank of matrix `a`. 

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

2218 Singular values of `a`. 

2219 

2220 Raises 

2221 ------ 

2222 LinAlgError 

2223 If computation does not converge. 

2224 

2225 See Also 

2226 -------- 

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

2228 

2229 Notes 

2230 ----- 

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

2232 

2233 Examples 

2234 -------- 

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

2236 

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

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

2239 

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

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

2242 

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

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

2245 

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

2247 >>> A 

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

2249 [ 1., 1.], 

2250 [ 2., 1.], 

2251 [ 3., 1.]]) 

2252 

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

2254 >>> m, c 

2255 (1.0 -0.95) # may vary 

2256 

2257 Plot the data along with the fitted line: 

2258 

2259 >>> import matplotlib.pyplot as plt 

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

2261 >>> _ = plt.plot(x, m*x + c, 'r', label='Fitted line') 

2262 >>> _ = plt.legend() 

2263 >>> plt.show() 

2264 

2265 """ 

2266 a, _ = _makearray(a) 

2267 b, wrap = _makearray(b) 

2268 is_1d = b.ndim == 1 

2269 if is_1d: 

2270 b = b[:, newaxis] 

2271 _assert_2d(a, b) 

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

2273 m2, n_rhs = b.shape[-2:] 

2274 if m != m2: 

2275 raise LinAlgError('Incompatible dimensions') 

2276 

2277 t, result_t = _commonType(a, b) 

2278 # FIXME: real_t is unused 

2279 real_t = _linalgRealType(t) 

2280 result_real_t = _realType(result_t) 

2281 

2282 # Determine default rcond value 

2283 if rcond == "warn": 

2284 # 2017-08-19, 1.14.0 

2285 warnings.warn("`rcond` parameter will change to the default of " 

2286 "machine precision times ``max(M, N)`` where M and N " 

2287 "are the input matrix dimensions.\n" 

2288 "To use the future default and silence this warning " 

2289 "we advise to pass `rcond=None`, to keep using the old, " 

2290 "explicitly pass `rcond=-1`.", 

2291 FutureWarning, stacklevel=3) 

2292 rcond = -1 

2293 if rcond is None: 

2294 rcond = finfo(t).eps * max(n, m) 

2295 

2296 if m <= n: 

2297 gufunc = _umath_linalg.lstsq_m 

2298 else: 

2299 gufunc = _umath_linalg.lstsq_n 

2300 

2301 signature = 'DDd->Ddid' if isComplexType(t) else 'ddd->ddid' 

2302 extobj = get_linalg_error_extobj(_raise_linalgerror_lstsq) 

2303 if n_rhs == 0: 

2304 # lapack can't handle n_rhs = 0 - so allocate the array one larger in that axis 

2305 b = zeros(b.shape[:-2] + (m, n_rhs + 1), dtype=b.dtype) 

2306 x, resids, rank, s = gufunc(a, b, rcond, signature=signature, extobj=extobj) 

2307 if m == 0: 

2308 x[...] = 0 

2309 if n_rhs == 0: 

2310 # remove the item we added 

2311 x = x[..., :n_rhs] 

2312 resids = resids[..., :n_rhs] 

2313 

2314 # remove the axis we added 

2315 if is_1d: 

2316 x = x.squeeze(axis=-1) 

2317 # we probably should squeeze resids too, but we can't 

2318 # without breaking compatibility. 

2319 

2320 # as documented 

2321 if rank != n or m <= n: 

2322 resids = array([], result_real_t) 

2323 

2324 # coerce output arrays 

2325 s = s.astype(result_real_t, copy=False) 

2326 resids = resids.astype(result_real_t, copy=False) 

2327 x = x.astype(result_t, copy=True) # Copying lets the memory in r_parts be freed 

2328 return wrap(x), wrap(resids), rank, s 

2329 

2330 

2331def _multi_svd_norm(x, row_axis, col_axis, op): 

2332 """Compute a function of the singular values of the 2-D matrices in `x`. 

2333 

2334 This is a private utility function used by `numpy.linalg.norm()`. 

2335 

2336 Parameters 

2337 ---------- 

2338 x : ndarray 

2339 row_axis, col_axis : int 

2340 The axes of `x` that hold the 2-D matrices. 

2341 op : callable 

2342 This should be either numpy.amin or `numpy.amax` or `numpy.sum`. 

2343 

2344 Returns 

2345 ------- 

2346 result : float or ndarray 

2347 If `x` is 2-D, the return values is a float. 

2348 Otherwise, it is an array with ``x.ndim - 2`` dimensions. 

2349 The return values are either the minimum or maximum or sum of the 

2350 singular values of the matrices, depending on whether `op` 

2351 is `numpy.amin` or `numpy.amax` or `numpy.sum`. 

2352 

2353 """ 

2354 y = moveaxis(x, (row_axis, col_axis), (-2, -1)) 

2355 result = op(svd(y, compute_uv=False), axis=-1) 

2356 return result 

2357 

2358 

2359def _norm_dispatcher(x, ord=None, axis=None, keepdims=None): 

2360 return (x,) 

2361 

2362 

2363@array_function_dispatch(_norm_dispatcher) 

2364def norm(x, ord=None, axis=None, keepdims=False): 

2365 """ 

2366 Matrix or vector norm. 

2367 

2368 This function is able to return one of eight different matrix norms, 

2369 or one of an infinite number of vector norms (described below), depending 

2370 on the value of the ``ord`` parameter. 

2371 

2372 Parameters 

2373 ---------- 

2374 x : array_like 

2375 Input array. If `axis` is None, `x` must be 1-D or 2-D, unless `ord` 

2376 is None. If both `axis` and `ord` are None, the 2-norm of 

2377 ``x.ravel`` will be returned. 

2378 ord : {non-zero int, inf, -inf, 'fro', 'nuc'}, optional 

2379 Order of the norm (see table under ``Notes``). inf means numpy's 

2380 `inf` object. The default is None. 

2381 axis : {None, int, 2-tuple of ints}, optional. 

2382 If `axis` is an integer, it specifies the axis of `x` along which to 

2383 compute the vector norms. If `axis` is a 2-tuple, it specifies the 

2384 axes that hold 2-D matrices, and the matrix norms of these matrices 

2385 are computed. If `axis` is None then either a vector norm (when `x` 

2386 is 1-D) or a matrix norm (when `x` is 2-D) is returned. The default 

2387 is None. 

2388 

2389 .. versionadded:: 1.8.0 

2390 

2391 keepdims : bool, optional 

2392 If this is set to True, the axes which are normed over are left in the 

2393 result as dimensions with size one. With this option the result will 

2394 broadcast correctly against the original `x`. 

2395 

2396 .. versionadded:: 1.10.0 

2397 

2398 Returns 

2399 ------- 

2400 n : float or ndarray 

2401 Norm of the matrix or vector(s). 

2402 

2403 See Also 

2404 -------- 

2405 scipy.linalg.norm : Similar function in SciPy. 

2406 

2407 Notes 

2408 ----- 

2409 For values of ``ord < 1``, the result is, strictly speaking, not a 

2410 mathematical 'norm', but it may still be useful for various numerical 

2411 purposes. 

2412 

2413 The following norms can be calculated: 

2414 

2415 ===== ============================ ========================== 

2416 ord norm for matrices norm for vectors 

2417 ===== ============================ ========================== 

2418 None Frobenius norm 2-norm 

2419 'fro' Frobenius norm -- 

2420 'nuc' nuclear norm -- 

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

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

2423 0 -- sum(x != 0) 

2424 1 max(sum(abs(x), axis=0)) as below 

2425 -1 min(sum(abs(x), axis=0)) as below 

2426 2 2-norm (largest sing. value) as below 

2427 -2 smallest singular value as below 

2428 other -- sum(abs(x)**ord)**(1./ord) 

2429 ===== ============================ ========================== 

2430 

2431 The Frobenius norm is given by [1]_: 

2432 

2433 :math:`||A||_F = [\\sum_{i,j} abs(a_{i,j})^2]^{1/2}` 

2434 

2435 The nuclear norm is the sum of the singular values. 

2436 

2437 Both the Frobenius and nuclear norm orders are only defined for 

2438 matrices and raise a ValueError when ``x.ndim != 2``. 

2439 

2440 References 

2441 ---------- 

2442 .. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*, 

2443 Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15 

2444 

2445 Examples 

2446 -------- 

2447 >>> from numpy import linalg as LA 

2448 >>> a = np.arange(9) - 4 

2449 >>> a 

2450 array([-4, -3, -2, ..., 2, 3, 4]) 

2451 >>> b = a.reshape((3, 3)) 

2452 >>> b 

2453 array([[-4, -3, -2], 

2454 [-1, 0, 1], 

2455 [ 2, 3, 4]]) 

2456 

2457 >>> LA.norm(a) 

2458 7.745966692414834 

2459 >>> LA.norm(b) 

2460 7.745966692414834 

2461 >>> LA.norm(b, 'fro') 

2462 7.745966692414834 

2463 >>> LA.norm(a, np.inf) 

2464 4.0 

2465 >>> LA.norm(b, np.inf) 

2466 9.0 

2467 >>> LA.norm(a, -np.inf) 

2468 0.0 

2469 >>> LA.norm(b, -np.inf) 

2470 2.0 

2471 

2472 >>> LA.norm(a, 1) 

2473 20.0 

2474 >>> LA.norm(b, 1) 

2475 7.0 

2476 >>> LA.norm(a, -1) 

2477 -4.6566128774142013e-010 

2478 >>> LA.norm(b, -1) 

2479 6.0 

2480 >>> LA.norm(a, 2) 

2481 7.745966692414834 

2482 >>> LA.norm(b, 2) 

2483 7.3484692283495345 

2484 

2485 >>> LA.norm(a, -2) 

2486 0.0 

2487 >>> LA.norm(b, -2) 

2488 1.8570331885190563e-016 # may vary 

2489 >>> LA.norm(a, 3) 

2490 5.8480354764257312 # may vary 

2491 >>> LA.norm(a, -3) 

2492 0.0 

2493 

2494 Using the `axis` argument to compute vector norms: 

2495 

2496 >>> c = np.array([[ 1, 2, 3], 

2497 ... [-1, 1, 4]]) 

2498 >>> LA.norm(c, axis=0) 

2499 array([ 1.41421356, 2.23606798, 5. ]) 

2500 >>> LA.norm(c, axis=1) 

2501 array([ 3.74165739, 4.24264069]) 

2502 >>> LA.norm(c, ord=1, axis=1) 

2503 array([ 6., 6.]) 

2504 

2505 Using the `axis` argument to compute matrix norms: 

2506 

2507 >>> m = np.arange(8).reshape(2,2,2) 

2508 >>> LA.norm(m, axis=(1,2)) 

2509 array([ 3.74165739, 11.22497216]) 

2510 >>> LA.norm(m[0, :, :]), LA.norm(m[1, :, :]) 

2511 (3.7416573867739413, 11.224972160321824) 

2512 

2513 """ 

2514 x = asarray(x) 

2515 

2516 if not issubclass(x.dtype.type, (inexact, object_)): 

2517 x = x.astype(float) 

2518 

2519 # Immediately handle some default, simple, fast, and common cases. 

2520 if axis is None: 

2521 ndim = x.ndim 

2522 if ((ord is None) or 

2523 (ord in ('f', 'fro') and ndim == 2) or 

2524 (ord == 2 and ndim == 1)): 

2525 

2526 x = x.ravel(order='K') 

2527 if isComplexType(x.dtype.type): 

2528 sqnorm = dot(x.real, x.real) + dot(x.imag, x.imag) 

2529 else: 

2530 sqnorm = dot(x, x) 

2531 ret = sqrt(sqnorm) 

2532 if keepdims: 

2533 ret = ret.reshape(ndim*[1]) 

2534 return ret 

2535 

2536 # Normalize the `axis` argument to a tuple. 

2537 nd = x.ndim 

2538 if axis is None: 

2539 axis = tuple(range(nd)) 

2540 elif not isinstance(axis, tuple): 

2541 try: 

2542 axis = int(axis) 

2543 except Exception as e: 

2544 raise TypeError("'axis' must be None, an integer or a tuple of integers") from e 

2545 axis = (axis,) 

2546 

2547 if len(axis) == 1: 

2548 if ord == Inf: 

2549 return abs(x).max(axis=axis, keepdims=keepdims) 

2550 elif ord == -Inf: 

2551 return abs(x).min(axis=axis, keepdims=keepdims) 

2552 elif ord == 0: 

2553 # Zero norm 

2554 return (x != 0).astype(x.real.dtype).sum(axis=axis, keepdims=keepdims) 

2555 elif ord == 1: 

2556 # special case for speedup 

2557 return add.reduce(abs(x), axis=axis, keepdims=keepdims) 

2558 elif ord is None or ord == 2: 

2559 # special case for speedup 

2560 s = (x.conj() * x).real 

2561 return sqrt(add.reduce(s, axis=axis, keepdims=keepdims)) 

2562 # None of the str-type keywords for ord ('fro', 'nuc')  

2563 # are valid for vectors 

2564 elif isinstance(ord, str): 

2565 raise ValueError(f"Invalid norm order '{ord}' for vectors") 

2566 else: 

2567 absx = abs(x) 

2568 absx **= ord 

2569 ret = add.reduce(absx, axis=axis, keepdims=keepdims) 

2570 ret **= (1 / ord) 

2571 return ret 

2572 elif len(axis) == 2: 

2573 row_axis, col_axis = axis 

2574 row_axis = normalize_axis_index(row_axis, nd) 

2575 col_axis = normalize_axis_index(col_axis, nd) 

2576 if row_axis == col_axis: 

2577 raise ValueError('Duplicate axes given.') 

2578 if ord == 2: 

2579 ret = _multi_svd_norm(x, row_axis, col_axis, amax) 

2580 elif ord == -2: 

2581 ret = _multi_svd_norm(x, row_axis, col_axis, amin) 

2582 elif ord == 1: 

2583 if col_axis > row_axis: 

2584 col_axis -= 1 

2585 ret = add.reduce(abs(x), axis=row_axis).max(axis=col_axis) 

2586 elif ord == Inf: 

2587 if row_axis > col_axis: 

2588 row_axis -= 1 

2589 ret = add.reduce(abs(x), axis=col_axis).max(axis=row_axis) 

2590 elif ord == -1: 

2591 if col_axis > row_axis: 

2592 col_axis -= 1 

2593 ret = add.reduce(abs(x), axis=row_axis).min(axis=col_axis) 

2594 elif ord == -Inf: 

2595 if row_axis > col_axis: 

2596 row_axis -= 1 

2597 ret = add.reduce(abs(x), axis=col_axis).min(axis=row_axis) 

2598 elif ord in [None, 'fro', 'f']: 

2599 ret = sqrt(add.reduce((x.conj() * x).real, axis=axis)) 

2600 elif ord == 'nuc': 

2601 ret = _multi_svd_norm(x, row_axis, col_axis, sum) 

2602 else: 

2603 raise ValueError("Invalid norm order for matrices.") 

2604 if keepdims: 

2605 ret_shape = list(x.shape) 

2606 ret_shape[axis[0]] = 1 

2607 ret_shape[axis[1]] = 1 

2608 ret = ret.reshape(ret_shape) 

2609 return ret 

2610 else: 

2611 raise ValueError("Improper number of dimensions to norm.") 

2612 

2613 

2614# multi_dot 

2615 

2616def _multidot_dispatcher(arrays, *, out=None): 

2617 yield from arrays 

2618 yield out 

2619 

2620 

2621@array_function_dispatch(_multidot_dispatcher) 

2622def multi_dot(arrays, *, out=None): 

2623 """ 

2624 Compute the dot product of two or more arrays in a single function call, 

2625 while automatically selecting the fastest evaluation order. 

2626 

2627 `multi_dot` chains `numpy.dot` and uses optimal parenthesization 

2628 of the matrices [1]_ [2]_. Depending on the shapes of the matrices, 

2629 this can speed up the multiplication a lot. 

2630 

2631 If the first argument is 1-D it is treated as a row vector. 

2632 If the last argument is 1-D it is treated as a column vector. 

2633 The other arguments must be 2-D. 

2634 

2635 Think of `multi_dot` as:: 

2636 

2637 def multi_dot(arrays): return functools.reduce(np.dot, arrays) 

2638 

2639 

2640 Parameters 

2641 ---------- 

2642 arrays : sequence of array_like 

2643 If the first argument is 1-D it is treated as row vector. 

2644 If the last argument is 1-D it is treated as column vector. 

2645 The other arguments must be 2-D. 

2646 out : ndarray, optional 

2647 Output argument. This must have the exact kind that would be returned 

2648 if it was not used. In particular, it must have the right type, must be 

2649 C-contiguous, and its dtype must be the dtype that would be returned 

2650 for `dot(a, b)`. This is a performance feature. Therefore, if these 

2651 conditions are not met, an exception is raised, instead of attempting 

2652 to be flexible. 

2653 

2654 .. versionadded:: 1.19.0 

2655 

2656 Returns 

2657 ------- 

2658 output : ndarray 

2659 Returns the dot product of the supplied arrays. 

2660 

2661 See Also 

2662 -------- 

2663 dot : dot multiplication with two arguments. 

2664 

2665 References 

2666 ---------- 

2667 

2668 .. [1] Cormen, "Introduction to Algorithms", Chapter 15.2, p. 370-378 

2669 .. [2] https://en.wikipedia.org/wiki/Matrix_chain_multiplication 

2670 

2671 Examples 

2672 -------- 

2673 `multi_dot` allows you to write:: 

2674 

2675 >>> from numpy.linalg import multi_dot 

2676 >>> # Prepare some data 

2677 >>> A = np.random.random((10000, 100)) 

2678 >>> B = np.random.random((100, 1000)) 

2679 >>> C = np.random.random((1000, 5)) 

2680 >>> D = np.random.random((5, 333)) 

2681 >>> # the actual dot multiplication 

2682 >>> _ = multi_dot([A, B, C, D]) 

2683 

2684 instead of:: 

2685 

2686 >>> _ = np.dot(np.dot(np.dot(A, B), C), D) 

2687 >>> # or 

2688 >>> _ = A.dot(B).dot(C).dot(D) 

2689 

2690 Notes 

2691 ----- 

2692 The cost for a matrix multiplication can be calculated with the 

2693 following function:: 

2694 

2695 def cost(A, B): 

2696 return A.shape[0] * A.shape[1] * B.shape[1] 

2697 

2698 Assume we have three matrices 

2699 :math:`A_{10x100}, B_{100x5}, C_{5x50}`. 

2700 

2701 The costs for the two different parenthesizations are as follows:: 

2702 

2703 cost((AB)C) = 10*100*5 + 10*5*50 = 5000 + 2500 = 7500 

2704 cost(A(BC)) = 10*100*50 + 100*5*50 = 50000 + 25000 = 75000 

2705 

2706 """ 

2707 n = len(arrays) 

2708 # optimization only makes sense for len(arrays) > 2 

2709 if n < 2: 

2710 raise ValueError("Expecting at least two arrays.") 

2711 elif n == 2: 

2712 return dot(arrays[0], arrays[1], out=out) 

2713 

2714 arrays = [asanyarray(a) for a in arrays] 

2715 

2716 # save original ndim to reshape the result array into the proper form later 

2717 ndim_first, ndim_last = arrays[0].ndim, arrays[-1].ndim 

2718 # Explicitly convert vectors to 2D arrays to keep the logic of the internal 

2719 # _multi_dot_* functions as simple as possible. 

2720 if arrays[0].ndim == 1: 

2721 arrays[0] = atleast_2d(arrays[0]) 

2722 if arrays[-1].ndim == 1: 

2723 arrays[-1] = atleast_2d(arrays[-1]).T 

2724 _assert_2d(*arrays) 

2725 

2726 # _multi_dot_three is much faster than _multi_dot_matrix_chain_order 

2727 if n == 3: 

2728 result = _multi_dot_three(arrays[0], arrays[1], arrays[2], out=out) 

2729 else: 

2730 order = _multi_dot_matrix_chain_order(arrays) 

2731 result = _multi_dot(arrays, order, 0, n - 1, out=out) 

2732 

2733 # return proper shape 

2734 if ndim_first == 1 and ndim_last == 1: 

2735 return result[0, 0] # scalar 

2736 elif ndim_first == 1 or ndim_last == 1: 

2737 return result.ravel() # 1-D 

2738 else: 

2739 return result 

2740 

2741 

2742def _multi_dot_three(A, B, C, out=None): 

2743 """ 

2744 Find the best order for three arrays and do the multiplication. 

2745 

2746 For three arguments `_multi_dot_three` is approximately 15 times faster 

2747 than `_multi_dot_matrix_chain_order` 

2748 

2749 """ 

2750 a0, a1b0 = A.shape 

2751 b1c0, c1 = C.shape 

2752 # cost1 = cost((AB)C) = a0*a1b0*b1c0 + a0*b1c0*c1 

2753 cost1 = a0 * b1c0 * (a1b0 + c1) 

2754 # cost2 = cost(A(BC)) = a1b0*b1c0*c1 + a0*a1b0*c1 

2755 cost2 = a1b0 * c1 * (a0 + b1c0) 

2756 

2757 if cost1 < cost2: 

2758 return dot(dot(A, B), C, out=out) 

2759 else: 

2760 return dot(A, dot(B, C), out=out) 

2761 

2762 

2763def _multi_dot_matrix_chain_order(arrays, return_costs=False): 

2764 """ 

2765 Return a np.array that encodes the optimal order of mutiplications. 

2766 

2767 The optimal order array is then used by `_multi_dot()` to do the 

2768 multiplication. 

2769 

2770 Also return the cost matrix if `return_costs` is `True` 

2771 

2772 The implementation CLOSELY follows Cormen, "Introduction to Algorithms", 

2773 Chapter 15.2, p. 370-378. Note that Cormen uses 1-based indices. 

2774 

2775 cost[i, j] = min([ 

2776 cost[prefix] + cost[suffix] + cost_mult(prefix, suffix) 

2777 for k in range(i, j)]) 

2778 

2779 """ 

2780 n = len(arrays) 

2781 # p stores the dimensions of the matrices 

2782 # Example for p: A_{10x100}, B_{100x5}, C_{5x50} --> p = [10, 100, 5, 50] 

2783 p = [a.shape[0] for a in arrays] + [arrays[-1].shape[1]] 

2784 # m is a matrix of costs of the subproblems 

2785 # m[i,j]: min number of scalar multiplications needed to compute A_{i..j} 

2786 m = zeros((n, n), dtype=double) 

2787 # s is the actual ordering 

2788 # s[i, j] is the value of k at which we split the product A_i..A_j 

2789 s = empty((n, n), dtype=intp) 

2790 

2791 for l in range(1, n): 

2792 for i in range(n - l): 

2793 j = i + l 

2794 m[i, j] = Inf 

2795 for k in range(i, j): 

2796 q = m[i, k] + m[k+1, j] + p[i]*p[k+1]*p[j+1] 

2797 if q < m[i, j]: 

2798 m[i, j] = q 

2799 s[i, j] = k # Note that Cormen uses 1-based index 

2800 

2801 return (s, m) if return_costs else s 

2802 

2803 

2804def _multi_dot(arrays, order, i, j, out=None): 

2805 """Actually do the multiplication with the given order.""" 

2806 if i == j: 

2807 # the initial call with non-None out should never get here 

2808 assert out is None 

2809 

2810 return arrays[i] 

2811 else: 

2812 return dot(_multi_dot(arrays, order, i, order[i, j]), 

2813 _multi_dot(arrays, order, order[i, j] + 1, j), 

2814 out=out)