Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/numpy/linalg/linalg.py: 18%
669 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-04-03 06:39 +0000
1"""Lite version of scipy.linalg.
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"""
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']
17import functools
18import operator
19import warnings
20from typing import NamedTuple, Any
22from .._utils import set_module
23from numpy.core import (
24 array, asarray, zeros, empty, empty_like, intc, single, double,
25 csingle, cdouble, inexact, complexfloating, newaxis, all, Inf, dot,
26 add, multiply, sqrt, sum, isfinite,
27 finfo, errstate, geterrobj, moveaxis, amin, amax, prod, abs,
28 atleast_2d, intp, asanyarray, object_, matmul,
29 swapaxes, divide, count_nonzero, isnan, sign, argsort, sort,
30 reciprocal
31)
32from numpy.core.multiarray import normalize_axis_index
33from numpy.core import overrides
34from numpy.lib.twodim_base import triu, eye
35from numpy.linalg import _umath_linalg
37from numpy._typing import NDArray
39class EigResult(NamedTuple):
40 eigenvalues: NDArray[Any]
41 eigenvectors: NDArray[Any]
43class EighResult(NamedTuple):
44 eigenvalues: NDArray[Any]
45 eigenvectors: NDArray[Any]
47class QRResult(NamedTuple):
48 Q: NDArray[Any]
49 R: NDArray[Any]
51class SlogdetResult(NamedTuple):
52 sign: NDArray[Any]
53 logabsdet: NDArray[Any]
55class SVDResult(NamedTuple):
56 U: NDArray[Any]
57 S: NDArray[Any]
58 Vh: NDArray[Any]
60array_function_dispatch = functools.partial(
61 overrides.array_function_dispatch, module='numpy.linalg')
64fortran_int = intc
67@set_module('numpy.linalg')
68class LinAlgError(ValueError):
69 """
70 Generic Python-exception-derived object raised by linalg functions.
72 General purpose exception class, derived from Python's ValueError
73 class, programmatically raised in linalg functions when a Linear
74 Algebra-related condition would prevent further correct execution of the
75 function.
77 Parameters
78 ----------
79 None
81 Examples
82 --------
83 >>> from numpy import linalg as LA
84 >>> LA.inv(np.zeros((2,2)))
85 Traceback (most recent call last):
86 File "<stdin>", line 1, in <module>
87 File "...linalg.py", line 350,
88 in inv return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))
89 File "...linalg.py", line 249,
90 in solve
91 raise LinAlgError('Singular matrix')
92 numpy.linalg.LinAlgError: Singular matrix
94 """
97def _determine_error_states():
98 errobj = geterrobj()
99 bufsize = errobj[0]
101 with errstate(invalid='call', over='ignore',
102 divide='ignore', under='ignore'):
103 invalid_call_errmask = geterrobj()[1]
105 return [bufsize, invalid_call_errmask, None]
107# Dealing with errors in _umath_linalg
108_linalg_error_extobj = _determine_error_states()
109del _determine_error_states
111def _raise_linalgerror_singular(err, flag):
112 raise LinAlgError("Singular matrix")
114def _raise_linalgerror_nonposdef(err, flag):
115 raise LinAlgError("Matrix is not positive definite")
117def _raise_linalgerror_eigenvalues_nonconvergence(err, flag):
118 raise LinAlgError("Eigenvalues did not converge")
120def _raise_linalgerror_svd_nonconvergence(err, flag):
121 raise LinAlgError("SVD did not converge")
123def _raise_linalgerror_lstsq(err, flag):
124 raise LinAlgError("SVD did not converge in Linear Least Squares")
126def _raise_linalgerror_qr(err, flag):
127 raise LinAlgError("Incorrect argument found while performing "
128 "QR factorization")
130def get_linalg_error_extobj(callback):
131 extobj = list(_linalg_error_extobj) # make a copy
132 extobj[2] = callback
133 return extobj
135def _makearray(a):
136 new = asarray(a)
137 wrap = getattr(a, "__array_prepare__", new.__array_wrap__)
138 return new, wrap
140def isComplexType(t):
141 return issubclass(t, complexfloating)
143_real_types_map = {single : single,
144 double : double,
145 csingle : single,
146 cdouble : double}
148_complex_types_map = {single : csingle,
149 double : cdouble,
150 csingle : csingle,
151 cdouble : cdouble}
153def _realType(t, default=double):
154 return _real_types_map.get(t, default)
156def _complexType(t, default=cdouble):
157 return _complex_types_map.get(t, default)
159def _commonType(*arrays):
160 # in lite version, use higher precision (always double or cdouble)
161 result_type = single
162 is_complex = False
163 for a in arrays:
164 type_ = a.dtype.type
165 if issubclass(type_, inexact):
166 if isComplexType(type_):
167 is_complex = True
168 rt = _realType(type_, default=None)
169 if rt is double:
170 result_type = double
171 elif rt is None:
172 # unsupported inexact scalar
173 raise TypeError("array type %s is unsupported in linalg" %
174 (a.dtype.name,))
175 else:
176 result_type = double
177 if is_complex:
178 result_type = _complex_types_map[result_type]
179 return cdouble, result_type
180 else:
181 return double, result_type
184def _to_native_byte_order(*arrays):
185 ret = []
186 for arr in arrays:
187 if arr.dtype.byteorder not in ('=', '|'):
188 ret.append(asarray(arr, dtype=arr.dtype.newbyteorder('=')))
189 else:
190 ret.append(arr)
191 if len(ret) == 1:
192 return ret[0]
193 else:
194 return ret
197def _assert_2d(*arrays):
198 for a in arrays:
199 if a.ndim != 2:
200 raise LinAlgError('%d-dimensional array given. Array must be '
201 'two-dimensional' % a.ndim)
203def _assert_stacked_2d(*arrays):
204 for a in arrays:
205 if a.ndim < 2:
206 raise LinAlgError('%d-dimensional array given. Array must be '
207 'at least two-dimensional' % a.ndim)
209def _assert_stacked_square(*arrays):
210 for a in arrays:
211 m, n = a.shape[-2:]
212 if m != n:
213 raise LinAlgError('Last 2 dimensions of the array must be square')
215def _assert_finite(*arrays):
216 for a in arrays:
217 if not isfinite(a).all():
218 raise LinAlgError("Array must not contain infs or NaNs")
220def _is_empty_2d(arr):
221 # check size first for efficiency
222 return arr.size == 0 and prod(arr.shape[-2:]) == 0
225def transpose(a):
226 """
227 Transpose each matrix in a stack of matrices.
229 Unlike np.transpose, this only swaps the last two axes, rather than all of
230 them
232 Parameters
233 ----------
234 a : (...,M,N) array_like
236 Returns
237 -------
238 aT : (...,N,M) ndarray
239 """
240 return swapaxes(a, -1, -2)
242# Linear equations
244def _tensorsolve_dispatcher(a, b, axes=None):
245 return (a, b)
248@array_function_dispatch(_tensorsolve_dispatcher)
249def tensorsolve(a, b, axes=None):
250 """
251 Solve the tensor equation ``a x = b`` for x.
253 It is assumed that all indices of `x` are summed over in the product,
254 together with the rightmost indices of `a`, as is done in, for example,
255 ``tensordot(a, x, axes=x.ndim)``.
257 Parameters
258 ----------
259 a : array_like
260 Coefficient tensor, of shape ``b.shape + Q``. `Q`, a tuple, equals
261 the shape of that sub-tensor of `a` consisting of the appropriate
262 number of its rightmost indices, and must be such that
263 ``prod(Q) == prod(b.shape)`` (in which sense `a` is said to be
264 'square').
265 b : array_like
266 Right-hand tensor, which can be of any shape.
267 axes : tuple of ints, optional
268 Axes in `a` to reorder to the right, before inversion.
269 If None (default), no reordering is done.
271 Returns
272 -------
273 x : ndarray, shape Q
275 Raises
276 ------
277 LinAlgError
278 If `a` is singular or not 'square' (in the above sense).
280 See Also
281 --------
282 numpy.tensordot, tensorinv, numpy.einsum
284 Examples
285 --------
286 >>> a = np.eye(2*3*4)
287 >>> a.shape = (2*3, 4, 2, 3, 4)
288 >>> b = np.random.randn(2*3, 4)
289 >>> x = np.linalg.tensorsolve(a, b)
290 >>> x.shape
291 (2, 3, 4)
292 >>> np.allclose(np.tensordot(a, x, axes=3), b)
293 True
295 """
296 a, wrap = _makearray(a)
297 b = asarray(b)
298 an = a.ndim
300 if axes is not None:
301 allaxes = list(range(0, an))
302 for k in axes:
303 allaxes.remove(k)
304 allaxes.insert(an, k)
305 a = a.transpose(allaxes)
307 oldshape = a.shape[-(an-b.ndim):]
308 prod = 1
309 for k in oldshape:
310 prod *= k
312 if a.size != prod ** 2:
313 raise LinAlgError(
314 "Input arrays must satisfy the requirement \
315 prod(a.shape[b.ndim:]) == prod(a.shape[:b.ndim])"
316 )
318 a = a.reshape(prod, prod)
319 b = b.ravel()
320 res = wrap(solve(a, b))
321 res.shape = oldshape
322 return res
325def _solve_dispatcher(a, b):
326 return (a, b)
329@array_function_dispatch(_solve_dispatcher)
330def solve(a, b):
331 """
332 Solve a linear matrix equation, or system of linear scalar equations.
334 Computes the "exact" solution, `x`, of the well-determined, i.e., full
335 rank, linear matrix equation `ax = b`.
337 Parameters
338 ----------
339 a : (..., M, M) array_like
340 Coefficient matrix.
341 b : {(..., M,), (..., M, K)}, array_like
342 Ordinate or "dependent variable" values.
344 Returns
345 -------
346 x : {(..., M,), (..., M, K)} ndarray
347 Solution to the system a x = b. Returned shape is identical to `b`.
349 Raises
350 ------
351 LinAlgError
352 If `a` is singular or not square.
354 See Also
355 --------
356 scipy.linalg.solve : Similar function in SciPy.
358 Notes
359 -----
361 .. versionadded:: 1.8.0
363 Broadcasting rules apply, see the `numpy.linalg` documentation for
364 details.
366 The solutions are computed using LAPACK routine ``_gesv``.
368 `a` must be square and of full-rank, i.e., all rows (or, equivalently,
369 columns) must be linearly independent; if either is not true, use
370 `lstsq` for the least-squares best "solution" of the
371 system/equation.
373 References
374 ----------
375 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
376 FL, Academic Press, Inc., 1980, pg. 22.
378 Examples
379 --------
380 Solve the system of equations ``x0 + 2 * x1 = 1`` and ``3 * x0 + 5 * x1 = 2``:
382 >>> a = np.array([[1, 2], [3, 5]])
383 >>> b = np.array([1, 2])
384 >>> x = np.linalg.solve(a, b)
385 >>> x
386 array([-1., 1.])
388 Check that the solution is correct:
390 >>> np.allclose(np.dot(a, x), b)
391 True
393 """
394 a, _ = _makearray(a)
395 _assert_stacked_2d(a)
396 _assert_stacked_square(a)
397 b, wrap = _makearray(b)
398 t, result_t = _commonType(a, b)
400 # We use the b = (..., M,) logic, only if the number of extra dimensions
401 # match exactly
402 if b.ndim == a.ndim - 1:
403 gufunc = _umath_linalg.solve1
404 else:
405 gufunc = _umath_linalg.solve
407 signature = 'DD->D' if isComplexType(t) else 'dd->d'
408 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
409 r = gufunc(a, b, signature=signature, extobj=extobj)
411 return wrap(r.astype(result_t, copy=False))
414def _tensorinv_dispatcher(a, ind=None):
415 return (a,)
418@array_function_dispatch(_tensorinv_dispatcher)
419def tensorinv(a, ind=2):
420 """
421 Compute the 'inverse' of an N-dimensional array.
423 The result is an inverse for `a` relative to the tensordot operation
424 ``tensordot(a, b, ind)``, i. e., up to floating-point accuracy,
425 ``tensordot(tensorinv(a), a, ind)`` is the "identity" tensor for the
426 tensordot operation.
428 Parameters
429 ----------
430 a : array_like
431 Tensor to 'invert'. Its shape must be 'square', i. e.,
432 ``prod(a.shape[:ind]) == prod(a.shape[ind:])``.
433 ind : int, optional
434 Number of first indices that are involved in the inverse sum.
435 Must be a positive integer, default is 2.
437 Returns
438 -------
439 b : ndarray
440 `a`'s tensordot inverse, shape ``a.shape[ind:] + a.shape[:ind]``.
442 Raises
443 ------
444 LinAlgError
445 If `a` is singular or not 'square' (in the above sense).
447 See Also
448 --------
449 numpy.tensordot, tensorsolve
451 Examples
452 --------
453 >>> a = np.eye(4*6)
454 >>> a.shape = (4, 6, 8, 3)
455 >>> ainv = np.linalg.tensorinv(a, ind=2)
456 >>> ainv.shape
457 (8, 3, 4, 6)
458 >>> b = np.random.randn(4, 6)
459 >>> np.allclose(np.tensordot(ainv, b), np.linalg.tensorsolve(a, b))
460 True
462 >>> a = np.eye(4*6)
463 >>> a.shape = (24, 8, 3)
464 >>> ainv = np.linalg.tensorinv(a, ind=1)
465 >>> ainv.shape
466 (8, 3, 24)
467 >>> b = np.random.randn(24)
468 >>> np.allclose(np.tensordot(ainv, b, 1), np.linalg.tensorsolve(a, b))
469 True
471 """
472 a = asarray(a)
473 oldshape = a.shape
474 prod = 1
475 if ind > 0:
476 invshape = oldshape[ind:] + oldshape[:ind]
477 for k in oldshape[ind:]:
478 prod *= k
479 else:
480 raise ValueError("Invalid ind argument.")
481 a = a.reshape(prod, -1)
482 ia = inv(a)
483 return ia.reshape(*invshape)
486# Matrix inversion
488def _unary_dispatcher(a):
489 return (a,)
492@array_function_dispatch(_unary_dispatcher)
493def inv(a):
494 """
495 Compute the (multiplicative) inverse of a matrix.
497 Given a square matrix `a`, return the matrix `ainv` satisfying
498 ``dot(a, ainv) = dot(ainv, a) = eye(a.shape[0])``.
500 Parameters
501 ----------
502 a : (..., M, M) array_like
503 Matrix to be inverted.
505 Returns
506 -------
507 ainv : (..., M, M) ndarray or matrix
508 (Multiplicative) inverse of the matrix `a`.
510 Raises
511 ------
512 LinAlgError
513 If `a` is not square or inversion fails.
515 See Also
516 --------
517 scipy.linalg.inv : Similar function in SciPy.
519 Notes
520 -----
522 .. versionadded:: 1.8.0
524 Broadcasting rules apply, see the `numpy.linalg` documentation for
525 details.
527 Examples
528 --------
529 >>> from numpy.linalg import inv
530 >>> a = np.array([[1., 2.], [3., 4.]])
531 >>> ainv = inv(a)
532 >>> np.allclose(np.dot(a, ainv), np.eye(2))
533 True
534 >>> np.allclose(np.dot(ainv, a), np.eye(2))
535 True
537 If a is a matrix object, then the return value is a matrix as well:
539 >>> ainv = inv(np.matrix(a))
540 >>> ainv
541 matrix([[-2. , 1. ],
542 [ 1.5, -0.5]])
544 Inverses of several matrices can be computed at once:
546 >>> a = np.array([[[1., 2.], [3., 4.]], [[1, 3], [3, 5]]])
547 >>> inv(a)
548 array([[[-2. , 1. ],
549 [ 1.5 , -0.5 ]],
550 [[-1.25, 0.75],
551 [ 0.75, -0.25]]])
553 """
554 a, wrap = _makearray(a)
555 _assert_stacked_2d(a)
556 _assert_stacked_square(a)
557 t, result_t = _commonType(a)
559 signature = 'D->D' if isComplexType(t) else 'd->d'
560 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
561 ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
562 return wrap(ainv.astype(result_t, copy=False))
565def _matrix_power_dispatcher(a, n):
566 return (a,)
569@array_function_dispatch(_matrix_power_dispatcher)
570def matrix_power(a, n):
571 """
572 Raise a square matrix to the (integer) power `n`.
574 For positive integers `n`, the power is computed by repeated matrix
575 squarings and matrix multiplications. If ``n == 0``, the identity matrix
576 of the same shape as M is returned. If ``n < 0``, the inverse
577 is computed and then raised to the ``abs(n)``.
579 .. note:: Stacks of object matrices are not currently supported.
581 Parameters
582 ----------
583 a : (..., M, M) array_like
584 Matrix to be "powered".
585 n : int
586 The exponent can be any integer or long integer, positive,
587 negative, or zero.
589 Returns
590 -------
591 a**n : (..., M, M) ndarray or matrix object
592 The return value is the same shape and type as `M`;
593 if the exponent is positive or zero then the type of the
594 elements is the same as those of `M`. If the exponent is
595 negative the elements are floating-point.
597 Raises
598 ------
599 LinAlgError
600 For matrices that are not square or that (for negative powers) cannot
601 be inverted numerically.
603 Examples
604 --------
605 >>> from numpy.linalg import matrix_power
606 >>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit
607 >>> matrix_power(i, 3) # should = -i
608 array([[ 0, -1],
609 [ 1, 0]])
610 >>> matrix_power(i, 0)
611 array([[1, 0],
612 [0, 1]])
613 >>> matrix_power(i, -3) # should = 1/(-i) = i, but w/ f.p. elements
614 array([[ 0., 1.],
615 [-1., 0.]])
617 Somewhat more sophisticated example
619 >>> q = np.zeros((4, 4))
620 >>> q[0:2, 0:2] = -i
621 >>> q[2:4, 2:4] = i
622 >>> q # one of the three quaternion units not equal to 1
623 array([[ 0., -1., 0., 0.],
624 [ 1., 0., 0., 0.],
625 [ 0., 0., 0., 1.],
626 [ 0., 0., -1., 0.]])
627 >>> matrix_power(q, 2) # = -np.eye(4)
628 array([[-1., 0., 0., 0.],
629 [ 0., -1., 0., 0.],
630 [ 0., 0., -1., 0.],
631 [ 0., 0., 0., -1.]])
633 """
634 a = asanyarray(a)
635 _assert_stacked_2d(a)
636 _assert_stacked_square(a)
638 try:
639 n = operator.index(n)
640 except TypeError as e:
641 raise TypeError("exponent must be an integer") from e
643 # Fall back on dot for object arrays. Object arrays are not supported by
644 # the current implementation of matmul using einsum
645 if a.dtype != object:
646 fmatmul = matmul
647 elif a.ndim == 2:
648 fmatmul = dot
649 else:
650 raise NotImplementedError(
651 "matrix_power not supported for stacks of object arrays")
653 if n == 0:
654 a = empty_like(a)
655 a[...] = eye(a.shape[-2], dtype=a.dtype)
656 return a
658 elif n < 0:
659 a = inv(a)
660 n = abs(n)
662 # short-cuts.
663 if n == 1:
664 return a
666 elif n == 2:
667 return fmatmul(a, a)
669 elif n == 3:
670 return fmatmul(fmatmul(a, a), a)
672 # Use binary decomposition to reduce the number of matrix multiplications.
673 # Here, we iterate over the bits of n, from LSB to MSB, raise `a` to
674 # increasing powers of 2, and multiply into the result as needed.
675 z = result = None
676 while n > 0:
677 z = a if z is None else fmatmul(z, z)
678 n, bit = divmod(n, 2)
679 if bit:
680 result = z if result is None else fmatmul(result, z)
682 return result
685# Cholesky decomposition
688@array_function_dispatch(_unary_dispatcher)
689def cholesky(a):
690 """
691 Cholesky decomposition.
693 Return the Cholesky decomposition, `L * L.H`, of the square matrix `a`,
694 where `L` is lower-triangular and .H is the conjugate transpose operator
695 (which is the ordinary transpose if `a` is real-valued). `a` must be
696 Hermitian (symmetric if real-valued) and positive-definite. No
697 checking is performed to verify whether `a` is Hermitian or not.
698 In addition, only the lower-triangular and diagonal elements of `a`
699 are used. Only `L` is actually returned.
701 Parameters
702 ----------
703 a : (..., M, M) array_like
704 Hermitian (symmetric if all elements are real), positive-definite
705 input matrix.
707 Returns
708 -------
709 L : (..., M, M) array_like
710 Lower-triangular Cholesky factor of `a`. Returns a matrix object if
711 `a` is a matrix object.
713 Raises
714 ------
715 LinAlgError
716 If the decomposition fails, for example, if `a` is not
717 positive-definite.
719 See Also
720 --------
721 scipy.linalg.cholesky : Similar function in SciPy.
722 scipy.linalg.cholesky_banded : Cholesky decompose a banded Hermitian
723 positive-definite matrix.
724 scipy.linalg.cho_factor : Cholesky decomposition of a matrix, to use in
725 `scipy.linalg.cho_solve`.
727 Notes
728 -----
730 .. versionadded:: 1.8.0
732 Broadcasting rules apply, see the `numpy.linalg` documentation for
733 details.
735 The Cholesky decomposition is often used as a fast way of solving
737 .. math:: A \\mathbf{x} = \\mathbf{b}
739 (when `A` is both Hermitian/symmetric and positive-definite).
741 First, we solve for :math:`\\mathbf{y}` in
743 .. math:: L \\mathbf{y} = \\mathbf{b},
745 and then for :math:`\\mathbf{x}` in
747 .. math:: L.H \\mathbf{x} = \\mathbf{y}.
749 Examples
750 --------
751 >>> A = np.array([[1,-2j],[2j,5]])
752 >>> A
753 array([[ 1.+0.j, -0.-2.j],
754 [ 0.+2.j, 5.+0.j]])
755 >>> L = np.linalg.cholesky(A)
756 >>> L
757 array([[1.+0.j, 0.+0.j],
758 [0.+2.j, 1.+0.j]])
759 >>> np.dot(L, L.T.conj()) # verify that L * L.H = A
760 array([[1.+0.j, 0.-2.j],
761 [0.+2.j, 5.+0.j]])
762 >>> A = [[1,-2j],[2j,5]] # what happens if A is only array_like?
763 >>> np.linalg.cholesky(A) # an ndarray object is returned
764 array([[1.+0.j, 0.+0.j],
765 [0.+2.j, 1.+0.j]])
766 >>> # But a matrix object is returned if A is a matrix object
767 >>> np.linalg.cholesky(np.matrix(A))
768 matrix([[ 1.+0.j, 0.+0.j],
769 [ 0.+2.j, 1.+0.j]])
771 """
772 extobj = get_linalg_error_extobj(_raise_linalgerror_nonposdef)
773 gufunc = _umath_linalg.cholesky_lo
774 a, wrap = _makearray(a)
775 _assert_stacked_2d(a)
776 _assert_stacked_square(a)
777 t, result_t = _commonType(a)
778 signature = 'D->D' if isComplexType(t) else 'd->d'
779 r = gufunc(a, signature=signature, extobj=extobj)
780 return wrap(r.astype(result_t, copy=False))
783# QR decomposition
785def _qr_dispatcher(a, mode=None):
786 return (a,)
789@array_function_dispatch(_qr_dispatcher)
790def qr(a, mode='reduced'):
791 """
792 Compute the qr factorization of a matrix.
794 Factor the matrix `a` as *qr*, where `q` is orthonormal and `r` is
795 upper-triangular.
797 Parameters
798 ----------
799 a : array_like, shape (..., M, N)
800 An array-like object with the dimensionality of at least 2.
801 mode : {'reduced', 'complete', 'r', 'raw'}, optional
802 If K = min(M, N), then
804 * 'reduced' : returns Q, R with dimensions (..., M, K), (..., K, N) (default)
805 * 'complete' : returns Q, R with dimensions (..., M, M), (..., M, N)
806 * 'r' : returns R only with dimensions (..., K, N)
807 * 'raw' : returns h, tau with dimensions (..., N, M), (..., K,)
809 The options 'reduced', 'complete, and 'raw' are new in numpy 1.8,
810 see the notes for more information. The default is 'reduced', and to
811 maintain backward compatibility with earlier versions of numpy both
812 it and the old default 'full' can be omitted. Note that array h
813 returned in 'raw' mode is transposed for calling Fortran. The
814 'economic' mode is deprecated. The modes 'full' and 'economic' may
815 be passed using only the first letter for backwards compatibility,
816 but all others must be spelled out. See the Notes for more
817 explanation.
820 Returns
821 -------
822 When mode is 'reduced' or 'complete', the result will be a namedtuple with
823 the attributes `Q` and `R`.
825 Q : ndarray of float or complex, optional
826 A matrix with orthonormal columns. When mode = 'complete' the
827 result is an orthogonal/unitary matrix depending on whether or not
828 a is real/complex. The determinant may be either +/- 1 in that
829 case. In case the number of dimensions in the input array is
830 greater than 2 then a stack of the matrices with above properties
831 is returned.
832 R : ndarray of float or complex, optional
833 The upper-triangular matrix or a stack of upper-triangular
834 matrices if the number of dimensions in the input array is greater
835 than 2.
836 (h, tau) : ndarrays of np.double or np.cdouble, optional
837 The array h contains the Householder reflectors that generate q
838 along with r. The tau array contains scaling factors for the
839 reflectors. In the deprecated 'economic' mode only h is returned.
841 Raises
842 ------
843 LinAlgError
844 If factoring fails.
846 See Also
847 --------
848 scipy.linalg.qr : Similar function in SciPy.
849 scipy.linalg.rq : Compute RQ decomposition of a matrix.
851 Notes
852 -----
853 This is an interface to the LAPACK routines ``dgeqrf``, ``zgeqrf``,
854 ``dorgqr``, and ``zungqr``.
856 For more information on the qr factorization, see for example:
857 https://en.wikipedia.org/wiki/QR_factorization
859 Subclasses of `ndarray` are preserved except for the 'raw' mode. So if
860 `a` is of type `matrix`, all the return values will be matrices too.
862 New 'reduced', 'complete', and 'raw' options for mode were added in
863 NumPy 1.8.0 and the old option 'full' was made an alias of 'reduced'. In
864 addition the options 'full' and 'economic' were deprecated. Because
865 'full' was the previous default and 'reduced' is the new default,
866 backward compatibility can be maintained by letting `mode` default.
867 The 'raw' option was added so that LAPACK routines that can multiply
868 arrays by q using the Householder reflectors can be used. Note that in
869 this case the returned arrays are of type np.double or np.cdouble and
870 the h array is transposed to be FORTRAN compatible. No routines using
871 the 'raw' return are currently exposed by numpy, but some are available
872 in lapack_lite and just await the necessary work.
874 Examples
875 --------
876 >>> a = np.random.randn(9, 6)
877 >>> Q, R = np.linalg.qr(a)
878 >>> np.allclose(a, np.dot(Q, R)) # a does equal QR
879 True
880 >>> R2 = np.linalg.qr(a, mode='r')
881 >>> np.allclose(R, R2) # mode='r' returns the same R as mode='full'
882 True
883 >>> a = np.random.normal(size=(3, 2, 2)) # Stack of 2 x 2 matrices as input
884 >>> Q, R = np.linalg.qr(a)
885 >>> Q.shape
886 (3, 2, 2)
887 >>> R.shape
888 (3, 2, 2)
889 >>> np.allclose(a, np.matmul(Q, R))
890 True
892 Example illustrating a common use of `qr`: solving of least squares
893 problems
895 What are the least-squares-best `m` and `y0` in ``y = y0 + mx`` for
896 the following data: {(0,1), (1,0), (1,2), (2,1)}. (Graph the points
897 and you'll see that it should be y0 = 0, m = 1.) The answer is provided
898 by solving the over-determined matrix equation ``Ax = b``, where::
900 A = array([[0, 1], [1, 1], [1, 1], [2, 1]])
901 x = array([[y0], [m]])
902 b = array([[1], [0], [2], [1]])
904 If A = QR such that Q is orthonormal (which is always possible via
905 Gram-Schmidt), then ``x = inv(R) * (Q.T) * b``. (In numpy practice,
906 however, we simply use `lstsq`.)
908 >>> A = np.array([[0, 1], [1, 1], [1, 1], [2, 1]])
909 >>> A
910 array([[0, 1],
911 [1, 1],
912 [1, 1],
913 [2, 1]])
914 >>> b = np.array([1, 2, 2, 3])
915 >>> Q, R = np.linalg.qr(A)
916 >>> p = np.dot(Q.T, b)
917 >>> np.dot(np.linalg.inv(R), p)
918 array([ 1., 1.])
920 """
921 if mode not in ('reduced', 'complete', 'r', 'raw'):
922 if mode in ('f', 'full'):
923 # 2013-04-01, 1.8
924 msg = "".join((
925 "The 'full' option is deprecated in favor of 'reduced'.\n",
926 "For backward compatibility let mode default."))
927 warnings.warn(msg, DeprecationWarning, stacklevel=2)
928 mode = 'reduced'
929 elif mode in ('e', 'economic'):
930 # 2013-04-01, 1.8
931 msg = "The 'economic' option is deprecated."
932 warnings.warn(msg, DeprecationWarning, stacklevel=2)
933 mode = 'economic'
934 else:
935 raise ValueError(f"Unrecognized mode '{mode}'")
937 a, wrap = _makearray(a)
938 _assert_stacked_2d(a)
939 m, n = a.shape[-2:]
940 t, result_t = _commonType(a)
941 a = a.astype(t, copy=True)
942 a = _to_native_byte_order(a)
943 mn = min(m, n)
945 if m <= n:
946 gufunc = _umath_linalg.qr_r_raw_m
947 else:
948 gufunc = _umath_linalg.qr_r_raw_n
950 signature = 'D->D' if isComplexType(t) else 'd->d'
951 extobj = get_linalg_error_extobj(_raise_linalgerror_qr)
952 tau = gufunc(a, signature=signature, extobj=extobj)
954 # handle modes that don't return q
955 if mode == 'r':
956 r = triu(a[..., :mn, :])
957 r = r.astype(result_t, copy=False)
958 return wrap(r)
960 if mode == 'raw':
961 q = transpose(a)
962 q = q.astype(result_t, copy=False)
963 tau = tau.astype(result_t, copy=False)
964 return wrap(q), tau
966 if mode == 'economic':
967 a = a.astype(result_t, copy=False)
968 return wrap(a)
970 # mc is the number of columns in the resulting q
971 # matrix. If the mode is complete then it is
972 # same as number of rows, and if the mode is reduced,
973 # then it is the minimum of number of rows and columns.
974 if mode == 'complete' and m > n:
975 mc = m
976 gufunc = _umath_linalg.qr_complete
977 else:
978 mc = mn
979 gufunc = _umath_linalg.qr_reduced
981 signature = 'DD->D' if isComplexType(t) else 'dd->d'
982 extobj = get_linalg_error_extobj(_raise_linalgerror_qr)
983 q = gufunc(a, tau, signature=signature, extobj=extobj)
984 r = triu(a[..., :mc, :])
986 q = q.astype(result_t, copy=False)
987 r = r.astype(result_t, copy=False)
989 return QRResult(wrap(q), wrap(r))
991# Eigenvalues
994@array_function_dispatch(_unary_dispatcher)
995def eigvals(a):
996 """
997 Compute the eigenvalues of a general matrix.
999 Main difference between `eigvals` and `eig`: the eigenvectors aren't
1000 returned.
1002 Parameters
1003 ----------
1004 a : (..., M, M) array_like
1005 A complex- or real-valued matrix whose eigenvalues will be computed.
1007 Returns
1008 -------
1009 w : (..., M,) ndarray
1010 The eigenvalues, each repeated according to its multiplicity.
1011 They are not necessarily ordered, nor are they necessarily
1012 real for real matrices.
1014 Raises
1015 ------
1016 LinAlgError
1017 If the eigenvalue computation does not converge.
1019 See Also
1020 --------
1021 eig : eigenvalues and right eigenvectors of general arrays
1022 eigvalsh : eigenvalues of real symmetric or complex Hermitian
1023 (conjugate symmetric) arrays.
1024 eigh : eigenvalues and eigenvectors of real symmetric or complex
1025 Hermitian (conjugate symmetric) arrays.
1026 scipy.linalg.eigvals : Similar function in SciPy.
1028 Notes
1029 -----
1031 .. versionadded:: 1.8.0
1033 Broadcasting rules apply, see the `numpy.linalg` documentation for
1034 details.
1036 This is implemented using the ``_geev`` LAPACK routines which compute
1037 the eigenvalues and eigenvectors of general square arrays.
1039 Examples
1040 --------
1041 Illustration, using the fact that the eigenvalues of a diagonal matrix
1042 are its diagonal elements, that multiplying a matrix on the left
1043 by an orthogonal matrix, `Q`, and on the right by `Q.T` (the transpose
1044 of `Q`), preserves the eigenvalues of the "middle" matrix. In other words,
1045 if `Q` is orthogonal, then ``Q * A * Q.T`` has the same eigenvalues as
1046 ``A``:
1048 >>> from numpy import linalg as LA
1049 >>> x = np.random.random()
1050 >>> Q = np.array([[np.cos(x), -np.sin(x)], [np.sin(x), np.cos(x)]])
1051 >>> LA.norm(Q[0, :]), LA.norm(Q[1, :]), np.dot(Q[0, :],Q[1, :])
1052 (1.0, 1.0, 0.0)
1054 Now multiply a diagonal matrix by ``Q`` on one side and by ``Q.T`` on the other:
1056 >>> D = np.diag((-1,1))
1057 >>> LA.eigvals(D)
1058 array([-1., 1.])
1059 >>> A = np.dot(Q, D)
1060 >>> A = np.dot(A, Q.T)
1061 >>> LA.eigvals(A)
1062 array([ 1., -1.]) # random
1064 """
1065 a, wrap = _makearray(a)
1066 _assert_stacked_2d(a)
1067 _assert_stacked_square(a)
1068 _assert_finite(a)
1069 t, result_t = _commonType(a)
1071 extobj = get_linalg_error_extobj(
1072 _raise_linalgerror_eigenvalues_nonconvergence)
1073 signature = 'D->D' if isComplexType(t) else 'd->D'
1074 w = _umath_linalg.eigvals(a, signature=signature, extobj=extobj)
1076 if not isComplexType(t):
1077 if all(w.imag == 0):
1078 w = w.real
1079 result_t = _realType(result_t)
1080 else:
1081 result_t = _complexType(result_t)
1083 return w.astype(result_t, copy=False)
1086def _eigvalsh_dispatcher(a, UPLO=None):
1087 return (a,)
1090@array_function_dispatch(_eigvalsh_dispatcher)
1091def eigvalsh(a, UPLO='L'):
1092 """
1093 Compute the eigenvalues of a complex Hermitian or real symmetric matrix.
1095 Main difference from eigh: the eigenvectors are not computed.
1097 Parameters
1098 ----------
1099 a : (..., M, M) array_like
1100 A complex- or real-valued matrix whose eigenvalues are to be
1101 computed.
1102 UPLO : {'L', 'U'}, optional
1103 Specifies whether the calculation is done with the lower triangular
1104 part of `a` ('L', default) or the upper triangular part ('U').
1105 Irrespective of this value only the real parts of the diagonal will
1106 be considered in the computation to preserve the notion of a Hermitian
1107 matrix. It therefore follows that the imaginary part of the diagonal
1108 will always be treated as zero.
1110 Returns
1111 -------
1112 w : (..., M,) ndarray
1113 The eigenvalues in ascending order, each repeated according to
1114 its multiplicity.
1116 Raises
1117 ------
1118 LinAlgError
1119 If the eigenvalue computation does not converge.
1121 See Also
1122 --------
1123 eigh : eigenvalues and eigenvectors of real symmetric or complex Hermitian
1124 (conjugate symmetric) arrays.
1125 eigvals : eigenvalues of general real or complex arrays.
1126 eig : eigenvalues and right eigenvectors of general real or complex
1127 arrays.
1128 scipy.linalg.eigvalsh : Similar function in SciPy.
1130 Notes
1131 -----
1133 .. versionadded:: 1.8.0
1135 Broadcasting rules apply, see the `numpy.linalg` documentation for
1136 details.
1138 The eigenvalues are computed using LAPACK routines ``_syevd``, ``_heevd``.
1140 Examples
1141 --------
1142 >>> from numpy import linalg as LA
1143 >>> a = np.array([[1, -2j], [2j, 5]])
1144 >>> LA.eigvalsh(a)
1145 array([ 0.17157288, 5.82842712]) # may vary
1147 >>> # demonstrate the treatment of the imaginary part of the diagonal
1148 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1149 >>> a
1150 array([[5.+2.j, 9.-2.j],
1151 [0.+2.j, 2.-1.j]])
1152 >>> # with UPLO='L' this is numerically equivalent to using LA.eigvals()
1153 >>> # with:
1154 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1155 >>> b
1156 array([[5.+0.j, 0.-2.j],
1157 [0.+2.j, 2.+0.j]])
1158 >>> wa = LA.eigvalsh(a)
1159 >>> wb = LA.eigvals(b)
1160 >>> wa; wb
1161 array([1., 6.])
1162 array([6.+0.j, 1.+0.j])
1164 """
1165 UPLO = UPLO.upper()
1166 if UPLO not in ('L', 'U'):
1167 raise ValueError("UPLO argument must be 'L' or 'U'")
1169 extobj = get_linalg_error_extobj(
1170 _raise_linalgerror_eigenvalues_nonconvergence)
1171 if UPLO == 'L':
1172 gufunc = _umath_linalg.eigvalsh_lo
1173 else:
1174 gufunc = _umath_linalg.eigvalsh_up
1176 a, wrap = _makearray(a)
1177 _assert_stacked_2d(a)
1178 _assert_stacked_square(a)
1179 t, result_t = _commonType(a)
1180 signature = 'D->d' if isComplexType(t) else 'd->d'
1181 w = gufunc(a, signature=signature, extobj=extobj)
1182 return w.astype(_realType(result_t), copy=False)
1184def _convertarray(a):
1185 t, result_t = _commonType(a)
1186 a = a.astype(t).T.copy()
1187 return a, t, result_t
1190# Eigenvectors
1193@array_function_dispatch(_unary_dispatcher)
1194def eig(a):
1195 """
1196 Compute the eigenvalues and right eigenvectors of a square array.
1198 Parameters
1199 ----------
1200 a : (..., M, M) array
1201 Matrices for which the eigenvalues and right eigenvectors will
1202 be computed
1204 Returns
1205 -------
1206 A namedtuple with the following attributes:
1208 eigenvalues : (..., M) array
1209 The eigenvalues, each repeated according to its multiplicity.
1210 The eigenvalues are not necessarily ordered. The resulting
1211 array will be of complex type, unless the imaginary part is
1212 zero in which case it will be cast to a real type. When `a`
1213 is real the resulting eigenvalues will be real (0 imaginary
1214 part) or occur in conjugate pairs
1216 eigenvectors : (..., M, M) array
1217 The normalized (unit "length") eigenvectors, such that the
1218 column ``eigenvectors[:,i]`` is the eigenvector corresponding to the
1219 eigenvalue ``eigenvalues[i]``.
1221 Raises
1222 ------
1223 LinAlgError
1224 If the eigenvalue computation does not converge.
1226 See Also
1227 --------
1228 eigvals : eigenvalues of a non-symmetric array.
1229 eigh : eigenvalues and eigenvectors of a real symmetric or complex
1230 Hermitian (conjugate symmetric) array.
1231 eigvalsh : eigenvalues of a real symmetric or complex Hermitian
1232 (conjugate symmetric) array.
1233 scipy.linalg.eig : Similar function in SciPy that also solves the
1234 generalized eigenvalue problem.
1235 scipy.linalg.schur : Best choice for unitary and other non-Hermitian
1236 normal matrices.
1238 Notes
1239 -----
1241 .. versionadded:: 1.8.0
1243 Broadcasting rules apply, see the `numpy.linalg` documentation for
1244 details.
1246 This is implemented using the ``_geev`` LAPACK routines which compute
1247 the eigenvalues and eigenvectors of general square arrays.
1249 The number `w` is an eigenvalue of `a` if there exists a vector `v` such
1250 that ``a @ v = w * v``. Thus, the arrays `a`, `eigenvalues`, and
1251 `eigenvectors` satisfy the equations ``a @ eigenvectors[:,i] =
1252 eigenvalues[i] * eigenvalues[:,i]`` for :math:`i \\in \\{0,...,M-1\\}`.
1254 The array `eigenvectors` may not be of maximum rank, that is, some of the
1255 columns may be linearly dependent, although round-off error may obscure
1256 that fact. If the eigenvalues are all different, then theoretically the
1257 eigenvectors are linearly independent and `a` can be diagonalized by a
1258 similarity transformation using `eigenvectors`, i.e, ``inv(eigenvectors) @
1259 a @ eigenvectors`` is diagonal.
1261 For non-Hermitian normal matrices the SciPy function `scipy.linalg.schur`
1262 is preferred because the matrix `eigenvectors` is guaranteed to be
1263 unitary, which is not the case when using `eig`. The Schur factorization
1264 produces an upper triangular matrix rather than a diagonal matrix, but for
1265 normal matrices only the diagonal of the upper triangular matrix is
1266 needed, the rest is roundoff error.
1268 Finally, it is emphasized that `eigenvectors` consists of the *right* (as
1269 in right-hand side) eigenvectors of `a`. A vector `y` satisfying ``y.T @ a
1270 = z * y.T`` for some number `z` is called a *left* eigenvector of `a`,
1271 and, in general, the left and right eigenvectors of a matrix are not
1272 necessarily the (perhaps conjugate) transposes of each other.
1274 References
1275 ----------
1276 G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, FL,
1277 Academic Press, Inc., 1980, Various pp.
1279 Examples
1280 --------
1281 >>> from numpy import linalg as LA
1283 (Almost) trivial example with real eigenvalues and eigenvectors.
1285 >>> eigenvalues, eigenvectors = LA.eig(np.diag((1, 2, 3)))
1286 >>> eigenvalues
1287 array([1., 2., 3.])
1288 >>> eigenvectors
1289 array([[1., 0., 0.],
1290 [0., 1., 0.],
1291 [0., 0., 1.]])
1293 Real matrix possessing complex eigenvalues and eigenvectors; note that the
1294 eigenvalues are complex conjugates of each other.
1296 >>> eigenvalues, eigenvectors = LA.eig(np.array([[1, -1], [1, 1]]))
1297 >>> eigenvalues
1298 array([1.+1.j, 1.-1.j])
1299 >>> eigenvectors
1300 array([[0.70710678+0.j , 0.70710678-0.j ],
1301 [0. -0.70710678j, 0. +0.70710678j]])
1303 Complex-valued matrix with real eigenvalues (but complex-valued eigenvectors);
1304 note that ``a.conj().T == a``, i.e., `a` is Hermitian.
1306 >>> a = np.array([[1, 1j], [-1j, 1]])
1307 >>> eigenvalues, eigenvectors = LA.eig(a)
1308 >>> eigenvalues
1309 array([2.+0.j, 0.+0.j])
1310 >>> eigenvectors
1311 array([[ 0. +0.70710678j, 0.70710678+0.j ], # may vary
1312 [ 0.70710678+0.j , -0. +0.70710678j]])
1314 Be careful about round-off error!
1316 >>> a = np.array([[1 + 1e-9, 0], [0, 1 - 1e-9]])
1317 >>> # Theor. eigenvalues are 1 +/- 1e-9
1318 >>> eigenvalues, eigenvectors = LA.eig(a)
1319 >>> eigenvalues
1320 array([1., 1.])
1321 >>> eigenvectors
1322 array([[1., 0.],
1323 [0., 1.]])
1325 """
1326 a, wrap = _makearray(a)
1327 _assert_stacked_2d(a)
1328 _assert_stacked_square(a)
1329 _assert_finite(a)
1330 t, result_t = _commonType(a)
1332 extobj = get_linalg_error_extobj(
1333 _raise_linalgerror_eigenvalues_nonconvergence)
1334 signature = 'D->DD' if isComplexType(t) else 'd->DD'
1335 w, vt = _umath_linalg.eig(a, signature=signature, extobj=extobj)
1337 if not isComplexType(t) and all(w.imag == 0.0):
1338 w = w.real
1339 vt = vt.real
1340 result_t = _realType(result_t)
1341 else:
1342 result_t = _complexType(result_t)
1344 vt = vt.astype(result_t, copy=False)
1345 return EigResult(w.astype(result_t, copy=False), wrap(vt))
1348@array_function_dispatch(_eigvalsh_dispatcher)
1349def eigh(a, UPLO='L'):
1350 """
1351 Return the eigenvalues and eigenvectors of a complex Hermitian
1352 (conjugate symmetric) or a real symmetric matrix.
1354 Returns two objects, a 1-D array containing the eigenvalues of `a`, and
1355 a 2-D square array or matrix (depending on the input type) of the
1356 corresponding eigenvectors (in columns).
1358 Parameters
1359 ----------
1360 a : (..., M, M) array
1361 Hermitian or real symmetric matrices whose eigenvalues and
1362 eigenvectors are to be computed.
1363 UPLO : {'L', 'U'}, optional
1364 Specifies whether the calculation is done with the lower triangular
1365 part of `a` ('L', default) or the upper triangular part ('U').
1366 Irrespective of this value only the real parts of the diagonal will
1367 be considered in the computation to preserve the notion of a Hermitian
1368 matrix. It therefore follows that the imaginary part of the diagonal
1369 will always be treated as zero.
1371 Returns
1372 -------
1373 A namedtuple with the following attributes:
1375 eigenvalues : (..., M) ndarray
1376 The eigenvalues in ascending order, each repeated according to
1377 its multiplicity.
1378 eigenvectors : {(..., M, M) ndarray, (..., M, M) matrix}
1379 The column ``eigenvectors[:, i]`` is the normalized eigenvector
1380 corresponding to the eigenvalue ``eigenvalues[i]``. Will return a
1381 matrix object if `a` is a matrix object.
1383 Raises
1384 ------
1385 LinAlgError
1386 If the eigenvalue computation does not converge.
1388 See Also
1389 --------
1390 eigvalsh : eigenvalues of real symmetric or complex Hermitian
1391 (conjugate symmetric) arrays.
1392 eig : eigenvalues and right eigenvectors for non-symmetric arrays.
1393 eigvals : eigenvalues of non-symmetric arrays.
1394 scipy.linalg.eigh : Similar function in SciPy (but also solves the
1395 generalized eigenvalue problem).
1397 Notes
1398 -----
1400 .. versionadded:: 1.8.0
1402 Broadcasting rules apply, see the `numpy.linalg` documentation for
1403 details.
1405 The eigenvalues/eigenvectors are computed using LAPACK routines ``_syevd``,
1406 ``_heevd``.
1408 The eigenvalues of real symmetric or complex Hermitian matrices are always
1409 real. [1]_ The array `eigenvalues` of (column) eigenvectors is unitary and
1410 `a`, `eigenvalues`, and `eigenvectors` satisfy the equations ``dot(a,
1411 eigenvectors[:, i]) = eigenvalues[i] * eigenvectors[:, i]``.
1413 References
1414 ----------
1415 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
1416 FL, Academic Press, Inc., 1980, pg. 222.
1418 Examples
1419 --------
1420 >>> from numpy import linalg as LA
1421 >>> a = np.array([[1, -2j], [2j, 5]])
1422 >>> a
1423 array([[ 1.+0.j, -0.-2.j],
1424 [ 0.+2.j, 5.+0.j]])
1425 >>> eigenvalues, eigenvectors = LA.eigh(a)
1426 >>> eigenvalues
1427 array([0.17157288, 5.82842712])
1428 >>> eigenvectors
1429 array([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1430 [ 0. +0.38268343j, 0. -0.92387953j]])
1432 >>> np.dot(a, eigenvectors[:, 0]) - eigenvalues[0] * eigenvectors[:, 0] # verify 1st eigenval/vec pair
1433 array([5.55111512e-17+0.0000000e+00j, 0.00000000e+00+1.2490009e-16j])
1434 >>> np.dot(a, eigenvectors[:, 1]) - eigenvalues[1] * eigenvectors[:, 1] # verify 2nd eigenval/vec pair
1435 array([0.+0.j, 0.+0.j])
1437 >>> A = np.matrix(a) # what happens if input is a matrix object
1438 >>> A
1439 matrix([[ 1.+0.j, -0.-2.j],
1440 [ 0.+2.j, 5.+0.j]])
1441 >>> eigenvalues, eigenvectors = LA.eigh(A)
1442 >>> eigenvalues
1443 array([0.17157288, 5.82842712])
1444 >>> eigenvectors
1445 matrix([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1446 [ 0. +0.38268343j, 0. -0.92387953j]])
1448 >>> # demonstrate the treatment of the imaginary part of the diagonal
1449 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1450 >>> a
1451 array([[5.+2.j, 9.-2.j],
1452 [0.+2.j, 2.-1.j]])
1453 >>> # with UPLO='L' this is numerically equivalent to using LA.eig() with:
1454 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1455 >>> b
1456 array([[5.+0.j, 0.-2.j],
1457 [0.+2.j, 2.+0.j]])
1458 >>> wa, va = LA.eigh(a)
1459 >>> wb, vb = LA.eig(b)
1460 >>> wa; wb
1461 array([1., 6.])
1462 array([6.+0.j, 1.+0.j])
1463 >>> va; vb
1464 array([[-0.4472136 +0.j , -0.89442719+0.j ], # may vary
1465 [ 0. +0.89442719j, 0. -0.4472136j ]])
1466 array([[ 0.89442719+0.j , -0. +0.4472136j],
1467 [-0. +0.4472136j, 0.89442719+0.j ]])
1469 """
1470 UPLO = UPLO.upper()
1471 if UPLO not in ('L', 'U'):
1472 raise ValueError("UPLO argument must be 'L' or 'U'")
1474 a, wrap = _makearray(a)
1475 _assert_stacked_2d(a)
1476 _assert_stacked_square(a)
1477 t, result_t = _commonType(a)
1479 extobj = get_linalg_error_extobj(
1480 _raise_linalgerror_eigenvalues_nonconvergence)
1481 if UPLO == 'L':
1482 gufunc = _umath_linalg.eigh_lo
1483 else:
1484 gufunc = _umath_linalg.eigh_up
1486 signature = 'D->dD' if isComplexType(t) else 'd->dd'
1487 w, vt = gufunc(a, signature=signature, extobj=extobj)
1488 w = w.astype(_realType(result_t), copy=False)
1489 vt = vt.astype(result_t, copy=False)
1490 return EighResult(w, wrap(vt))
1493# Singular value decomposition
1495def _svd_dispatcher(a, full_matrices=None, compute_uv=None, hermitian=None):
1496 return (a,)
1499@array_function_dispatch(_svd_dispatcher)
1500def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
1501 """
1502 Singular Value Decomposition.
1504 When `a` is a 2D array, and ``full_matrices=False``, then it is
1505 factorized as ``u @ np.diag(s) @ vh = (u * s) @ vh``, where
1506 `u` and the Hermitian transpose of `vh` are 2D arrays with
1507 orthonormal columns and `s` is a 1D array of `a`'s singular
1508 values. When `a` is higher-dimensional, SVD is applied in
1509 stacked mode as explained below.
1511 Parameters
1512 ----------
1513 a : (..., M, N) array_like
1514 A real or complex array with ``a.ndim >= 2``.
1515 full_matrices : bool, optional
1516 If True (default), `u` and `vh` have the shapes ``(..., M, M)`` and
1517 ``(..., N, N)``, respectively. Otherwise, the shapes are
1518 ``(..., M, K)`` and ``(..., K, N)``, respectively, where
1519 ``K = min(M, N)``.
1520 compute_uv : bool, optional
1521 Whether or not to compute `u` and `vh` in addition to `s`. True
1522 by default.
1523 hermitian : bool, optional
1524 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
1525 enabling a more efficient method for finding singular values.
1526 Defaults to False.
1528 .. versionadded:: 1.17.0
1530 Returns
1531 -------
1532 When `compute_uv` is True, the result is a namedtuple with the following
1533 attribute names:
1535 U : { (..., M, M), (..., M, K) } array
1536 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1537 size as those of the input `a`. The size of the last two dimensions
1538 depends on the value of `full_matrices`. Only returned when
1539 `compute_uv` is True.
1540 S : (..., K) array
1541 Vector(s) with the singular values, within each vector sorted in
1542 descending order. The first ``a.ndim - 2`` dimensions have the same
1543 size as those of the input `a`.
1544 Vh : { (..., N, N), (..., K, N) } array
1545 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1546 size as those of the input `a`. The size of the last two dimensions
1547 depends on the value of `full_matrices`. Only returned when
1548 `compute_uv` is True.
1550 Raises
1551 ------
1552 LinAlgError
1553 If SVD computation does not converge.
1555 See Also
1556 --------
1557 scipy.linalg.svd : Similar function in SciPy.
1558 scipy.linalg.svdvals : Compute singular values of a matrix.
1560 Notes
1561 -----
1563 .. versionchanged:: 1.8.0
1564 Broadcasting rules apply, see the `numpy.linalg` documentation for
1565 details.
1567 The decomposition is performed using LAPACK routine ``_gesdd``.
1569 SVD is usually described for the factorization of a 2D matrix :math:`A`.
1570 The higher-dimensional case will be discussed below. In the 2D case, SVD is
1571 written as :math:`A = U S V^H`, where :math:`A = a`, :math:`U= u`,
1572 :math:`S= \\mathtt{np.diag}(s)` and :math:`V^H = vh`. The 1D array `s`
1573 contains the singular values of `a` and `u` and `vh` are unitary. The rows
1574 of `vh` are the eigenvectors of :math:`A^H A` and the columns of `u` are
1575 the eigenvectors of :math:`A A^H`. In both cases the corresponding
1576 (possibly non-zero) eigenvalues are given by ``s**2``.
1578 If `a` has more than two dimensions, then broadcasting rules apply, as
1579 explained in :ref:`routines.linalg-broadcasting`. This means that SVD is
1580 working in "stacked" mode: it iterates over all indices of the first
1581 ``a.ndim - 2`` dimensions and for each combination SVD is applied to the
1582 last two indices. The matrix `a` can be reconstructed from the
1583 decomposition with either ``(u * s[..., None, :]) @ vh`` or
1584 ``u @ (s[..., None] * vh)``. (The ``@`` operator can be replaced by the
1585 function ``np.matmul`` for python versions below 3.5.)
1587 If `a` is a ``matrix`` object (as opposed to an ``ndarray``), then so are
1588 all the return values.
1590 Examples
1591 --------
1592 >>> a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6)
1593 >>> b = np.random.randn(2, 7, 8, 3) + 1j*np.random.randn(2, 7, 8, 3)
1595 Reconstruction based on full SVD, 2D case:
1597 >>> U, S, Vh = np.linalg.svd(a, full_matrices=True)
1598 >>> U.shape, S.shape, Vh.shape
1599 ((9, 9), (6,), (6, 6))
1600 >>> np.allclose(a, np.dot(U[:, :6] * S, Vh))
1601 True
1602 >>> smat = np.zeros((9, 6), dtype=complex)
1603 >>> smat[:6, :6] = np.diag(S)
1604 >>> np.allclose(a, np.dot(U, np.dot(smat, Vh)))
1605 True
1607 Reconstruction based on reduced SVD, 2D case:
1609 >>> U, S, Vh = np.linalg.svd(a, full_matrices=False)
1610 >>> U.shape, S.shape, Vh.shape
1611 ((9, 6), (6,), (6, 6))
1612 >>> np.allclose(a, np.dot(U * S, Vh))
1613 True
1614 >>> smat = np.diag(S)
1615 >>> np.allclose(a, np.dot(U, np.dot(smat, Vh)))
1616 True
1618 Reconstruction based on full SVD, 4D case:
1620 >>> U, S, Vh = np.linalg.svd(b, full_matrices=True)
1621 >>> U.shape, S.shape, Vh.shape
1622 ((2, 7, 8, 8), (2, 7, 3), (2, 7, 3, 3))
1623 >>> np.allclose(b, np.matmul(U[..., :3] * S[..., None, :], Vh))
1624 True
1625 >>> np.allclose(b, np.matmul(U[..., :3], S[..., None] * Vh))
1626 True
1628 Reconstruction based on reduced SVD, 4D case:
1630 >>> U, S, Vh = np.linalg.svd(b, full_matrices=False)
1631 >>> U.shape, S.shape, Vh.shape
1632 ((2, 7, 8, 3), (2, 7, 3), (2, 7, 3, 3))
1633 >>> np.allclose(b, np.matmul(U * S[..., None, :], Vh))
1634 True
1635 >>> np.allclose(b, np.matmul(U, S[..., None] * Vh))
1636 True
1638 """
1639 import numpy as _nx
1640 a, wrap = _makearray(a)
1642 if hermitian:
1643 # note: lapack svd returns eigenvalues with s ** 2 sorted descending,
1644 # but eig returns s sorted ascending, so we re-order the eigenvalues
1645 # and related arrays to have the correct order
1646 if compute_uv:
1647 s, u = eigh(a)
1648 sgn = sign(s)
1649 s = abs(s)
1650 sidx = argsort(s)[..., ::-1]
1651 sgn = _nx.take_along_axis(sgn, sidx, axis=-1)
1652 s = _nx.take_along_axis(s, sidx, axis=-1)
1653 u = _nx.take_along_axis(u, sidx[..., None, :], axis=-1)
1654 # singular values are unsigned, move the sign into v
1655 vt = transpose(u * sgn[..., None, :]).conjugate()
1656 return SVDResult(wrap(u), s, wrap(vt))
1657 else:
1658 s = eigvalsh(a)
1659 s = abs(s)
1660 return sort(s)[..., ::-1]
1662 _assert_stacked_2d(a)
1663 t, result_t = _commonType(a)
1665 extobj = get_linalg_error_extobj(_raise_linalgerror_svd_nonconvergence)
1667 m, n = a.shape[-2:]
1668 if compute_uv:
1669 if full_matrices:
1670 if m < n:
1671 gufunc = _umath_linalg.svd_m_f
1672 else:
1673 gufunc = _umath_linalg.svd_n_f
1674 else:
1675 if m < n:
1676 gufunc = _umath_linalg.svd_m_s
1677 else:
1678 gufunc = _umath_linalg.svd_n_s
1680 signature = 'D->DdD' if isComplexType(t) else 'd->ddd'
1681 u, s, vh = gufunc(a, signature=signature, extobj=extobj)
1682 u = u.astype(result_t, copy=False)
1683 s = s.astype(_realType(result_t), copy=False)
1684 vh = vh.astype(result_t, copy=False)
1685 return SVDResult(wrap(u), s, wrap(vh))
1686 else:
1687 if m < n:
1688 gufunc = _umath_linalg.svd_m
1689 else:
1690 gufunc = _umath_linalg.svd_n
1692 signature = 'D->d' if isComplexType(t) else 'd->d'
1693 s = gufunc(a, signature=signature, extobj=extobj)
1694 s = s.astype(_realType(result_t), copy=False)
1695 return s
1698def _cond_dispatcher(x, p=None):
1699 return (x,)
1702@array_function_dispatch(_cond_dispatcher)
1703def cond(x, p=None):
1704 """
1705 Compute the condition number of a matrix.
1707 This function is capable of returning the condition number using
1708 one of seven different norms, depending on the value of `p` (see
1709 Parameters below).
1711 Parameters
1712 ----------
1713 x : (..., M, N) array_like
1714 The matrix whose condition number is sought.
1715 p : {None, 1, -1, 2, -2, inf, -inf, 'fro'}, optional
1716 Order of the norm used in the condition number computation:
1718 ===== ============================
1719 p norm for matrices
1720 ===== ============================
1721 None 2-norm, computed directly using the ``SVD``
1722 'fro' Frobenius norm
1723 inf max(sum(abs(x), axis=1))
1724 -inf min(sum(abs(x), axis=1))
1725 1 max(sum(abs(x), axis=0))
1726 -1 min(sum(abs(x), axis=0))
1727 2 2-norm (largest sing. value)
1728 -2 smallest singular value
1729 ===== ============================
1731 inf means the `numpy.inf` object, and the Frobenius norm is
1732 the root-of-sum-of-squares norm.
1734 Returns
1735 -------
1736 c : {float, inf}
1737 The condition number of the matrix. May be infinite.
1739 See Also
1740 --------
1741 numpy.linalg.norm
1743 Notes
1744 -----
1745 The condition number of `x` is defined as the norm of `x` times the
1746 norm of the inverse of `x` [1]_; the norm can be the usual L2-norm
1747 (root-of-sum-of-squares) or one of a number of other matrix norms.
1749 References
1750 ----------
1751 .. [1] G. Strang, *Linear Algebra and Its Applications*, Orlando, FL,
1752 Academic Press, Inc., 1980, pg. 285.
1754 Examples
1755 --------
1756 >>> from numpy import linalg as LA
1757 >>> a = np.array([[1, 0, -1], [0, 1, 0], [1, 0, 1]])
1758 >>> a
1759 array([[ 1, 0, -1],
1760 [ 0, 1, 0],
1761 [ 1, 0, 1]])
1762 >>> LA.cond(a)
1763 1.4142135623730951
1764 >>> LA.cond(a, 'fro')
1765 3.1622776601683795
1766 >>> LA.cond(a, np.inf)
1767 2.0
1768 >>> LA.cond(a, -np.inf)
1769 1.0
1770 >>> LA.cond(a, 1)
1771 2.0
1772 >>> LA.cond(a, -1)
1773 1.0
1774 >>> LA.cond(a, 2)
1775 1.4142135623730951
1776 >>> LA.cond(a, -2)
1777 0.70710678118654746 # may vary
1778 >>> min(LA.svd(a, compute_uv=False))*min(LA.svd(LA.inv(a), compute_uv=False))
1779 0.70710678118654746 # may vary
1781 """
1782 x = asarray(x) # in case we have a matrix
1783 if _is_empty_2d(x):
1784 raise LinAlgError("cond is not defined on empty arrays")
1785 if p is None or p == 2 or p == -2:
1786 s = svd(x, compute_uv=False)
1787 with errstate(all='ignore'):
1788 if p == -2:
1789 r = s[..., -1] / s[..., 0]
1790 else:
1791 r = s[..., 0] / s[..., -1]
1792 else:
1793 # Call inv(x) ignoring errors. The result array will
1794 # contain nans in the entries where inversion failed.
1795 _assert_stacked_2d(x)
1796 _assert_stacked_square(x)
1797 t, result_t = _commonType(x)
1798 signature = 'D->D' if isComplexType(t) else 'd->d'
1799 with errstate(all='ignore'):
1800 invx = _umath_linalg.inv(x, signature=signature)
1801 r = norm(x, p, axis=(-2, -1)) * norm(invx, p, axis=(-2, -1))
1802 r = r.astype(result_t, copy=False)
1804 # Convert nans to infs unless the original array had nan entries
1805 r = asarray(r)
1806 nan_mask = isnan(r)
1807 if nan_mask.any():
1808 nan_mask &= ~isnan(x).any(axis=(-2, -1))
1809 if r.ndim > 0:
1810 r[nan_mask] = Inf
1811 elif nan_mask:
1812 r[()] = Inf
1814 # Convention is to return scalars instead of 0d arrays
1815 if r.ndim == 0:
1816 r = r[()]
1818 return r
1821def _matrix_rank_dispatcher(A, tol=None, hermitian=None):
1822 return (A,)
1825@array_function_dispatch(_matrix_rank_dispatcher)
1826def matrix_rank(A, tol=None, hermitian=False):
1827 """
1828 Return matrix rank of array using SVD method
1830 Rank of the array is the number of singular values of the array that are
1831 greater than `tol`.
1833 .. versionchanged:: 1.14
1834 Can now operate on stacks of matrices
1836 Parameters
1837 ----------
1838 A : {(M,), (..., M, N)} array_like
1839 Input vector or stack of matrices.
1840 tol : (...) array_like, float, optional
1841 Threshold below which SVD values are considered zero. If `tol` is
1842 None, and ``S`` is an array with singular values for `M`, and
1843 ``eps`` is the epsilon value for datatype of ``S``, then `tol` is
1844 set to ``S.max() * max(M, N) * eps``.
1846 .. versionchanged:: 1.14
1847 Broadcasted against the stack of matrices
1848 hermitian : bool, optional
1849 If True, `A` is assumed to be Hermitian (symmetric if real-valued),
1850 enabling a more efficient method for finding singular values.
1851 Defaults to False.
1853 .. versionadded:: 1.14
1855 Returns
1856 -------
1857 rank : (...) array_like
1858 Rank of A.
1860 Notes
1861 -----
1862 The default threshold to detect rank deficiency is a test on the magnitude
1863 of the singular values of `A`. By default, we identify singular values less
1864 than ``S.max() * max(M, N) * eps`` as indicating rank deficiency (with
1865 the symbols defined above). This is the algorithm MATLAB uses [1]. It also
1866 appears in *Numerical recipes* in the discussion of SVD solutions for linear
1867 least squares [2].
1869 This default threshold is designed to detect rank deficiency accounting for
1870 the numerical errors of the SVD computation. Imagine that there is a column
1871 in `A` that is an exact (in floating point) linear combination of other
1872 columns in `A`. Computing the SVD on `A` will not produce a singular value
1873 exactly equal to 0 in general: any difference of the smallest SVD value from
1874 0 will be caused by numerical imprecision in the calculation of the SVD.
1875 Our threshold for small SVD values takes this numerical imprecision into
1876 account, and the default threshold will detect such numerical rank
1877 deficiency. The threshold may declare a matrix `A` rank deficient even if
1878 the linear combination of some columns of `A` is not exactly equal to
1879 another column of `A` but only numerically very close to another column of
1880 `A`.
1882 We chose our default threshold because it is in wide use. Other thresholds
1883 are possible. For example, elsewhere in the 2007 edition of *Numerical
1884 recipes* there is an alternative threshold of ``S.max() *
1885 np.finfo(A.dtype).eps / 2. * np.sqrt(m + n + 1.)``. The authors describe
1886 this threshold as being based on "expected roundoff error" (p 71).
1888 The thresholds above deal with floating point roundoff error in the
1889 calculation of the SVD. However, you may have more information about the
1890 sources of error in `A` that would make you consider other tolerance values
1891 to detect *effective* rank deficiency. The most useful measure of the
1892 tolerance depends on the operations you intend to use on your matrix. For
1893 example, if your data come from uncertain measurements with uncertainties
1894 greater than floating point epsilon, choosing a tolerance near that
1895 uncertainty may be preferable. The tolerance may be absolute if the
1896 uncertainties are absolute rather than relative.
1898 References
1899 ----------
1900 .. [1] MATLAB reference documentation, "Rank"
1901 https://www.mathworks.com/help/techdoc/ref/rank.html
1902 .. [2] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery,
1903 "Numerical Recipes (3rd edition)", Cambridge University Press, 2007,
1904 page 795.
1906 Examples
1907 --------
1908 >>> from numpy.linalg import matrix_rank
1909 >>> matrix_rank(np.eye(4)) # Full rank matrix
1910 4
1911 >>> I=np.eye(4); I[-1,-1] = 0. # rank deficient matrix
1912 >>> matrix_rank(I)
1913 3
1914 >>> matrix_rank(np.ones((4,))) # 1 dimension - rank 1 unless all 0
1915 1
1916 >>> matrix_rank(np.zeros((4,)))
1917 0
1918 """
1919 A = asarray(A)
1920 if A.ndim < 2:
1921 return int(not all(A==0))
1922 S = svd(A, compute_uv=False, hermitian=hermitian)
1923 if tol is None:
1924 tol = S.max(axis=-1, keepdims=True) * max(A.shape[-2:]) * finfo(S.dtype).eps
1925 else:
1926 tol = asarray(tol)[..., newaxis]
1927 return count_nonzero(S > tol, axis=-1)
1930# Generalized inverse
1932def _pinv_dispatcher(a, rcond=None, hermitian=None):
1933 return (a,)
1936@array_function_dispatch(_pinv_dispatcher)
1937def pinv(a, rcond=1e-15, hermitian=False):
1938 """
1939 Compute the (Moore-Penrose) pseudo-inverse of a matrix.
1941 Calculate the generalized inverse of a matrix using its
1942 singular-value decomposition (SVD) and including all
1943 *large* singular values.
1945 .. versionchanged:: 1.14
1946 Can now operate on stacks of matrices
1948 Parameters
1949 ----------
1950 a : (..., M, N) array_like
1951 Matrix or stack of matrices to be pseudo-inverted.
1952 rcond : (...) array_like of float
1953 Cutoff for small singular values.
1954 Singular values less than or equal to
1955 ``rcond * largest_singular_value`` are set to zero.
1956 Broadcasts against the stack of matrices.
1957 hermitian : bool, optional
1958 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
1959 enabling a more efficient method for finding singular values.
1960 Defaults to False.
1962 .. versionadded:: 1.17.0
1964 Returns
1965 -------
1966 B : (..., N, M) ndarray
1967 The pseudo-inverse of `a`. If `a` is a `matrix` instance, then so
1968 is `B`.
1970 Raises
1971 ------
1972 LinAlgError
1973 If the SVD computation does not converge.
1975 See Also
1976 --------
1977 scipy.linalg.pinv : Similar function in SciPy.
1978 scipy.linalg.pinvh : Compute the (Moore-Penrose) pseudo-inverse of a
1979 Hermitian matrix.
1981 Notes
1982 -----
1983 The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
1984 defined as: "the matrix that 'solves' [the least-squares problem]
1985 :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then
1986 :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`.
1988 It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular
1989 value decomposition of A, then
1990 :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are
1991 orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting
1992 of A's so-called singular values, (followed, typically, by
1993 zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix
1994 consisting of the reciprocals of A's singular values
1995 (again, followed by zeros). [1]_
1997 References
1998 ----------
1999 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
2000 FL, Academic Press, Inc., 1980, pp. 139-142.
2002 Examples
2003 --------
2004 The following example checks that ``a * a+ * a == a`` and
2005 ``a+ * a * a+ == a+``:
2007 >>> a = np.random.randn(9, 6)
2008 >>> B = np.linalg.pinv(a)
2009 >>> np.allclose(a, np.dot(a, np.dot(B, a)))
2010 True
2011 >>> np.allclose(B, np.dot(B, np.dot(a, B)))
2012 True
2014 """
2015 a, wrap = _makearray(a)
2016 rcond = asarray(rcond)
2017 if _is_empty_2d(a):
2018 m, n = a.shape[-2:]
2019 res = empty(a.shape[:-2] + (n, m), dtype=a.dtype)
2020 return wrap(res)
2021 a = a.conjugate()
2022 u, s, vt = svd(a, full_matrices=False, hermitian=hermitian)
2024 # discard small singular values
2025 cutoff = rcond[..., newaxis] * amax(s, axis=-1, keepdims=True)
2026 large = s > cutoff
2027 s = divide(1, s, where=large, out=s)
2028 s[~large] = 0
2030 res = matmul(transpose(vt), multiply(s[..., newaxis], transpose(u)))
2031 return wrap(res)
2034# Determinant
2037@array_function_dispatch(_unary_dispatcher)
2038def slogdet(a):
2039 """
2040 Compute the sign and (natural) logarithm of the determinant of an array.
2042 If an array has a very small or very large determinant, then a call to
2043 `det` may overflow or underflow. This routine is more robust against such
2044 issues, because it computes the logarithm of the determinant rather than
2045 the determinant itself.
2047 Parameters
2048 ----------
2049 a : (..., M, M) array_like
2050 Input array, has to be a square 2-D array.
2052 Returns
2053 -------
2054 A namedtuple with the following attributes:
2056 sign : (...) array_like
2057 A number representing the sign of the determinant. For a real matrix,
2058 this is 1, 0, or -1. For a complex matrix, this is a complex number
2059 with absolute value 1 (i.e., it is on the unit circle), or else 0.
2060 logabsdet : (...) array_like
2061 The natural log of the absolute value of the determinant.
2063 If the determinant is zero, then `sign` will be 0 and `logabsdet` will be
2064 -Inf. In all cases, the determinant is equal to ``sign * np.exp(logabsdet)``.
2066 See Also
2067 --------
2068 det
2070 Notes
2071 -----
2073 .. versionadded:: 1.8.0
2075 Broadcasting rules apply, see the `numpy.linalg` documentation for
2076 details.
2078 .. versionadded:: 1.6.0
2080 The determinant is computed via LU factorization using the LAPACK
2081 routine ``z/dgetrf``.
2084 Examples
2085 --------
2086 The determinant of a 2-D array ``[[a, b], [c, d]]`` is ``ad - bc``:
2088 >>> a = np.array([[1, 2], [3, 4]])
2089 >>> (sign, logabsdet) = np.linalg.slogdet(a)
2090 >>> (sign, logabsdet)
2091 (-1, 0.69314718055994529) # may vary
2092 >>> sign * np.exp(logabsdet)
2093 -2.0
2095 Computing log-determinants for a stack of matrices:
2097 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2098 >>> a.shape
2099 (3, 2, 2)
2100 >>> sign, logabsdet = np.linalg.slogdet(a)
2101 >>> (sign, logabsdet)
2102 (array([-1., -1., -1.]), array([ 0.69314718, 1.09861229, 2.07944154]))
2103 >>> sign * np.exp(logabsdet)
2104 array([-2., -3., -8.])
2106 This routine succeeds where ordinary `det` does not:
2108 >>> np.linalg.det(np.eye(500) * 0.1)
2109 0.0
2110 >>> np.linalg.slogdet(np.eye(500) * 0.1)
2111 (1, -1151.2925464970228)
2113 """
2114 a = asarray(a)
2115 _assert_stacked_2d(a)
2116 _assert_stacked_square(a)
2117 t, result_t = _commonType(a)
2118 real_t = _realType(result_t)
2119 signature = 'D->Dd' if isComplexType(t) else 'd->dd'
2120 sign, logdet = _umath_linalg.slogdet(a, signature=signature)
2121 sign = sign.astype(result_t, copy=False)
2122 logdet = logdet.astype(real_t, copy=False)
2123 return SlogdetResult(sign, logdet)
2126@array_function_dispatch(_unary_dispatcher)
2127def det(a):
2128 """
2129 Compute the determinant of an array.
2131 Parameters
2132 ----------
2133 a : (..., M, M) array_like
2134 Input array to compute determinants for.
2136 Returns
2137 -------
2138 det : (...) array_like
2139 Determinant of `a`.
2141 See Also
2142 --------
2143 slogdet : Another way to represent the determinant, more suitable
2144 for large matrices where underflow/overflow may occur.
2145 scipy.linalg.det : Similar function in SciPy.
2147 Notes
2148 -----
2150 .. versionadded:: 1.8.0
2152 Broadcasting rules apply, see the `numpy.linalg` documentation for
2153 details.
2155 The determinant is computed via LU factorization using the LAPACK
2156 routine ``z/dgetrf``.
2158 Examples
2159 --------
2160 The determinant of a 2-D array [[a, b], [c, d]] is ad - bc:
2162 >>> a = np.array([[1, 2], [3, 4]])
2163 >>> np.linalg.det(a)
2164 -2.0 # may vary
2166 Computing determinants for a stack of matrices:
2168 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2169 >>> a.shape
2170 (3, 2, 2)
2171 >>> np.linalg.det(a)
2172 array([-2., -3., -8.])
2174 """
2175 a = asarray(a)
2176 _assert_stacked_2d(a)
2177 _assert_stacked_square(a)
2178 t, result_t = _commonType(a)
2179 signature = 'D->D' if isComplexType(t) else 'd->d'
2180 r = _umath_linalg.det(a, signature=signature)
2181 r = r.astype(result_t, copy=False)
2182 return r
2185# Linear Least Squares
2187def _lstsq_dispatcher(a, b, rcond=None):
2188 return (a, b)
2191@array_function_dispatch(_lstsq_dispatcher)
2192def lstsq(a, b, rcond="warn"):
2193 r"""
2194 Return the least-squares solution to a linear matrix equation.
2196 Computes the vector `x` that approximately solves the equation
2197 ``a @ x = b``. The equation may be under-, well-, or over-determined
2198 (i.e., the number of linearly independent rows of `a` can be less than,
2199 equal to, or greater than its number of linearly independent columns).
2200 If `a` is square and of full rank, then `x` (but for round-off error)
2201 is the "exact" solution of the equation. Else, `x` minimizes the
2202 Euclidean 2-norm :math:`||b - ax||`. If there are multiple minimizing
2203 solutions, the one with the smallest 2-norm :math:`||x||` is returned.
2205 Parameters
2206 ----------
2207 a : (M, N) array_like
2208 "Coefficient" matrix.
2209 b : {(M,), (M, K)} array_like
2210 Ordinate or "dependent variable" values. If `b` is two-dimensional,
2211 the least-squares solution is calculated for each of the `K` columns
2212 of `b`.
2213 rcond : float, optional
2214 Cut-off ratio for small singular values of `a`.
2215 For the purposes of rank determination, singular values are treated
2216 as zero if they are smaller than `rcond` times the largest singular
2217 value of `a`.
2219 .. versionchanged:: 1.14.0
2220 If not set, a FutureWarning is given. The previous default
2221 of ``-1`` will use the machine precision as `rcond` parameter,
2222 the new default will use the machine precision times `max(M, N)`.
2223 To silence the warning and use the new default, use ``rcond=None``,
2224 to keep using the old behavior, use ``rcond=-1``.
2226 Returns
2227 -------
2228 x : {(N,), (N, K)} ndarray
2229 Least-squares solution. If `b` is two-dimensional,
2230 the solutions are in the `K` columns of `x`.
2231 residuals : {(1,), (K,), (0,)} ndarray
2232 Sums of squared residuals: Squared Euclidean 2-norm for each column in
2233 ``b - a @ x``.
2234 If the rank of `a` is < N or M <= N, this is an empty array.
2235 If `b` is 1-dimensional, this is a (1,) shape array.
2236 Otherwise the shape is (K,).
2237 rank : int
2238 Rank of matrix `a`.
2239 s : (min(M, N),) ndarray
2240 Singular values of `a`.
2242 Raises
2243 ------
2244 LinAlgError
2245 If computation does not converge.
2247 See Also
2248 --------
2249 scipy.linalg.lstsq : Similar function in SciPy.
2251 Notes
2252 -----
2253 If `b` is a matrix, then all array results are returned as matrices.
2255 Examples
2256 --------
2257 Fit a line, ``y = mx + c``, through some noisy data-points:
2259 >>> x = np.array([0, 1, 2, 3])
2260 >>> y = np.array([-1, 0.2, 0.9, 2.1])
2262 By examining the coefficients, we see that the line should have a
2263 gradient of roughly 1 and cut the y-axis at, more or less, -1.
2265 We can rewrite the line equation as ``y = Ap``, where ``A = [[x 1]]``
2266 and ``p = [[m], [c]]``. Now use `lstsq` to solve for `p`:
2268 >>> A = np.vstack([x, np.ones(len(x))]).T
2269 >>> A
2270 array([[ 0., 1.],
2271 [ 1., 1.],
2272 [ 2., 1.],
2273 [ 3., 1.]])
2275 >>> m, c = np.linalg.lstsq(A, y, rcond=None)[0]
2276 >>> m, c
2277 (1.0 -0.95) # may vary
2279 Plot the data along with the fitted line:
2281 >>> import matplotlib.pyplot as plt
2282 >>> _ = plt.plot(x, y, 'o', label='Original data', markersize=10)
2283 >>> _ = plt.plot(x, m*x + c, 'r', label='Fitted line')
2284 >>> _ = plt.legend()
2285 >>> plt.show()
2287 """
2288 a, _ = _makearray(a)
2289 b, wrap = _makearray(b)
2290 is_1d = b.ndim == 1
2291 if is_1d:
2292 b = b[:, newaxis]
2293 _assert_2d(a, b)
2294 m, n = a.shape[-2:]
2295 m2, n_rhs = b.shape[-2:]
2296 if m != m2:
2297 raise LinAlgError('Incompatible dimensions')
2299 t, result_t = _commonType(a, b)
2300 result_real_t = _realType(result_t)
2302 # Determine default rcond value
2303 if rcond == "warn":
2304 # 2017-08-19, 1.14.0
2305 warnings.warn("`rcond` parameter will change to the default of "
2306 "machine precision times ``max(M, N)`` where M and N "
2307 "are the input matrix dimensions.\n"
2308 "To use the future default and silence this warning "
2309 "we advise to pass `rcond=None`, to keep using the old, "
2310 "explicitly pass `rcond=-1`.",
2311 FutureWarning, stacklevel=2)
2312 rcond = -1
2313 if rcond is None:
2314 rcond = finfo(t).eps * max(n, m)
2316 if m <= n:
2317 gufunc = _umath_linalg.lstsq_m
2318 else:
2319 gufunc = _umath_linalg.lstsq_n
2321 signature = 'DDd->Ddid' if isComplexType(t) else 'ddd->ddid'
2322 extobj = get_linalg_error_extobj(_raise_linalgerror_lstsq)
2323 if n_rhs == 0:
2324 # lapack can't handle n_rhs = 0 - so allocate the array one larger in that axis
2325 b = zeros(b.shape[:-2] + (m, n_rhs + 1), dtype=b.dtype)
2326 x, resids, rank, s = gufunc(a, b, rcond, signature=signature, extobj=extobj)
2327 if m == 0:
2328 x[...] = 0
2329 if n_rhs == 0:
2330 # remove the item we added
2331 x = x[..., :n_rhs]
2332 resids = resids[..., :n_rhs]
2334 # remove the axis we added
2335 if is_1d:
2336 x = x.squeeze(axis=-1)
2337 # we probably should squeeze resids too, but we can't
2338 # without breaking compatibility.
2340 # as documented
2341 if rank != n or m <= n:
2342 resids = array([], result_real_t)
2344 # coerce output arrays
2345 s = s.astype(result_real_t, copy=False)
2346 resids = resids.astype(result_real_t, copy=False)
2347 x = x.astype(result_t, copy=True) # Copying lets the memory in r_parts be freed
2348 return wrap(x), wrap(resids), rank, s
2351def _multi_svd_norm(x, row_axis, col_axis, op):
2352 """Compute a function of the singular values of the 2-D matrices in `x`.
2354 This is a private utility function used by `numpy.linalg.norm()`.
2356 Parameters
2357 ----------
2358 x : ndarray
2359 row_axis, col_axis : int
2360 The axes of `x` that hold the 2-D matrices.
2361 op : callable
2362 This should be either numpy.amin or `numpy.amax` or `numpy.sum`.
2364 Returns
2365 -------
2366 result : float or ndarray
2367 If `x` is 2-D, the return values is a float.
2368 Otherwise, it is an array with ``x.ndim - 2`` dimensions.
2369 The return values are either the minimum or maximum or sum of the
2370 singular values of the matrices, depending on whether `op`
2371 is `numpy.amin` or `numpy.amax` or `numpy.sum`.
2373 """
2374 y = moveaxis(x, (row_axis, col_axis), (-2, -1))
2375 result = op(svd(y, compute_uv=False), axis=-1)
2376 return result
2379def _norm_dispatcher(x, ord=None, axis=None, keepdims=None):
2380 return (x,)
2383@array_function_dispatch(_norm_dispatcher)
2384def norm(x, ord=None, axis=None, keepdims=False):
2385 """
2386 Matrix or vector norm.
2388 This function is able to return one of eight different matrix norms,
2389 or one of an infinite number of vector norms (described below), depending
2390 on the value of the ``ord`` parameter.
2392 Parameters
2393 ----------
2394 x : array_like
2395 Input array. If `axis` is None, `x` must be 1-D or 2-D, unless `ord`
2396 is None. If both `axis` and `ord` are None, the 2-norm of
2397 ``x.ravel`` will be returned.
2398 ord : {non-zero int, inf, -inf, 'fro', 'nuc'}, optional
2399 Order of the norm (see table under ``Notes``). inf means numpy's
2400 `inf` object. The default is None.
2401 axis : {None, int, 2-tuple of ints}, optional.
2402 If `axis` is an integer, it specifies the axis of `x` along which to
2403 compute the vector norms. If `axis` is a 2-tuple, it specifies the
2404 axes that hold 2-D matrices, and the matrix norms of these matrices
2405 are computed. If `axis` is None then either a vector norm (when `x`
2406 is 1-D) or a matrix norm (when `x` is 2-D) is returned. The default
2407 is None.
2409 .. versionadded:: 1.8.0
2411 keepdims : bool, optional
2412 If this is set to True, the axes which are normed over are left in the
2413 result as dimensions with size one. With this option the result will
2414 broadcast correctly against the original `x`.
2416 .. versionadded:: 1.10.0
2418 Returns
2419 -------
2420 n : float or ndarray
2421 Norm of the matrix or vector(s).
2423 See Also
2424 --------
2425 scipy.linalg.norm : Similar function in SciPy.
2427 Notes
2428 -----
2429 For values of ``ord < 1``, the result is, strictly speaking, not a
2430 mathematical 'norm', but it may still be useful for various numerical
2431 purposes.
2433 The following norms can be calculated:
2435 ===== ============================ ==========================
2436 ord norm for matrices norm for vectors
2437 ===== ============================ ==========================
2438 None Frobenius norm 2-norm
2439 'fro' Frobenius norm --
2440 'nuc' nuclear norm --
2441 inf max(sum(abs(x), axis=1)) max(abs(x))
2442 -inf min(sum(abs(x), axis=1)) min(abs(x))
2443 0 -- sum(x != 0)
2444 1 max(sum(abs(x), axis=0)) as below
2445 -1 min(sum(abs(x), axis=0)) as below
2446 2 2-norm (largest sing. value) as below
2447 -2 smallest singular value as below
2448 other -- sum(abs(x)**ord)**(1./ord)
2449 ===== ============================ ==========================
2451 The Frobenius norm is given by [1]_:
2453 :math:`||A||_F = [\\sum_{i,j} abs(a_{i,j})^2]^{1/2}`
2455 The nuclear norm is the sum of the singular values.
2457 Both the Frobenius and nuclear norm orders are only defined for
2458 matrices and raise a ValueError when ``x.ndim != 2``.
2460 References
2461 ----------
2462 .. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*,
2463 Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15
2465 Examples
2466 --------
2467 >>> from numpy import linalg as LA
2468 >>> a = np.arange(9) - 4
2469 >>> a
2470 array([-4, -3, -2, ..., 2, 3, 4])
2471 >>> b = a.reshape((3, 3))
2472 >>> b
2473 array([[-4, -3, -2],
2474 [-1, 0, 1],
2475 [ 2, 3, 4]])
2477 >>> LA.norm(a)
2478 7.745966692414834
2479 >>> LA.norm(b)
2480 7.745966692414834
2481 >>> LA.norm(b, 'fro')
2482 7.745966692414834
2483 >>> LA.norm(a, np.inf)
2484 4.0
2485 >>> LA.norm(b, np.inf)
2486 9.0
2487 >>> LA.norm(a, -np.inf)
2488 0.0
2489 >>> LA.norm(b, -np.inf)
2490 2.0
2492 >>> LA.norm(a, 1)
2493 20.0
2494 >>> LA.norm(b, 1)
2495 7.0
2496 >>> LA.norm(a, -1)
2497 -4.6566128774142013e-010
2498 >>> LA.norm(b, -1)
2499 6.0
2500 >>> LA.norm(a, 2)
2501 7.745966692414834
2502 >>> LA.norm(b, 2)
2503 7.3484692283495345
2505 >>> LA.norm(a, -2)
2506 0.0
2507 >>> LA.norm(b, -2)
2508 1.8570331885190563e-016 # may vary
2509 >>> LA.norm(a, 3)
2510 5.8480354764257312 # may vary
2511 >>> LA.norm(a, -3)
2512 0.0
2514 Using the `axis` argument to compute vector norms:
2516 >>> c = np.array([[ 1, 2, 3],
2517 ... [-1, 1, 4]])
2518 >>> LA.norm(c, axis=0)
2519 array([ 1.41421356, 2.23606798, 5. ])
2520 >>> LA.norm(c, axis=1)
2521 array([ 3.74165739, 4.24264069])
2522 >>> LA.norm(c, ord=1, axis=1)
2523 array([ 6., 6.])
2525 Using the `axis` argument to compute matrix norms:
2527 >>> m = np.arange(8).reshape(2,2,2)
2528 >>> LA.norm(m, axis=(1,2))
2529 array([ 3.74165739, 11.22497216])
2530 >>> LA.norm(m[0, :, :]), LA.norm(m[1, :, :])
2531 (3.7416573867739413, 11.224972160321824)
2533 """
2534 x = asarray(x)
2536 if not issubclass(x.dtype.type, (inexact, object_)):
2537 x = x.astype(float)
2539 # Immediately handle some default, simple, fast, and common cases.
2540 if axis is None:
2541 ndim = x.ndim
2542 if ((ord is None) or
2543 (ord in ('f', 'fro') and ndim == 2) or
2544 (ord == 2 and ndim == 1)):
2546 x = x.ravel(order='K')
2547 if isComplexType(x.dtype.type):
2548 x_real = x.real
2549 x_imag = x.imag
2550 sqnorm = x_real.dot(x_real) + x_imag.dot(x_imag)
2551 else:
2552 sqnorm = x.dot(x)
2553 ret = sqrt(sqnorm)
2554 if keepdims:
2555 ret = ret.reshape(ndim*[1])
2556 return ret
2558 # Normalize the `axis` argument to a tuple.
2559 nd = x.ndim
2560 if axis is None:
2561 axis = tuple(range(nd))
2562 elif not isinstance(axis, tuple):
2563 try:
2564 axis = int(axis)
2565 except Exception as e:
2566 raise TypeError("'axis' must be None, an integer or a tuple of integers") from e
2567 axis = (axis,)
2569 if len(axis) == 1:
2570 if ord == Inf:
2571 return abs(x).max(axis=axis, keepdims=keepdims)
2572 elif ord == -Inf:
2573 return abs(x).min(axis=axis, keepdims=keepdims)
2574 elif ord == 0:
2575 # Zero norm
2576 return (x != 0).astype(x.real.dtype).sum(axis=axis, keepdims=keepdims)
2577 elif ord == 1:
2578 # special case for speedup
2579 return add.reduce(abs(x), axis=axis, keepdims=keepdims)
2580 elif ord is None or ord == 2:
2581 # special case for speedup
2582 s = (x.conj() * x).real
2583 return sqrt(add.reduce(s, axis=axis, keepdims=keepdims))
2584 # None of the str-type keywords for ord ('fro', 'nuc')
2585 # are valid for vectors
2586 elif isinstance(ord, str):
2587 raise ValueError(f"Invalid norm order '{ord}' for vectors")
2588 else:
2589 absx = abs(x)
2590 absx **= ord
2591 ret = add.reduce(absx, axis=axis, keepdims=keepdims)
2592 ret **= reciprocal(ord, dtype=ret.dtype)
2593 return ret
2594 elif len(axis) == 2:
2595 row_axis, col_axis = axis
2596 row_axis = normalize_axis_index(row_axis, nd)
2597 col_axis = normalize_axis_index(col_axis, nd)
2598 if row_axis == col_axis:
2599 raise ValueError('Duplicate axes given.')
2600 if ord == 2:
2601 ret = _multi_svd_norm(x, row_axis, col_axis, amax)
2602 elif ord == -2:
2603 ret = _multi_svd_norm(x, row_axis, col_axis, amin)
2604 elif ord == 1:
2605 if col_axis > row_axis:
2606 col_axis -= 1
2607 ret = add.reduce(abs(x), axis=row_axis).max(axis=col_axis)
2608 elif ord == Inf:
2609 if row_axis > col_axis:
2610 row_axis -= 1
2611 ret = add.reduce(abs(x), axis=col_axis).max(axis=row_axis)
2612 elif ord == -1:
2613 if col_axis > row_axis:
2614 col_axis -= 1
2615 ret = add.reduce(abs(x), axis=row_axis).min(axis=col_axis)
2616 elif ord == -Inf:
2617 if row_axis > col_axis:
2618 row_axis -= 1
2619 ret = add.reduce(abs(x), axis=col_axis).min(axis=row_axis)
2620 elif ord in [None, 'fro', 'f']:
2621 ret = sqrt(add.reduce((x.conj() * x).real, axis=axis))
2622 elif ord == 'nuc':
2623 ret = _multi_svd_norm(x, row_axis, col_axis, sum)
2624 else:
2625 raise ValueError("Invalid norm order for matrices.")
2626 if keepdims:
2627 ret_shape = list(x.shape)
2628 ret_shape[axis[0]] = 1
2629 ret_shape[axis[1]] = 1
2630 ret = ret.reshape(ret_shape)
2631 return ret
2632 else:
2633 raise ValueError("Improper number of dimensions to norm.")
2636# multi_dot
2638def _multidot_dispatcher(arrays, *, out=None):
2639 yield from arrays
2640 yield out
2643@array_function_dispatch(_multidot_dispatcher)
2644def multi_dot(arrays, *, out=None):
2645 """
2646 Compute the dot product of two or more arrays in a single function call,
2647 while automatically selecting the fastest evaluation order.
2649 `multi_dot` chains `numpy.dot` and uses optimal parenthesization
2650 of the matrices [1]_ [2]_. Depending on the shapes of the matrices,
2651 this can speed up the multiplication a lot.
2653 If the first argument is 1-D it is treated as a row vector.
2654 If the last argument is 1-D it is treated as a column vector.
2655 The other arguments must be 2-D.
2657 Think of `multi_dot` as::
2659 def multi_dot(arrays): return functools.reduce(np.dot, arrays)
2662 Parameters
2663 ----------
2664 arrays : sequence of array_like
2665 If the first argument is 1-D it is treated as row vector.
2666 If the last argument is 1-D it is treated as column vector.
2667 The other arguments must be 2-D.
2668 out : ndarray, optional
2669 Output argument. This must have the exact kind that would be returned
2670 if it was not used. In particular, it must have the right type, must be
2671 C-contiguous, and its dtype must be the dtype that would be returned
2672 for `dot(a, b)`. This is a performance feature. Therefore, if these
2673 conditions are not met, an exception is raised, instead of attempting
2674 to be flexible.
2676 .. versionadded:: 1.19.0
2678 Returns
2679 -------
2680 output : ndarray
2681 Returns the dot product of the supplied arrays.
2683 See Also
2684 --------
2685 numpy.dot : dot multiplication with two arguments.
2687 References
2688 ----------
2690 .. [1] Cormen, "Introduction to Algorithms", Chapter 15.2, p. 370-378
2691 .. [2] https://en.wikipedia.org/wiki/Matrix_chain_multiplication
2693 Examples
2694 --------
2695 `multi_dot` allows you to write::
2697 >>> from numpy.linalg import multi_dot
2698 >>> # Prepare some data
2699 >>> A = np.random.random((10000, 100))
2700 >>> B = np.random.random((100, 1000))
2701 >>> C = np.random.random((1000, 5))
2702 >>> D = np.random.random((5, 333))
2703 >>> # the actual dot multiplication
2704 >>> _ = multi_dot([A, B, C, D])
2706 instead of::
2708 >>> _ = np.dot(np.dot(np.dot(A, B), C), D)
2709 >>> # or
2710 >>> _ = A.dot(B).dot(C).dot(D)
2712 Notes
2713 -----
2714 The cost for a matrix multiplication can be calculated with the
2715 following function::
2717 def cost(A, B):
2718 return A.shape[0] * A.shape[1] * B.shape[1]
2720 Assume we have three matrices
2721 :math:`A_{10x100}, B_{100x5}, C_{5x50}`.
2723 The costs for the two different parenthesizations are as follows::
2725 cost((AB)C) = 10*100*5 + 10*5*50 = 5000 + 2500 = 7500
2726 cost(A(BC)) = 10*100*50 + 100*5*50 = 50000 + 25000 = 75000
2728 """
2729 n = len(arrays)
2730 # optimization only makes sense for len(arrays) > 2
2731 if n < 2:
2732 raise ValueError("Expecting at least two arrays.")
2733 elif n == 2:
2734 return dot(arrays[0], arrays[1], out=out)
2736 arrays = [asanyarray(a) for a in arrays]
2738 # save original ndim to reshape the result array into the proper form later
2739 ndim_first, ndim_last = arrays[0].ndim, arrays[-1].ndim
2740 # Explicitly convert vectors to 2D arrays to keep the logic of the internal
2741 # _multi_dot_* functions as simple as possible.
2742 if arrays[0].ndim == 1:
2743 arrays[0] = atleast_2d(arrays[0])
2744 if arrays[-1].ndim == 1:
2745 arrays[-1] = atleast_2d(arrays[-1]).T
2746 _assert_2d(*arrays)
2748 # _multi_dot_three is much faster than _multi_dot_matrix_chain_order
2749 if n == 3:
2750 result = _multi_dot_three(arrays[0], arrays[1], arrays[2], out=out)
2751 else:
2752 order = _multi_dot_matrix_chain_order(arrays)
2753 result = _multi_dot(arrays, order, 0, n - 1, out=out)
2755 # return proper shape
2756 if ndim_first == 1 and ndim_last == 1:
2757 return result[0, 0] # scalar
2758 elif ndim_first == 1 or ndim_last == 1:
2759 return result.ravel() # 1-D
2760 else:
2761 return result
2764def _multi_dot_three(A, B, C, out=None):
2765 """
2766 Find the best order for three arrays and do the multiplication.
2768 For three arguments `_multi_dot_three` is approximately 15 times faster
2769 than `_multi_dot_matrix_chain_order`
2771 """
2772 a0, a1b0 = A.shape
2773 b1c0, c1 = C.shape
2774 # cost1 = cost((AB)C) = a0*a1b0*b1c0 + a0*b1c0*c1
2775 cost1 = a0 * b1c0 * (a1b0 + c1)
2776 # cost2 = cost(A(BC)) = a1b0*b1c0*c1 + a0*a1b0*c1
2777 cost2 = a1b0 * c1 * (a0 + b1c0)
2779 if cost1 < cost2:
2780 return dot(dot(A, B), C, out=out)
2781 else:
2782 return dot(A, dot(B, C), out=out)
2785def _multi_dot_matrix_chain_order(arrays, return_costs=False):
2786 """
2787 Return a np.array that encodes the optimal order of mutiplications.
2789 The optimal order array is then used by `_multi_dot()` to do the
2790 multiplication.
2792 Also return the cost matrix if `return_costs` is `True`
2794 The implementation CLOSELY follows Cormen, "Introduction to Algorithms",
2795 Chapter 15.2, p. 370-378. Note that Cormen uses 1-based indices.
2797 cost[i, j] = min([
2798 cost[prefix] + cost[suffix] + cost_mult(prefix, suffix)
2799 for k in range(i, j)])
2801 """
2802 n = len(arrays)
2803 # p stores the dimensions of the matrices
2804 # Example for p: A_{10x100}, B_{100x5}, C_{5x50} --> p = [10, 100, 5, 50]
2805 p = [a.shape[0] for a in arrays] + [arrays[-1].shape[1]]
2806 # m is a matrix of costs of the subproblems
2807 # m[i,j]: min number of scalar multiplications needed to compute A_{i..j}
2808 m = zeros((n, n), dtype=double)
2809 # s is the actual ordering
2810 # s[i, j] is the value of k at which we split the product A_i..A_j
2811 s = empty((n, n), dtype=intp)
2813 for l in range(1, n):
2814 for i in range(n - l):
2815 j = i + l
2816 m[i, j] = Inf
2817 for k in range(i, j):
2818 q = m[i, k] + m[k+1, j] + p[i]*p[k+1]*p[j+1]
2819 if q < m[i, j]:
2820 m[i, j] = q
2821 s[i, j] = k # Note that Cormen uses 1-based index
2823 return (s, m) if return_costs else s
2826def _multi_dot(arrays, order, i, j, out=None):
2827 """Actually do the multiplication with the given order."""
2828 if i == j:
2829 # the initial call with non-None out should never get here
2830 assert out is None
2832 return arrays[i]
2833 else:
2834 return dot(_multi_dot(arrays, order, i, order[i, j]),
2835 _multi_dot(arrays, order, order[i, j] + 1, j),
2836 out=out)