Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/numpy/linalg/linalg.py: 16%
651 statements
« prev ^ index » next coverage.py v7.4.0, created at 2024-01-03 07:57 +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
21from numpy.core import (
22 array, asarray, zeros, empty, empty_like, intc, single, double,
23 csingle, cdouble, inexact, complexfloating, newaxis, all, Inf, dot,
24 add, multiply, sqrt, sum, isfinite,
25 finfo, errstate, geterrobj, moveaxis, amin, amax, product, abs,
26 atleast_2d, intp, asanyarray, object_, matmul,
27 swapaxes, divide, count_nonzero, isnan, sign, argsort, sort,
28 reciprocal
29)
30from numpy.core.multiarray import normalize_axis_index
31from numpy.core.overrides import set_module
32from numpy.core import overrides
33from numpy.lib.twodim_base import triu, eye
34from numpy.linalg import _umath_linalg
37array_function_dispatch = functools.partial(
38 overrides.array_function_dispatch, module='numpy.linalg')
41fortran_int = intc
44@set_module('numpy.linalg')
45class LinAlgError(Exception):
46 """
47 Generic Python-exception-derived object raised by linalg functions.
49 General purpose exception class, derived from Python's exception.Exception
50 class, programmatically raised in linalg functions when a Linear
51 Algebra-related condition would prevent further correct execution of the
52 function.
54 Parameters
55 ----------
56 None
58 Examples
59 --------
60 >>> from numpy import linalg as LA
61 >>> LA.inv(np.zeros((2,2)))
62 Traceback (most recent call last):
63 File "<stdin>", line 1, in <module>
64 File "...linalg.py", line 350,
65 in inv return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))
66 File "...linalg.py", line 249,
67 in solve
68 raise LinAlgError('Singular matrix')
69 numpy.linalg.LinAlgError: Singular matrix
71 """
74def _determine_error_states():
75 errobj = geterrobj()
76 bufsize = errobj[0]
78 with errstate(invalid='call', over='ignore',
79 divide='ignore', under='ignore'):
80 invalid_call_errmask = geterrobj()[1]
82 return [bufsize, invalid_call_errmask, None]
84# Dealing with errors in _umath_linalg
85_linalg_error_extobj = _determine_error_states()
86del _determine_error_states
88def _raise_linalgerror_singular(err, flag):
89 raise LinAlgError("Singular matrix")
91def _raise_linalgerror_nonposdef(err, flag):
92 raise LinAlgError("Matrix is not positive definite")
94def _raise_linalgerror_eigenvalues_nonconvergence(err, flag):
95 raise LinAlgError("Eigenvalues did not converge")
97def _raise_linalgerror_svd_nonconvergence(err, flag):
98 raise LinAlgError("SVD did not converge")
100def _raise_linalgerror_lstsq(err, flag):
101 raise LinAlgError("SVD did not converge in Linear Least Squares")
103def _raise_linalgerror_qr(err, flag):
104 raise LinAlgError("Incorrect argument found while performing "
105 "QR factorization")
107def get_linalg_error_extobj(callback):
108 extobj = list(_linalg_error_extobj) # make a copy
109 extobj[2] = callback
110 return extobj
112def _makearray(a):
113 new = asarray(a)
114 wrap = getattr(a, "__array_prepare__", new.__array_wrap__)
115 return new, wrap
117def isComplexType(t):
118 return issubclass(t, complexfloating)
120_real_types_map = {single : single,
121 double : double,
122 csingle : single,
123 cdouble : double}
125_complex_types_map = {single : csingle,
126 double : cdouble,
127 csingle : csingle,
128 cdouble : cdouble}
130def _realType(t, default=double):
131 return _real_types_map.get(t, default)
133def _complexType(t, default=cdouble):
134 return _complex_types_map.get(t, default)
136def _commonType(*arrays):
137 # in lite version, use higher precision (always double or cdouble)
138 result_type = single
139 is_complex = False
140 for a in arrays:
141 if issubclass(a.dtype.type, inexact):
142 if isComplexType(a.dtype.type):
143 is_complex = True
144 rt = _realType(a.dtype.type, default=None)
145 if rt is None:
146 # unsupported inexact scalar
147 raise TypeError("array type %s is unsupported in linalg" %
148 (a.dtype.name,))
149 else:
150 rt = double
151 if rt is double:
152 result_type = double
153 if is_complex:
154 t = cdouble
155 result_type = _complex_types_map[result_type]
156 else:
157 t = double
158 return t, result_type
161def _to_native_byte_order(*arrays):
162 ret = []
163 for arr in arrays:
164 if arr.dtype.byteorder not in ('=', '|'):
165 ret.append(asarray(arr, dtype=arr.dtype.newbyteorder('=')))
166 else:
167 ret.append(arr)
168 if len(ret) == 1:
169 return ret[0]
170 else:
171 return ret
174def _assert_2d(*arrays):
175 for a in arrays:
176 if a.ndim != 2:
177 raise LinAlgError('%d-dimensional array given. Array must be '
178 'two-dimensional' % a.ndim)
180def _assert_stacked_2d(*arrays):
181 for a in arrays:
182 if a.ndim < 2:
183 raise LinAlgError('%d-dimensional array given. Array must be '
184 'at least two-dimensional' % a.ndim)
186def _assert_stacked_square(*arrays):
187 for a in arrays:
188 m, n = a.shape[-2:]
189 if m != n:
190 raise LinAlgError('Last 2 dimensions of the array must be square')
192def _assert_finite(*arrays):
193 for a in arrays:
194 if not isfinite(a).all():
195 raise LinAlgError("Array must not contain infs or NaNs")
197def _is_empty_2d(arr):
198 # check size first for efficiency
199 return arr.size == 0 and product(arr.shape[-2:]) == 0
202def transpose(a):
203 """
204 Transpose each matrix in a stack of matrices.
206 Unlike np.transpose, this only swaps the last two axes, rather than all of
207 them
209 Parameters
210 ----------
211 a : (...,M,N) array_like
213 Returns
214 -------
215 aT : (...,N,M) ndarray
216 """
217 return swapaxes(a, -1, -2)
219# Linear equations
221def _tensorsolve_dispatcher(a, b, axes=None):
222 return (a, b)
225@array_function_dispatch(_tensorsolve_dispatcher)
226def tensorsolve(a, b, axes=None):
227 """
228 Solve the tensor equation ``a x = b`` for x.
230 It is assumed that all indices of `x` are summed over in the product,
231 together with the rightmost indices of `a`, as is done in, for example,
232 ``tensordot(a, x, axes=x.ndim)``.
234 Parameters
235 ----------
236 a : array_like
237 Coefficient tensor, of shape ``b.shape + Q``. `Q`, a tuple, equals
238 the shape of that sub-tensor of `a` consisting of the appropriate
239 number of its rightmost indices, and must be such that
240 ``prod(Q) == prod(b.shape)`` (in which sense `a` is said to be
241 'square').
242 b : array_like
243 Right-hand tensor, which can be of any shape.
244 axes : tuple of ints, optional
245 Axes in `a` to reorder to the right, before inversion.
246 If None (default), no reordering is done.
248 Returns
249 -------
250 x : ndarray, shape Q
252 Raises
253 ------
254 LinAlgError
255 If `a` is singular or not 'square' (in the above sense).
257 See Also
258 --------
259 numpy.tensordot, tensorinv, numpy.einsum
261 Examples
262 --------
263 >>> a = np.eye(2*3*4)
264 >>> a.shape = (2*3, 4, 2, 3, 4)
265 >>> b = np.random.randn(2*3, 4)
266 >>> x = np.linalg.tensorsolve(a, b)
267 >>> x.shape
268 (2, 3, 4)
269 >>> np.allclose(np.tensordot(a, x, axes=3), b)
270 True
272 """
273 a, wrap = _makearray(a)
274 b = asarray(b)
275 an = a.ndim
277 if axes is not None:
278 allaxes = list(range(0, an))
279 for k in axes:
280 allaxes.remove(k)
281 allaxes.insert(an, k)
282 a = a.transpose(allaxes)
284 oldshape = a.shape[-(an-b.ndim):]
285 prod = 1
286 for k in oldshape:
287 prod *= k
289 if a.size != prod ** 2:
290 raise LinAlgError(
291 "Input arrays must satisfy the requirement \
292 prod(a.shape[b.ndim:]) == prod(a.shape[:b.ndim])"
293 )
295 a = a.reshape(prod, prod)
296 b = b.ravel()
297 res = wrap(solve(a, b))
298 res.shape = oldshape
299 return res
302def _solve_dispatcher(a, b):
303 return (a, b)
306@array_function_dispatch(_solve_dispatcher)
307def solve(a, b):
308 """
309 Solve a linear matrix equation, or system of linear scalar equations.
311 Computes the "exact" solution, `x`, of the well-determined, i.e., full
312 rank, linear matrix equation `ax = b`.
314 Parameters
315 ----------
316 a : (..., M, M) array_like
317 Coefficient matrix.
318 b : {(..., M,), (..., M, K)}, array_like
319 Ordinate or "dependent variable" values.
321 Returns
322 -------
323 x : {(..., M,), (..., M, K)} ndarray
324 Solution to the system a x = b. Returned shape is identical to `b`.
326 Raises
327 ------
328 LinAlgError
329 If `a` is singular or not square.
331 See Also
332 --------
333 scipy.linalg.solve : Similar function in SciPy.
335 Notes
336 -----
338 .. versionadded:: 1.8.0
340 Broadcasting rules apply, see the `numpy.linalg` documentation for
341 details.
343 The solutions are computed using LAPACK routine ``_gesv``.
345 `a` must be square and of full-rank, i.e., all rows (or, equivalently,
346 columns) must be linearly independent; if either is not true, use
347 `lstsq` for the least-squares best "solution" of the
348 system/equation.
350 References
351 ----------
352 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
353 FL, Academic Press, Inc., 1980, pg. 22.
355 Examples
356 --------
357 Solve the system of equations ``x0 + 2 * x1 = 1`` and ``3 * x0 + 5 * x1 = 2``:
359 >>> a = np.array([[1, 2], [3, 5]])
360 >>> b = np.array([1, 2])
361 >>> x = np.linalg.solve(a, b)
362 >>> x
363 array([-1., 1.])
365 Check that the solution is correct:
367 >>> np.allclose(np.dot(a, x), b)
368 True
370 """
371 a, _ = _makearray(a)
372 _assert_stacked_2d(a)
373 _assert_stacked_square(a)
374 b, wrap = _makearray(b)
375 t, result_t = _commonType(a, b)
377 # We use the b = (..., M,) logic, only if the number of extra dimensions
378 # match exactly
379 if b.ndim == a.ndim - 1:
380 gufunc = _umath_linalg.solve1
381 else:
382 gufunc = _umath_linalg.solve
384 signature = 'DD->D' if isComplexType(t) else 'dd->d'
385 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
386 r = gufunc(a, b, signature=signature, extobj=extobj)
388 return wrap(r.astype(result_t, copy=False))
391def _tensorinv_dispatcher(a, ind=None):
392 return (a,)
395@array_function_dispatch(_tensorinv_dispatcher)
396def tensorinv(a, ind=2):
397 """
398 Compute the 'inverse' of an N-dimensional array.
400 The result is an inverse for `a` relative to the tensordot operation
401 ``tensordot(a, b, ind)``, i. e., up to floating-point accuracy,
402 ``tensordot(tensorinv(a), a, ind)`` is the "identity" tensor for the
403 tensordot operation.
405 Parameters
406 ----------
407 a : array_like
408 Tensor to 'invert'. Its shape must be 'square', i. e.,
409 ``prod(a.shape[:ind]) == prod(a.shape[ind:])``.
410 ind : int, optional
411 Number of first indices that are involved in the inverse sum.
412 Must be a positive integer, default is 2.
414 Returns
415 -------
416 b : ndarray
417 `a`'s tensordot inverse, shape ``a.shape[ind:] + a.shape[:ind]``.
419 Raises
420 ------
421 LinAlgError
422 If `a` is singular or not 'square' (in the above sense).
424 See Also
425 --------
426 numpy.tensordot, tensorsolve
428 Examples
429 --------
430 >>> a = np.eye(4*6)
431 >>> a.shape = (4, 6, 8, 3)
432 >>> ainv = np.linalg.tensorinv(a, ind=2)
433 >>> ainv.shape
434 (8, 3, 4, 6)
435 >>> b = np.random.randn(4, 6)
436 >>> np.allclose(np.tensordot(ainv, b), np.linalg.tensorsolve(a, b))
437 True
439 >>> a = np.eye(4*6)
440 >>> a.shape = (24, 8, 3)
441 >>> ainv = np.linalg.tensorinv(a, ind=1)
442 >>> ainv.shape
443 (8, 3, 24)
444 >>> b = np.random.randn(24)
445 >>> np.allclose(np.tensordot(ainv, b, 1), np.linalg.tensorsolve(a, b))
446 True
448 """
449 a = asarray(a)
450 oldshape = a.shape
451 prod = 1
452 if ind > 0:
453 invshape = oldshape[ind:] + oldshape[:ind]
454 for k in oldshape[ind:]:
455 prod *= k
456 else:
457 raise ValueError("Invalid ind argument.")
458 a = a.reshape(prod, -1)
459 ia = inv(a)
460 return ia.reshape(*invshape)
463# Matrix inversion
465def _unary_dispatcher(a):
466 return (a,)
469@array_function_dispatch(_unary_dispatcher)
470def inv(a):
471 """
472 Compute the (multiplicative) inverse of a matrix.
474 Given a square matrix `a`, return the matrix `ainv` satisfying
475 ``dot(a, ainv) = dot(ainv, a) = eye(a.shape[0])``.
477 Parameters
478 ----------
479 a : (..., M, M) array_like
480 Matrix to be inverted.
482 Returns
483 -------
484 ainv : (..., M, M) ndarray or matrix
485 (Multiplicative) inverse of the matrix `a`.
487 Raises
488 ------
489 LinAlgError
490 If `a` is not square or inversion fails.
492 See Also
493 --------
494 scipy.linalg.inv : Similar function in SciPy.
496 Notes
497 -----
499 .. versionadded:: 1.8.0
501 Broadcasting rules apply, see the `numpy.linalg` documentation for
502 details.
504 Examples
505 --------
506 >>> from numpy.linalg import inv
507 >>> a = np.array([[1., 2.], [3., 4.]])
508 >>> ainv = inv(a)
509 >>> np.allclose(np.dot(a, ainv), np.eye(2))
510 True
511 >>> np.allclose(np.dot(ainv, a), np.eye(2))
512 True
514 If a is a matrix object, then the return value is a matrix as well:
516 >>> ainv = inv(np.matrix(a))
517 >>> ainv
518 matrix([[-2. , 1. ],
519 [ 1.5, -0.5]])
521 Inverses of several matrices can be computed at once:
523 >>> a = np.array([[[1., 2.], [3., 4.]], [[1, 3], [3, 5]]])
524 >>> inv(a)
525 array([[[-2. , 1. ],
526 [ 1.5 , -0.5 ]],
527 [[-1.25, 0.75],
528 [ 0.75, -0.25]]])
530 """
531 a, wrap = _makearray(a)
532 _assert_stacked_2d(a)
533 _assert_stacked_square(a)
534 t, result_t = _commonType(a)
536 signature = 'D->D' if isComplexType(t) else 'd->d'
537 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
538 ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
539 return wrap(ainv.astype(result_t, copy=False))
542def _matrix_power_dispatcher(a, n):
543 return (a,)
546@array_function_dispatch(_matrix_power_dispatcher)
547def matrix_power(a, n):
548 """
549 Raise a square matrix to the (integer) power `n`.
551 For positive integers `n`, the power is computed by repeated matrix
552 squarings and matrix multiplications. If ``n == 0``, the identity matrix
553 of the same shape as M is returned. If ``n < 0``, the inverse
554 is computed and then raised to the ``abs(n)``.
556 .. note:: Stacks of object matrices are not currently supported.
558 Parameters
559 ----------
560 a : (..., M, M) array_like
561 Matrix to be "powered".
562 n : int
563 The exponent can be any integer or long integer, positive,
564 negative, or zero.
566 Returns
567 -------
568 a**n : (..., M, M) ndarray or matrix object
569 The return value is the same shape and type as `M`;
570 if the exponent is positive or zero then the type of the
571 elements is the same as those of `M`. If the exponent is
572 negative the elements are floating-point.
574 Raises
575 ------
576 LinAlgError
577 For matrices that are not square or that (for negative powers) cannot
578 be inverted numerically.
580 Examples
581 --------
582 >>> from numpy.linalg import matrix_power
583 >>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit
584 >>> matrix_power(i, 3) # should = -i
585 array([[ 0, -1],
586 [ 1, 0]])
587 >>> matrix_power(i, 0)
588 array([[1, 0],
589 [0, 1]])
590 >>> matrix_power(i, -3) # should = 1/(-i) = i, but w/ f.p. elements
591 array([[ 0., 1.],
592 [-1., 0.]])
594 Somewhat more sophisticated example
596 >>> q = np.zeros((4, 4))
597 >>> q[0:2, 0:2] = -i
598 >>> q[2:4, 2:4] = i
599 >>> q # one of the three quaternion units not equal to 1
600 array([[ 0., -1., 0., 0.],
601 [ 1., 0., 0., 0.],
602 [ 0., 0., 0., 1.],
603 [ 0., 0., -1., 0.]])
604 >>> matrix_power(q, 2) # = -np.eye(4)
605 array([[-1., 0., 0., 0.],
606 [ 0., -1., 0., 0.],
607 [ 0., 0., -1., 0.],
608 [ 0., 0., 0., -1.]])
610 """
611 a = asanyarray(a)
612 _assert_stacked_2d(a)
613 _assert_stacked_square(a)
615 try:
616 n = operator.index(n)
617 except TypeError as e:
618 raise TypeError("exponent must be an integer") from e
620 # Fall back on dot for object arrays. Object arrays are not supported by
621 # the current implementation of matmul using einsum
622 if a.dtype != object:
623 fmatmul = matmul
624 elif a.ndim == 2:
625 fmatmul = dot
626 else:
627 raise NotImplementedError(
628 "matrix_power not supported for stacks of object arrays")
630 if n == 0:
631 a = empty_like(a)
632 a[...] = eye(a.shape[-2], dtype=a.dtype)
633 return a
635 elif n < 0:
636 a = inv(a)
637 n = abs(n)
639 # short-cuts.
640 if n == 1:
641 return a
643 elif n == 2:
644 return fmatmul(a, a)
646 elif n == 3:
647 return fmatmul(fmatmul(a, a), a)
649 # Use binary decomposition to reduce the number of matrix multiplications.
650 # Here, we iterate over the bits of n, from LSB to MSB, raise `a` to
651 # increasing powers of 2, and multiply into the result as needed.
652 z = result = None
653 while n > 0:
654 z = a if z is None else fmatmul(z, z)
655 n, bit = divmod(n, 2)
656 if bit:
657 result = z if result is None else fmatmul(result, z)
659 return result
662# Cholesky decomposition
665@array_function_dispatch(_unary_dispatcher)
666def cholesky(a):
667 """
668 Cholesky decomposition.
670 Return the Cholesky decomposition, `L * L.H`, of the square matrix `a`,
671 where `L` is lower-triangular and .H is the conjugate transpose operator
672 (which is the ordinary transpose if `a` is real-valued). `a` must be
673 Hermitian (symmetric if real-valued) and positive-definite. No
674 checking is performed to verify whether `a` is Hermitian or not.
675 In addition, only the lower-triangular and diagonal elements of `a`
676 are used. Only `L` is actually returned.
678 Parameters
679 ----------
680 a : (..., M, M) array_like
681 Hermitian (symmetric if all elements are real), positive-definite
682 input matrix.
684 Returns
685 -------
686 L : (..., M, M) array_like
687 Lower-triangular Cholesky factor of `a`. Returns a matrix object if
688 `a` is a matrix object.
690 Raises
691 ------
692 LinAlgError
693 If the decomposition fails, for example, if `a` is not
694 positive-definite.
696 See Also
697 --------
698 scipy.linalg.cholesky : Similar function in SciPy.
699 scipy.linalg.cholesky_banded : Cholesky decompose a banded Hermitian
700 positive-definite matrix.
701 scipy.linalg.cho_factor : Cholesky decomposition of a matrix, to use in
702 `scipy.linalg.cho_solve`.
704 Notes
705 -----
707 .. versionadded:: 1.8.0
709 Broadcasting rules apply, see the `numpy.linalg` documentation for
710 details.
712 The Cholesky decomposition is often used as a fast way of solving
714 .. math:: A \\mathbf{x} = \\mathbf{b}
716 (when `A` is both Hermitian/symmetric and positive-definite).
718 First, we solve for :math:`\\mathbf{y}` in
720 .. math:: L \\mathbf{y} = \\mathbf{b},
722 and then for :math:`\\mathbf{x}` in
724 .. math:: L.H \\mathbf{x} = \\mathbf{y}.
726 Examples
727 --------
728 >>> A = np.array([[1,-2j],[2j,5]])
729 >>> A
730 array([[ 1.+0.j, -0.-2.j],
731 [ 0.+2.j, 5.+0.j]])
732 >>> L = np.linalg.cholesky(A)
733 >>> L
734 array([[1.+0.j, 0.+0.j],
735 [0.+2.j, 1.+0.j]])
736 >>> np.dot(L, L.T.conj()) # verify that L * L.H = A
737 array([[1.+0.j, 0.-2.j],
738 [0.+2.j, 5.+0.j]])
739 >>> A = [[1,-2j],[2j,5]] # what happens if A is only array_like?
740 >>> np.linalg.cholesky(A) # an ndarray object is returned
741 array([[1.+0.j, 0.+0.j],
742 [0.+2.j, 1.+0.j]])
743 >>> # But a matrix object is returned if A is a matrix object
744 >>> np.linalg.cholesky(np.matrix(A))
745 matrix([[ 1.+0.j, 0.+0.j],
746 [ 0.+2.j, 1.+0.j]])
748 """
749 extobj = get_linalg_error_extobj(_raise_linalgerror_nonposdef)
750 gufunc = _umath_linalg.cholesky_lo
751 a, wrap = _makearray(a)
752 _assert_stacked_2d(a)
753 _assert_stacked_square(a)
754 t, result_t = _commonType(a)
755 signature = 'D->D' if isComplexType(t) else 'd->d'
756 r = gufunc(a, signature=signature, extobj=extobj)
757 return wrap(r.astype(result_t, copy=False))
760# QR decomposition
762def _qr_dispatcher(a, mode=None):
763 return (a,)
766@array_function_dispatch(_qr_dispatcher)
767def qr(a, mode='reduced'):
768 """
769 Compute the qr factorization of a matrix.
771 Factor the matrix `a` as *qr*, where `q` is orthonormal and `r` is
772 upper-triangular.
774 Parameters
775 ----------
776 a : array_like, shape (..., M, N)
777 An array-like object with the dimensionality of at least 2.
778 mode : {'reduced', 'complete', 'r', 'raw'}, optional
779 If K = min(M, N), then
781 * 'reduced' : returns q, r with dimensions
782 (..., M, K), (..., K, N) (default)
783 * 'complete' : returns q, r with dimensions (..., M, M), (..., M, N)
784 * 'r' : returns r only with dimensions (..., K, N)
785 * 'raw' : returns h, tau with dimensions (..., N, M), (..., K,)
787 The options 'reduced', 'complete, and 'raw' are new in numpy 1.8,
788 see the notes for more information. The default is 'reduced', and to
789 maintain backward compatibility with earlier versions of numpy both
790 it and the old default 'full' can be omitted. Note that array h
791 returned in 'raw' mode is transposed for calling Fortran. The
792 'economic' mode is deprecated. The modes 'full' and 'economic' may
793 be passed using only the first letter for backwards compatibility,
794 but all others must be spelled out. See the Notes for more
795 explanation.
798 Returns
799 -------
800 q : ndarray of float or complex, optional
801 A matrix with orthonormal columns. When mode = 'complete' the
802 result is an orthogonal/unitary matrix depending on whether or not
803 a is real/complex. The determinant may be either +/- 1 in that
804 case. In case the number of dimensions in the input array is
805 greater than 2 then a stack of the matrices with above properties
806 is returned.
807 r : ndarray of float or complex, optional
808 The upper-triangular matrix or a stack of upper-triangular
809 matrices if the number of dimensions in the input array is greater
810 than 2.
811 (h, tau) : ndarrays of np.double or np.cdouble, optional
812 The array h contains the Householder reflectors that generate q
813 along with r. The tau array contains scaling factors for the
814 reflectors. In the deprecated 'economic' mode only h is returned.
816 Raises
817 ------
818 LinAlgError
819 If factoring fails.
821 See Also
822 --------
823 scipy.linalg.qr : Similar function in SciPy.
824 scipy.linalg.rq : Compute RQ decomposition of a matrix.
826 Notes
827 -----
828 This is an interface to the LAPACK routines ``dgeqrf``, ``zgeqrf``,
829 ``dorgqr``, and ``zungqr``.
831 For more information on the qr factorization, see for example:
832 https://en.wikipedia.org/wiki/QR_factorization
834 Subclasses of `ndarray` are preserved except for the 'raw' mode. So if
835 `a` is of type `matrix`, all the return values will be matrices too.
837 New 'reduced', 'complete', and 'raw' options for mode were added in
838 NumPy 1.8.0 and the old option 'full' was made an alias of 'reduced'. In
839 addition the options 'full' and 'economic' were deprecated. Because
840 'full' was the previous default and 'reduced' is the new default,
841 backward compatibility can be maintained by letting `mode` default.
842 The 'raw' option was added so that LAPACK routines that can multiply
843 arrays by q using the Householder reflectors can be used. Note that in
844 this case the returned arrays are of type np.double or np.cdouble and
845 the h array is transposed to be FORTRAN compatible. No routines using
846 the 'raw' return are currently exposed by numpy, but some are available
847 in lapack_lite and just await the necessary work.
849 Examples
850 --------
851 >>> a = np.random.randn(9, 6)
852 >>> q, r = np.linalg.qr(a)
853 >>> np.allclose(a, np.dot(q, r)) # a does equal qr
854 True
855 >>> r2 = np.linalg.qr(a, mode='r')
856 >>> np.allclose(r, r2) # mode='r' returns the same r as mode='full'
857 True
858 >>> a = np.random.normal(size=(3, 2, 2)) # Stack of 2 x 2 matrices as input
859 >>> q, r = np.linalg.qr(a)
860 >>> q.shape
861 (3, 2, 2)
862 >>> r.shape
863 (3, 2, 2)
864 >>> np.allclose(a, np.matmul(q, r))
865 True
867 Example illustrating a common use of `qr`: solving of least squares
868 problems
870 What are the least-squares-best `m` and `y0` in ``y = y0 + mx`` for
871 the following data: {(0,1), (1,0), (1,2), (2,1)}. (Graph the points
872 and you'll see that it should be y0 = 0, m = 1.) The answer is provided
873 by solving the over-determined matrix equation ``Ax = b``, where::
875 A = array([[0, 1], [1, 1], [1, 1], [2, 1]])
876 x = array([[y0], [m]])
877 b = array([[1], [0], [2], [1]])
879 If A = qr such that q is orthonormal (which is always possible via
880 Gram-Schmidt), then ``x = inv(r) * (q.T) * b``. (In numpy practice,
881 however, we simply use `lstsq`.)
883 >>> A = np.array([[0, 1], [1, 1], [1, 1], [2, 1]])
884 >>> A
885 array([[0, 1],
886 [1, 1],
887 [1, 1],
888 [2, 1]])
889 >>> b = np.array([1, 2, 2, 3])
890 >>> q, r = np.linalg.qr(A)
891 >>> p = np.dot(q.T, b)
892 >>> np.dot(np.linalg.inv(r), p)
893 array([ 1., 1.])
895 """
896 if mode not in ('reduced', 'complete', 'r', 'raw'):
897 if mode in ('f', 'full'):
898 # 2013-04-01, 1.8
899 msg = "".join((
900 "The 'full' option is deprecated in favor of 'reduced'.\n",
901 "For backward compatibility let mode default."))
902 warnings.warn(msg, DeprecationWarning, stacklevel=3)
903 mode = 'reduced'
904 elif mode in ('e', 'economic'):
905 # 2013-04-01, 1.8
906 msg = "The 'economic' option is deprecated."
907 warnings.warn(msg, DeprecationWarning, stacklevel=3)
908 mode = 'economic'
909 else:
910 raise ValueError(f"Unrecognized mode '{mode}'")
912 a, wrap = _makearray(a)
913 _assert_stacked_2d(a)
914 m, n = a.shape[-2:]
915 t, result_t = _commonType(a)
916 a = a.astype(t, copy=True)
917 a = _to_native_byte_order(a)
918 mn = min(m, n)
920 if m <= n:
921 gufunc = _umath_linalg.qr_r_raw_m
922 else:
923 gufunc = _umath_linalg.qr_r_raw_n
925 signature = 'D->D' if isComplexType(t) else 'd->d'
926 extobj = get_linalg_error_extobj(_raise_linalgerror_qr)
927 tau = gufunc(a, signature=signature, extobj=extobj)
929 # handle modes that don't return q
930 if mode == 'r':
931 r = triu(a[..., :mn, :])
932 r = r.astype(result_t, copy=False)
933 return wrap(r)
935 if mode == 'raw':
936 q = transpose(a)
937 q = q.astype(result_t, copy=False)
938 tau = tau.astype(result_t, copy=False)
939 return wrap(q), tau
941 if mode == 'economic':
942 a = a.astype(result_t, copy=False)
943 return wrap(a)
945 # mc is the number of columns in the resulting q
946 # matrix. If the mode is complete then it is
947 # same as number of rows, and if the mode is reduced,
948 # then it is the minimum of number of rows and columns.
949 if mode == 'complete' and m > n:
950 mc = m
951 gufunc = _umath_linalg.qr_complete
952 else:
953 mc = mn
954 gufunc = _umath_linalg.qr_reduced
956 signature = 'DD->D' if isComplexType(t) else 'dd->d'
957 extobj = get_linalg_error_extobj(_raise_linalgerror_qr)
958 q = gufunc(a, tau, signature=signature, extobj=extobj)
959 r = triu(a[..., :mc, :])
961 q = q.astype(result_t, copy=False)
962 r = r.astype(result_t, copy=False)
964 return wrap(q), wrap(r)
966# Eigenvalues
969@array_function_dispatch(_unary_dispatcher)
970def eigvals(a):
971 """
972 Compute the eigenvalues of a general matrix.
974 Main difference between `eigvals` and `eig`: the eigenvectors aren't
975 returned.
977 Parameters
978 ----------
979 a : (..., M, M) array_like
980 A complex- or real-valued matrix whose eigenvalues will be computed.
982 Returns
983 -------
984 w : (..., M,) ndarray
985 The eigenvalues, each repeated according to its multiplicity.
986 They are not necessarily ordered, nor are they necessarily
987 real for real matrices.
989 Raises
990 ------
991 LinAlgError
992 If the eigenvalue computation does not converge.
994 See Also
995 --------
996 eig : eigenvalues and right eigenvectors of general arrays
997 eigvalsh : eigenvalues of real symmetric or complex Hermitian
998 (conjugate symmetric) arrays.
999 eigh : eigenvalues and eigenvectors of real symmetric or complex
1000 Hermitian (conjugate symmetric) arrays.
1001 scipy.linalg.eigvals : Similar function in SciPy.
1003 Notes
1004 -----
1006 .. versionadded:: 1.8.0
1008 Broadcasting rules apply, see the `numpy.linalg` documentation for
1009 details.
1011 This is implemented using the ``_geev`` LAPACK routines which compute
1012 the eigenvalues and eigenvectors of general square arrays.
1014 Examples
1015 --------
1016 Illustration, using the fact that the eigenvalues of a diagonal matrix
1017 are its diagonal elements, that multiplying a matrix on the left
1018 by an orthogonal matrix, `Q`, and on the right by `Q.T` (the transpose
1019 of `Q`), preserves the eigenvalues of the "middle" matrix. In other words,
1020 if `Q` is orthogonal, then ``Q * A * Q.T`` has the same eigenvalues as
1021 ``A``:
1023 >>> from numpy import linalg as LA
1024 >>> x = np.random.random()
1025 >>> Q = np.array([[np.cos(x), -np.sin(x)], [np.sin(x), np.cos(x)]])
1026 >>> LA.norm(Q[0, :]), LA.norm(Q[1, :]), np.dot(Q[0, :],Q[1, :])
1027 (1.0, 1.0, 0.0)
1029 Now multiply a diagonal matrix by ``Q`` on one side and by ``Q.T`` on the other:
1031 >>> D = np.diag((-1,1))
1032 >>> LA.eigvals(D)
1033 array([-1., 1.])
1034 >>> A = np.dot(Q, D)
1035 >>> A = np.dot(A, Q.T)
1036 >>> LA.eigvals(A)
1037 array([ 1., -1.]) # random
1039 """
1040 a, wrap = _makearray(a)
1041 _assert_stacked_2d(a)
1042 _assert_stacked_square(a)
1043 _assert_finite(a)
1044 t, result_t = _commonType(a)
1046 extobj = get_linalg_error_extobj(
1047 _raise_linalgerror_eigenvalues_nonconvergence)
1048 signature = 'D->D' if isComplexType(t) else 'd->D'
1049 w = _umath_linalg.eigvals(a, signature=signature, extobj=extobj)
1051 if not isComplexType(t):
1052 if all(w.imag == 0):
1053 w = w.real
1054 result_t = _realType(result_t)
1055 else:
1056 result_t = _complexType(result_t)
1058 return w.astype(result_t, copy=False)
1061def _eigvalsh_dispatcher(a, UPLO=None):
1062 return (a,)
1065@array_function_dispatch(_eigvalsh_dispatcher)
1066def eigvalsh(a, UPLO='L'):
1067 """
1068 Compute the eigenvalues of a complex Hermitian or real symmetric matrix.
1070 Main difference from eigh: the eigenvectors are not computed.
1072 Parameters
1073 ----------
1074 a : (..., M, M) array_like
1075 A complex- or real-valued matrix whose eigenvalues are to be
1076 computed.
1077 UPLO : {'L', 'U'}, optional
1078 Specifies whether the calculation is done with the lower triangular
1079 part of `a` ('L', default) or the upper triangular part ('U').
1080 Irrespective of this value only the real parts of the diagonal will
1081 be considered in the computation to preserve the notion of a Hermitian
1082 matrix. It therefore follows that the imaginary part of the diagonal
1083 will always be treated as zero.
1085 Returns
1086 -------
1087 w : (..., M,) ndarray
1088 The eigenvalues in ascending order, each repeated according to
1089 its multiplicity.
1091 Raises
1092 ------
1093 LinAlgError
1094 If the eigenvalue computation does not converge.
1096 See Also
1097 --------
1098 eigh : eigenvalues and eigenvectors of real symmetric or complex Hermitian
1099 (conjugate symmetric) arrays.
1100 eigvals : eigenvalues of general real or complex arrays.
1101 eig : eigenvalues and right eigenvectors of general real or complex
1102 arrays.
1103 scipy.linalg.eigvalsh : Similar function in SciPy.
1105 Notes
1106 -----
1108 .. versionadded:: 1.8.0
1110 Broadcasting rules apply, see the `numpy.linalg` documentation for
1111 details.
1113 The eigenvalues are computed using LAPACK routines ``_syevd``, ``_heevd``.
1115 Examples
1116 --------
1117 >>> from numpy import linalg as LA
1118 >>> a = np.array([[1, -2j], [2j, 5]])
1119 >>> LA.eigvalsh(a)
1120 array([ 0.17157288, 5.82842712]) # may vary
1122 >>> # demonstrate the treatment of the imaginary part of the diagonal
1123 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1124 >>> a
1125 array([[5.+2.j, 9.-2.j],
1126 [0.+2.j, 2.-1.j]])
1127 >>> # with UPLO='L' this is numerically equivalent to using LA.eigvals()
1128 >>> # with:
1129 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1130 >>> b
1131 array([[5.+0.j, 0.-2.j],
1132 [0.+2.j, 2.+0.j]])
1133 >>> wa = LA.eigvalsh(a)
1134 >>> wb = LA.eigvals(b)
1135 >>> wa; wb
1136 array([1., 6.])
1137 array([6.+0.j, 1.+0.j])
1139 """
1140 UPLO = UPLO.upper()
1141 if UPLO not in ('L', 'U'):
1142 raise ValueError("UPLO argument must be 'L' or 'U'")
1144 extobj = get_linalg_error_extobj(
1145 _raise_linalgerror_eigenvalues_nonconvergence)
1146 if UPLO == 'L':
1147 gufunc = _umath_linalg.eigvalsh_lo
1148 else:
1149 gufunc = _umath_linalg.eigvalsh_up
1151 a, wrap = _makearray(a)
1152 _assert_stacked_2d(a)
1153 _assert_stacked_square(a)
1154 t, result_t = _commonType(a)
1155 signature = 'D->d' if isComplexType(t) else 'd->d'
1156 w = gufunc(a, signature=signature, extobj=extobj)
1157 return w.astype(_realType(result_t), copy=False)
1159def _convertarray(a):
1160 t, result_t = _commonType(a)
1161 a = a.astype(t).T.copy()
1162 return a, t, result_t
1165# Eigenvectors
1168@array_function_dispatch(_unary_dispatcher)
1169def eig(a):
1170 """
1171 Compute the eigenvalues and right eigenvectors of a square array.
1173 Parameters
1174 ----------
1175 a : (..., M, M) array
1176 Matrices for which the eigenvalues and right eigenvectors will
1177 be computed
1179 Returns
1180 -------
1181 w : (..., M) array
1182 The eigenvalues, each repeated according to its multiplicity.
1183 The eigenvalues are not necessarily ordered. The resulting
1184 array will be of complex type, unless the imaginary part is
1185 zero in which case it will be cast to a real type. When `a`
1186 is real the resulting eigenvalues will be real (0 imaginary
1187 part) or occur in conjugate pairs
1189 v : (..., M, M) array
1190 The normalized (unit "length") eigenvectors, such that the
1191 column ``v[:,i]`` is the eigenvector corresponding to the
1192 eigenvalue ``w[i]``.
1194 Raises
1195 ------
1196 LinAlgError
1197 If the eigenvalue computation does not converge.
1199 See Also
1200 --------
1201 eigvals : eigenvalues of a non-symmetric array.
1202 eigh : eigenvalues and eigenvectors of a real symmetric or complex
1203 Hermitian (conjugate symmetric) array.
1204 eigvalsh : eigenvalues of a real symmetric or complex Hermitian
1205 (conjugate symmetric) array.
1206 scipy.linalg.eig : Similar function in SciPy that also solves the
1207 generalized eigenvalue problem.
1208 scipy.linalg.schur : Best choice for unitary and other non-Hermitian
1209 normal matrices.
1211 Notes
1212 -----
1214 .. versionadded:: 1.8.0
1216 Broadcasting rules apply, see the `numpy.linalg` documentation for
1217 details.
1219 This is implemented using the ``_geev`` LAPACK routines which compute
1220 the eigenvalues and eigenvectors of general square arrays.
1222 The number `w` is an eigenvalue of `a` if there exists a vector
1223 `v` such that ``a @ v = w * v``. Thus, the arrays `a`, `w`, and
1224 `v` satisfy the equations ``a @ v[:,i] = w[i] * v[:,i]``
1225 for :math:`i \\in \\{0,...,M-1\\}`.
1227 The array `v` of eigenvectors may not be of maximum rank, that is, some
1228 of the columns may be linearly dependent, although round-off error may
1229 obscure that fact. If the eigenvalues are all different, then theoretically
1230 the eigenvectors are linearly independent and `a` can be diagonalized by
1231 a similarity transformation using `v`, i.e, ``inv(v) @ a @ v`` is diagonal.
1233 For non-Hermitian normal matrices the SciPy function `scipy.linalg.schur`
1234 is preferred because the matrix `v` is guaranteed to be unitary, which is
1235 not the case when using `eig`. The Schur factorization produces an
1236 upper triangular matrix rather than a diagonal matrix, but for normal
1237 matrices only the diagonal of the upper triangular matrix is needed, the
1238 rest is roundoff error.
1240 Finally, it is emphasized that `v` consists of the *right* (as in
1241 right-hand side) eigenvectors of `a`. A vector `y` satisfying
1242 ``y.T @ a = z * y.T`` for some number `z` is called a *left*
1243 eigenvector of `a`, and, in general, the left and right eigenvectors
1244 of a matrix are not necessarily the (perhaps conjugate) transposes
1245 of each other.
1247 References
1248 ----------
1249 G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, FL,
1250 Academic Press, Inc., 1980, Various pp.
1252 Examples
1253 --------
1254 >>> from numpy import linalg as LA
1256 (Almost) trivial example with real e-values and e-vectors.
1258 >>> w, v = LA.eig(np.diag((1, 2, 3)))
1259 >>> w; v
1260 array([1., 2., 3.])
1261 array([[1., 0., 0.],
1262 [0., 1., 0.],
1263 [0., 0., 1.]])
1265 Real matrix possessing complex e-values and e-vectors; note that the
1266 e-values are complex conjugates of each other.
1268 >>> w, v = LA.eig(np.array([[1, -1], [1, 1]]))
1269 >>> w; v
1270 array([1.+1.j, 1.-1.j])
1271 array([[0.70710678+0.j , 0.70710678-0.j ],
1272 [0. -0.70710678j, 0. +0.70710678j]])
1274 Complex-valued matrix with real e-values (but complex-valued e-vectors);
1275 note that ``a.conj().T == a``, i.e., `a` is Hermitian.
1277 >>> a = np.array([[1, 1j], [-1j, 1]])
1278 >>> w, v = LA.eig(a)
1279 >>> w; v
1280 array([2.+0.j, 0.+0.j])
1281 array([[ 0. +0.70710678j, 0.70710678+0.j ], # may vary
1282 [ 0.70710678+0.j , -0. +0.70710678j]])
1284 Be careful about round-off error!
1286 >>> a = np.array([[1 + 1e-9, 0], [0, 1 - 1e-9]])
1287 >>> # Theor. e-values are 1 +/- 1e-9
1288 >>> w, v = LA.eig(a)
1289 >>> w; v
1290 array([1., 1.])
1291 array([[1., 0.],
1292 [0., 1.]])
1294 """
1295 a, wrap = _makearray(a)
1296 _assert_stacked_2d(a)
1297 _assert_stacked_square(a)
1298 _assert_finite(a)
1299 t, result_t = _commonType(a)
1301 extobj = get_linalg_error_extobj(
1302 _raise_linalgerror_eigenvalues_nonconvergence)
1303 signature = 'D->DD' if isComplexType(t) else 'd->DD'
1304 w, vt = _umath_linalg.eig(a, signature=signature, extobj=extobj)
1306 if not isComplexType(t) and all(w.imag == 0.0):
1307 w = w.real
1308 vt = vt.real
1309 result_t = _realType(result_t)
1310 else:
1311 result_t = _complexType(result_t)
1313 vt = vt.astype(result_t, copy=False)
1314 return w.astype(result_t, copy=False), wrap(vt)
1317@array_function_dispatch(_eigvalsh_dispatcher)
1318def eigh(a, UPLO='L'):
1319 """
1320 Return the eigenvalues and eigenvectors of a complex Hermitian
1321 (conjugate symmetric) or a real symmetric matrix.
1323 Returns two objects, a 1-D array containing the eigenvalues of `a`, and
1324 a 2-D square array or matrix (depending on the input type) of the
1325 corresponding eigenvectors (in columns).
1327 Parameters
1328 ----------
1329 a : (..., M, M) array
1330 Hermitian or real symmetric matrices whose eigenvalues and
1331 eigenvectors are to be computed.
1332 UPLO : {'L', 'U'}, optional
1333 Specifies whether the calculation is done with the lower triangular
1334 part of `a` ('L', default) or the upper triangular part ('U').
1335 Irrespective of this value only the real parts of the diagonal will
1336 be considered in the computation to preserve the notion of a Hermitian
1337 matrix. It therefore follows that the imaginary part of the diagonal
1338 will always be treated as zero.
1340 Returns
1341 -------
1342 w : (..., M) ndarray
1343 The eigenvalues in ascending order, each repeated according to
1344 its multiplicity.
1345 v : {(..., M, M) ndarray, (..., M, M) matrix}
1346 The column ``v[:, i]`` is the normalized eigenvector corresponding
1347 to the eigenvalue ``w[i]``. Will return a matrix object if `a` is
1348 a matrix object.
1350 Raises
1351 ------
1352 LinAlgError
1353 If the eigenvalue computation does not converge.
1355 See Also
1356 --------
1357 eigvalsh : eigenvalues of real symmetric or complex Hermitian
1358 (conjugate symmetric) arrays.
1359 eig : eigenvalues and right eigenvectors for non-symmetric arrays.
1360 eigvals : eigenvalues of non-symmetric arrays.
1361 scipy.linalg.eigh : Similar function in SciPy (but also solves the
1362 generalized eigenvalue problem).
1364 Notes
1365 -----
1367 .. versionadded:: 1.8.0
1369 Broadcasting rules apply, see the `numpy.linalg` documentation for
1370 details.
1372 The eigenvalues/eigenvectors are computed using LAPACK routines ``_syevd``,
1373 ``_heevd``.
1375 The eigenvalues of real symmetric or complex Hermitian matrices are
1376 always real. [1]_ The array `v` of (column) eigenvectors is unitary
1377 and `a`, `w`, and `v` satisfy the equations
1378 ``dot(a, v[:, i]) = w[i] * v[:, i]``.
1380 References
1381 ----------
1382 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
1383 FL, Academic Press, Inc., 1980, pg. 222.
1385 Examples
1386 --------
1387 >>> from numpy import linalg as LA
1388 >>> a = np.array([[1, -2j], [2j, 5]])
1389 >>> a
1390 array([[ 1.+0.j, -0.-2.j],
1391 [ 0.+2.j, 5.+0.j]])
1392 >>> w, v = LA.eigh(a)
1393 >>> w; v
1394 array([0.17157288, 5.82842712])
1395 array([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1396 [ 0. +0.38268343j, 0. -0.92387953j]])
1398 >>> np.dot(a, v[:, 0]) - w[0] * v[:, 0] # verify 1st e-val/vec pair
1399 array([5.55111512e-17+0.0000000e+00j, 0.00000000e+00+1.2490009e-16j])
1400 >>> np.dot(a, v[:, 1]) - w[1] * v[:, 1] # verify 2nd e-val/vec pair
1401 array([0.+0.j, 0.+0.j])
1403 >>> A = np.matrix(a) # what happens if input is a matrix object
1404 >>> A
1405 matrix([[ 1.+0.j, -0.-2.j],
1406 [ 0.+2.j, 5.+0.j]])
1407 >>> w, v = LA.eigh(A)
1408 >>> w; v
1409 array([0.17157288, 5.82842712])
1410 matrix([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1411 [ 0. +0.38268343j, 0. -0.92387953j]])
1413 >>> # demonstrate the treatment of the imaginary part of the diagonal
1414 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1415 >>> a
1416 array([[5.+2.j, 9.-2.j],
1417 [0.+2.j, 2.-1.j]])
1418 >>> # with UPLO='L' this is numerically equivalent to using LA.eig() with:
1419 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1420 >>> b
1421 array([[5.+0.j, 0.-2.j],
1422 [0.+2.j, 2.+0.j]])
1423 >>> wa, va = LA.eigh(a)
1424 >>> wb, vb = LA.eig(b)
1425 >>> wa; wb
1426 array([1., 6.])
1427 array([6.+0.j, 1.+0.j])
1428 >>> va; vb
1429 array([[-0.4472136 +0.j , -0.89442719+0.j ], # may vary
1430 [ 0. +0.89442719j, 0. -0.4472136j ]])
1431 array([[ 0.89442719+0.j , -0. +0.4472136j],
1432 [-0. +0.4472136j, 0.89442719+0.j ]])
1433 """
1434 UPLO = UPLO.upper()
1435 if UPLO not in ('L', 'U'):
1436 raise ValueError("UPLO argument must be 'L' or 'U'")
1438 a, wrap = _makearray(a)
1439 _assert_stacked_2d(a)
1440 _assert_stacked_square(a)
1441 t, result_t = _commonType(a)
1443 extobj = get_linalg_error_extobj(
1444 _raise_linalgerror_eigenvalues_nonconvergence)
1445 if UPLO == 'L':
1446 gufunc = _umath_linalg.eigh_lo
1447 else:
1448 gufunc = _umath_linalg.eigh_up
1450 signature = 'D->dD' if isComplexType(t) else 'd->dd'
1451 w, vt = gufunc(a, signature=signature, extobj=extobj)
1452 w = w.astype(_realType(result_t), copy=False)
1453 vt = vt.astype(result_t, copy=False)
1454 return w, wrap(vt)
1457# Singular value decomposition
1459def _svd_dispatcher(a, full_matrices=None, compute_uv=None, hermitian=None):
1460 return (a,)
1463@array_function_dispatch(_svd_dispatcher)
1464def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
1465 """
1466 Singular Value Decomposition.
1468 When `a` is a 2D array, and ``full_matrices=False``, then it is
1469 factorized as ``u @ np.diag(s) @ vh = (u * s) @ vh``, where
1470 `u` and the Hermitian transpose of `vh` are 2D arrays with
1471 orthonormal columns and `s` is a 1D array of `a`'s singular
1472 values. When `a` is higher-dimensional, SVD is applied in
1473 stacked mode as explained below.
1475 Parameters
1476 ----------
1477 a : (..., M, N) array_like
1478 A real or complex array with ``a.ndim >= 2``.
1479 full_matrices : bool, optional
1480 If True (default), `u` and `vh` have the shapes ``(..., M, M)`` and
1481 ``(..., N, N)``, respectively. Otherwise, the shapes are
1482 ``(..., M, K)`` and ``(..., K, N)``, respectively, where
1483 ``K = min(M, N)``.
1484 compute_uv : bool, optional
1485 Whether or not to compute `u` and `vh` in addition to `s`. True
1486 by default.
1487 hermitian : bool, optional
1488 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
1489 enabling a more efficient method for finding singular values.
1490 Defaults to False.
1492 .. versionadded:: 1.17.0
1494 Returns
1495 -------
1496 u : { (..., M, M), (..., M, K) } array
1497 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1498 size as those of the input `a`. The size of the last two dimensions
1499 depends on the value of `full_matrices`. Only returned when
1500 `compute_uv` is True.
1501 s : (..., K) array
1502 Vector(s) with the singular values, within each vector sorted in
1503 descending order. The first ``a.ndim - 2`` dimensions have the same
1504 size as those of the input `a`.
1505 vh : { (..., N, N), (..., K, N) } array
1506 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1507 size as those of the input `a`. The size of the last two dimensions
1508 depends on the value of `full_matrices`. Only returned when
1509 `compute_uv` is True.
1511 Raises
1512 ------
1513 LinAlgError
1514 If SVD computation does not converge.
1516 See Also
1517 --------
1518 scipy.linalg.svd : Similar function in SciPy.
1519 scipy.linalg.svdvals : Compute singular values of a matrix.
1521 Notes
1522 -----
1524 .. versionchanged:: 1.8.0
1525 Broadcasting rules apply, see the `numpy.linalg` documentation for
1526 details.
1528 The decomposition is performed using LAPACK routine ``_gesdd``.
1530 SVD is usually described for the factorization of a 2D matrix :math:`A`.
1531 The higher-dimensional case will be discussed below. In the 2D case, SVD is
1532 written as :math:`A = U S V^H`, where :math:`A = a`, :math:`U= u`,
1533 :math:`S= \\mathtt{np.diag}(s)` and :math:`V^H = vh`. The 1D array `s`
1534 contains the singular values of `a` and `u` and `vh` are unitary. The rows
1535 of `vh` are the eigenvectors of :math:`A^H A` and the columns of `u` are
1536 the eigenvectors of :math:`A A^H`. In both cases the corresponding
1537 (possibly non-zero) eigenvalues are given by ``s**2``.
1539 If `a` has more than two dimensions, then broadcasting rules apply, as
1540 explained in :ref:`routines.linalg-broadcasting`. This means that SVD is
1541 working in "stacked" mode: it iterates over all indices of the first
1542 ``a.ndim - 2`` dimensions and for each combination SVD is applied to the
1543 last two indices. The matrix `a` can be reconstructed from the
1544 decomposition with either ``(u * s[..., None, :]) @ vh`` or
1545 ``u @ (s[..., None] * vh)``. (The ``@`` operator can be replaced by the
1546 function ``np.matmul`` for python versions below 3.5.)
1548 If `a` is a ``matrix`` object (as opposed to an ``ndarray``), then so are
1549 all the return values.
1551 Examples
1552 --------
1553 >>> a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6)
1554 >>> b = np.random.randn(2, 7, 8, 3) + 1j*np.random.randn(2, 7, 8, 3)
1556 Reconstruction based on full SVD, 2D case:
1558 >>> u, s, vh = np.linalg.svd(a, full_matrices=True)
1559 >>> u.shape, s.shape, vh.shape
1560 ((9, 9), (6,), (6, 6))
1561 >>> np.allclose(a, np.dot(u[:, :6] * s, vh))
1562 True
1563 >>> smat = np.zeros((9, 6), dtype=complex)
1564 >>> smat[:6, :6] = np.diag(s)
1565 >>> np.allclose(a, np.dot(u, np.dot(smat, vh)))
1566 True
1568 Reconstruction based on reduced SVD, 2D case:
1570 >>> u, s, vh = np.linalg.svd(a, full_matrices=False)
1571 >>> u.shape, s.shape, vh.shape
1572 ((9, 6), (6,), (6, 6))
1573 >>> np.allclose(a, np.dot(u * s, vh))
1574 True
1575 >>> smat = np.diag(s)
1576 >>> np.allclose(a, np.dot(u, np.dot(smat, vh)))
1577 True
1579 Reconstruction based on full SVD, 4D case:
1581 >>> u, s, vh = np.linalg.svd(b, full_matrices=True)
1582 >>> u.shape, s.shape, vh.shape
1583 ((2, 7, 8, 8), (2, 7, 3), (2, 7, 3, 3))
1584 >>> np.allclose(b, np.matmul(u[..., :3] * s[..., None, :], vh))
1585 True
1586 >>> np.allclose(b, np.matmul(u[..., :3], s[..., None] * vh))
1587 True
1589 Reconstruction based on reduced SVD, 4D case:
1591 >>> u, s, vh = np.linalg.svd(b, full_matrices=False)
1592 >>> u.shape, s.shape, vh.shape
1593 ((2, 7, 8, 3), (2, 7, 3), (2, 7, 3, 3))
1594 >>> np.allclose(b, np.matmul(u * s[..., None, :], vh))
1595 True
1596 >>> np.allclose(b, np.matmul(u, s[..., None] * vh))
1597 True
1599 """
1600 import numpy as _nx
1601 a, wrap = _makearray(a)
1603 if hermitian:
1604 # note: lapack svd returns eigenvalues with s ** 2 sorted descending,
1605 # but eig returns s sorted ascending, so we re-order the eigenvalues
1606 # and related arrays to have the correct order
1607 if compute_uv:
1608 s, u = eigh(a)
1609 sgn = sign(s)
1610 s = abs(s)
1611 sidx = argsort(s)[..., ::-1]
1612 sgn = _nx.take_along_axis(sgn, sidx, axis=-1)
1613 s = _nx.take_along_axis(s, sidx, axis=-1)
1614 u = _nx.take_along_axis(u, sidx[..., None, :], axis=-1)
1615 # singular values are unsigned, move the sign into v
1616 vt = transpose(u * sgn[..., None, :]).conjugate()
1617 return wrap(u), s, wrap(vt)
1618 else:
1619 s = eigvalsh(a)
1620 s = abs(s)
1621 return sort(s)[..., ::-1]
1623 _assert_stacked_2d(a)
1624 t, result_t = _commonType(a)
1626 extobj = get_linalg_error_extobj(_raise_linalgerror_svd_nonconvergence)
1628 m, n = a.shape[-2:]
1629 if compute_uv:
1630 if full_matrices:
1631 if m < n:
1632 gufunc = _umath_linalg.svd_m_f
1633 else:
1634 gufunc = _umath_linalg.svd_n_f
1635 else:
1636 if m < n:
1637 gufunc = _umath_linalg.svd_m_s
1638 else:
1639 gufunc = _umath_linalg.svd_n_s
1641 signature = 'D->DdD' if isComplexType(t) else 'd->ddd'
1642 u, s, vh = gufunc(a, signature=signature, extobj=extobj)
1643 u = u.astype(result_t, copy=False)
1644 s = s.astype(_realType(result_t), copy=False)
1645 vh = vh.astype(result_t, copy=False)
1646 return wrap(u), s, wrap(vh)
1647 else:
1648 if m < n:
1649 gufunc = _umath_linalg.svd_m
1650 else:
1651 gufunc = _umath_linalg.svd_n
1653 signature = 'D->d' if isComplexType(t) else 'd->d'
1654 s = gufunc(a, signature=signature, extobj=extobj)
1655 s = s.astype(_realType(result_t), copy=False)
1656 return s
1659def _cond_dispatcher(x, p=None):
1660 return (x,)
1663@array_function_dispatch(_cond_dispatcher)
1664def cond(x, p=None):
1665 """
1666 Compute the condition number of a matrix.
1668 This function is capable of returning the condition number using
1669 one of seven different norms, depending on the value of `p` (see
1670 Parameters below).
1672 Parameters
1673 ----------
1674 x : (..., M, N) array_like
1675 The matrix whose condition number is sought.
1676 p : {None, 1, -1, 2, -2, inf, -inf, 'fro'}, optional
1677 Order of the norm used in the condition number computation:
1679 ===== ============================
1680 p norm for matrices
1681 ===== ============================
1682 None 2-norm, computed directly using the ``SVD``
1683 'fro' Frobenius norm
1684 inf max(sum(abs(x), axis=1))
1685 -inf min(sum(abs(x), axis=1))
1686 1 max(sum(abs(x), axis=0))
1687 -1 min(sum(abs(x), axis=0))
1688 2 2-norm (largest sing. value)
1689 -2 smallest singular value
1690 ===== ============================
1692 inf means the `numpy.inf` object, and the Frobenius norm is
1693 the root-of-sum-of-squares norm.
1695 Returns
1696 -------
1697 c : {float, inf}
1698 The condition number of the matrix. May be infinite.
1700 See Also
1701 --------
1702 numpy.linalg.norm
1704 Notes
1705 -----
1706 The condition number of `x` is defined as the norm of `x` times the
1707 norm of the inverse of `x` [1]_; the norm can be the usual L2-norm
1708 (root-of-sum-of-squares) or one of a number of other matrix norms.
1710 References
1711 ----------
1712 .. [1] G. Strang, *Linear Algebra and Its Applications*, Orlando, FL,
1713 Academic Press, Inc., 1980, pg. 285.
1715 Examples
1716 --------
1717 >>> from numpy import linalg as LA
1718 >>> a = np.array([[1, 0, -1], [0, 1, 0], [1, 0, 1]])
1719 >>> a
1720 array([[ 1, 0, -1],
1721 [ 0, 1, 0],
1722 [ 1, 0, 1]])
1723 >>> LA.cond(a)
1724 1.4142135623730951
1725 >>> LA.cond(a, 'fro')
1726 3.1622776601683795
1727 >>> LA.cond(a, np.inf)
1728 2.0
1729 >>> LA.cond(a, -np.inf)
1730 1.0
1731 >>> LA.cond(a, 1)
1732 2.0
1733 >>> LA.cond(a, -1)
1734 1.0
1735 >>> LA.cond(a, 2)
1736 1.4142135623730951
1737 >>> LA.cond(a, -2)
1738 0.70710678118654746 # may vary
1739 >>> min(LA.svd(a, compute_uv=False))*min(LA.svd(LA.inv(a), compute_uv=False))
1740 0.70710678118654746 # may vary
1742 """
1743 x = asarray(x) # in case we have a matrix
1744 if _is_empty_2d(x):
1745 raise LinAlgError("cond is not defined on empty arrays")
1746 if p is None or p == 2 or p == -2:
1747 s = svd(x, compute_uv=False)
1748 with errstate(all='ignore'):
1749 if p == -2:
1750 r = s[..., -1] / s[..., 0]
1751 else:
1752 r = s[..., 0] / s[..., -1]
1753 else:
1754 # Call inv(x) ignoring errors. The result array will
1755 # contain nans in the entries where inversion failed.
1756 _assert_stacked_2d(x)
1757 _assert_stacked_square(x)
1758 t, result_t = _commonType(x)
1759 signature = 'D->D' if isComplexType(t) else 'd->d'
1760 with errstate(all='ignore'):
1761 invx = _umath_linalg.inv(x, signature=signature)
1762 r = norm(x, p, axis=(-2, -1)) * norm(invx, p, axis=(-2, -1))
1763 r = r.astype(result_t, copy=False)
1765 # Convert nans to infs unless the original array had nan entries
1766 r = asarray(r)
1767 nan_mask = isnan(r)
1768 if nan_mask.any():
1769 nan_mask &= ~isnan(x).any(axis=(-2, -1))
1770 if r.ndim > 0:
1771 r[nan_mask] = Inf
1772 elif nan_mask:
1773 r[()] = Inf
1775 # Convention is to return scalars instead of 0d arrays
1776 if r.ndim == 0:
1777 r = r[()]
1779 return r
1782def _matrix_rank_dispatcher(A, tol=None, hermitian=None):
1783 return (A,)
1786@array_function_dispatch(_matrix_rank_dispatcher)
1787def matrix_rank(A, tol=None, hermitian=False):
1788 """
1789 Return matrix rank of array using SVD method
1791 Rank of the array is the number of singular values of the array that are
1792 greater than `tol`.
1794 .. versionchanged:: 1.14
1795 Can now operate on stacks of matrices
1797 Parameters
1798 ----------
1799 A : {(M,), (..., M, N)} array_like
1800 Input vector or stack of matrices.
1801 tol : (...) array_like, float, optional
1802 Threshold below which SVD values are considered zero. If `tol` is
1803 None, and ``S`` is an array with singular values for `M`, and
1804 ``eps`` is the epsilon value for datatype of ``S``, then `tol` is
1805 set to ``S.max() * max(M, N) * eps``.
1807 .. versionchanged:: 1.14
1808 Broadcasted against the stack of matrices
1809 hermitian : bool, optional
1810 If True, `A` is assumed to be Hermitian (symmetric if real-valued),
1811 enabling a more efficient method for finding singular values.
1812 Defaults to False.
1814 .. versionadded:: 1.14
1816 Returns
1817 -------
1818 rank : (...) array_like
1819 Rank of A.
1821 Notes
1822 -----
1823 The default threshold to detect rank deficiency is a test on the magnitude
1824 of the singular values of `A`. By default, we identify singular values less
1825 than ``S.max() * max(M, N) * eps`` as indicating rank deficiency (with
1826 the symbols defined above). This is the algorithm MATLAB uses [1]. It also
1827 appears in *Numerical recipes* in the discussion of SVD solutions for linear
1828 least squares [2].
1830 This default threshold is designed to detect rank deficiency accounting for
1831 the numerical errors of the SVD computation. Imagine that there is a column
1832 in `A` that is an exact (in floating point) linear combination of other
1833 columns in `A`. Computing the SVD on `A` will not produce a singular value
1834 exactly equal to 0 in general: any difference of the smallest SVD value from
1835 0 will be caused by numerical imprecision in the calculation of the SVD.
1836 Our threshold for small SVD values takes this numerical imprecision into
1837 account, and the default threshold will detect such numerical rank
1838 deficiency. The threshold may declare a matrix `A` rank deficient even if
1839 the linear combination of some columns of `A` is not exactly equal to
1840 another column of `A` but only numerically very close to another column of
1841 `A`.
1843 We chose our default threshold because it is in wide use. Other thresholds
1844 are possible. For example, elsewhere in the 2007 edition of *Numerical
1845 recipes* there is an alternative threshold of ``S.max() *
1846 np.finfo(A.dtype).eps / 2. * np.sqrt(m + n + 1.)``. The authors describe
1847 this threshold as being based on "expected roundoff error" (p 71).
1849 The thresholds above deal with floating point roundoff error in the
1850 calculation of the SVD. However, you may have more information about the
1851 sources of error in `A` that would make you consider other tolerance values
1852 to detect *effective* rank deficiency. The most useful measure of the
1853 tolerance depends on the operations you intend to use on your matrix. For
1854 example, if your data come from uncertain measurements with uncertainties
1855 greater than floating point epsilon, choosing a tolerance near that
1856 uncertainty may be preferable. The tolerance may be absolute if the
1857 uncertainties are absolute rather than relative.
1859 References
1860 ----------
1861 .. [1] MATLAB reference documentation, "Rank"
1862 https://www.mathworks.com/help/techdoc/ref/rank.html
1863 .. [2] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery,
1864 "Numerical Recipes (3rd edition)", Cambridge University Press, 2007,
1865 page 795.
1867 Examples
1868 --------
1869 >>> from numpy.linalg import matrix_rank
1870 >>> matrix_rank(np.eye(4)) # Full rank matrix
1871 4
1872 >>> I=np.eye(4); I[-1,-1] = 0. # rank deficient matrix
1873 >>> matrix_rank(I)
1874 3
1875 >>> matrix_rank(np.ones((4,))) # 1 dimension - rank 1 unless all 0
1876 1
1877 >>> matrix_rank(np.zeros((4,)))
1878 0
1879 """
1880 A = asarray(A)
1881 if A.ndim < 2:
1882 return int(not all(A==0))
1883 S = svd(A, compute_uv=False, hermitian=hermitian)
1884 if tol is None:
1885 tol = S.max(axis=-1, keepdims=True) * max(A.shape[-2:]) * finfo(S.dtype).eps
1886 else:
1887 tol = asarray(tol)[..., newaxis]
1888 return count_nonzero(S > tol, axis=-1)
1891# Generalized inverse
1893def _pinv_dispatcher(a, rcond=None, hermitian=None):
1894 return (a,)
1897@array_function_dispatch(_pinv_dispatcher)
1898def pinv(a, rcond=1e-15, hermitian=False):
1899 """
1900 Compute the (Moore-Penrose) pseudo-inverse of a matrix.
1902 Calculate the generalized inverse of a matrix using its
1903 singular-value decomposition (SVD) and including all
1904 *large* singular values.
1906 .. versionchanged:: 1.14
1907 Can now operate on stacks of matrices
1909 Parameters
1910 ----------
1911 a : (..., M, N) array_like
1912 Matrix or stack of matrices to be pseudo-inverted.
1913 rcond : (...) array_like of float
1914 Cutoff for small singular values.
1915 Singular values less than or equal to
1916 ``rcond * largest_singular_value`` are set to zero.
1917 Broadcasts against the stack of matrices.
1918 hermitian : bool, optional
1919 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
1920 enabling a more efficient method for finding singular values.
1921 Defaults to False.
1923 .. versionadded:: 1.17.0
1925 Returns
1926 -------
1927 B : (..., N, M) ndarray
1928 The pseudo-inverse of `a`. If `a` is a `matrix` instance, then so
1929 is `B`.
1931 Raises
1932 ------
1933 LinAlgError
1934 If the SVD computation does not converge.
1936 See Also
1937 --------
1938 scipy.linalg.pinv : Similar function in SciPy.
1939 scipy.linalg.pinvh : Compute the (Moore-Penrose) pseudo-inverse of a
1940 Hermitian matrix.
1942 Notes
1943 -----
1944 The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
1945 defined as: "the matrix that 'solves' [the least-squares problem]
1946 :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then
1947 :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`.
1949 It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular
1950 value decomposition of A, then
1951 :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are
1952 orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting
1953 of A's so-called singular values, (followed, typically, by
1954 zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix
1955 consisting of the reciprocals of A's singular values
1956 (again, followed by zeros). [1]_
1958 References
1959 ----------
1960 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
1961 FL, Academic Press, Inc., 1980, pp. 139-142.
1963 Examples
1964 --------
1965 The following example checks that ``a * a+ * a == a`` and
1966 ``a+ * a * a+ == a+``:
1968 >>> a = np.random.randn(9, 6)
1969 >>> B = np.linalg.pinv(a)
1970 >>> np.allclose(a, np.dot(a, np.dot(B, a)))
1971 True
1972 >>> np.allclose(B, np.dot(B, np.dot(a, B)))
1973 True
1975 """
1976 a, wrap = _makearray(a)
1977 rcond = asarray(rcond)
1978 if _is_empty_2d(a):
1979 m, n = a.shape[-2:]
1980 res = empty(a.shape[:-2] + (n, m), dtype=a.dtype)
1981 return wrap(res)
1982 a = a.conjugate()
1983 u, s, vt = svd(a, full_matrices=False, hermitian=hermitian)
1985 # discard small singular values
1986 cutoff = rcond[..., newaxis] * amax(s, axis=-1, keepdims=True)
1987 large = s > cutoff
1988 s = divide(1, s, where=large, out=s)
1989 s[~large] = 0
1991 res = matmul(transpose(vt), multiply(s[..., newaxis], transpose(u)))
1992 return wrap(res)
1995# Determinant
1998@array_function_dispatch(_unary_dispatcher)
1999def slogdet(a):
2000 """
2001 Compute the sign and (natural) logarithm of the determinant of an array.
2003 If an array has a very small or very large determinant, then a call to
2004 `det` may overflow or underflow. This routine is more robust against such
2005 issues, because it computes the logarithm of the determinant rather than
2006 the determinant itself.
2008 Parameters
2009 ----------
2010 a : (..., M, M) array_like
2011 Input array, has to be a square 2-D array.
2013 Returns
2014 -------
2015 sign : (...) array_like
2016 A number representing the sign of the determinant. For a real matrix,
2017 this is 1, 0, or -1. For a complex matrix, this is a complex number
2018 with absolute value 1 (i.e., it is on the unit circle), or else 0.
2019 logdet : (...) array_like
2020 The natural log of the absolute value of the determinant.
2022 If the determinant is zero, then `sign` will be 0 and `logdet` will be
2023 -Inf. In all cases, the determinant is equal to ``sign * np.exp(logdet)``.
2025 See Also
2026 --------
2027 det
2029 Notes
2030 -----
2032 .. versionadded:: 1.8.0
2034 Broadcasting rules apply, see the `numpy.linalg` documentation for
2035 details.
2037 .. versionadded:: 1.6.0
2039 The determinant is computed via LU factorization using the LAPACK
2040 routine ``z/dgetrf``.
2043 Examples
2044 --------
2045 The determinant of a 2-D array ``[[a, b], [c, d]]`` is ``ad - bc``:
2047 >>> a = np.array([[1, 2], [3, 4]])
2048 >>> (sign, logdet) = np.linalg.slogdet(a)
2049 >>> (sign, logdet)
2050 (-1, 0.69314718055994529) # may vary
2051 >>> sign * np.exp(logdet)
2052 -2.0
2054 Computing log-determinants for a stack of matrices:
2056 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2057 >>> a.shape
2058 (3, 2, 2)
2059 >>> sign, logdet = np.linalg.slogdet(a)
2060 >>> (sign, logdet)
2061 (array([-1., -1., -1.]), array([ 0.69314718, 1.09861229, 2.07944154]))
2062 >>> sign * np.exp(logdet)
2063 array([-2., -3., -8.])
2065 This routine succeeds where ordinary `det` does not:
2067 >>> np.linalg.det(np.eye(500) * 0.1)
2068 0.0
2069 >>> np.linalg.slogdet(np.eye(500) * 0.1)
2070 (1, -1151.2925464970228)
2072 """
2073 a = asarray(a)
2074 _assert_stacked_2d(a)
2075 _assert_stacked_square(a)
2076 t, result_t = _commonType(a)
2077 real_t = _realType(result_t)
2078 signature = 'D->Dd' if isComplexType(t) else 'd->dd'
2079 sign, logdet = _umath_linalg.slogdet(a, signature=signature)
2080 sign = sign.astype(result_t, copy=False)
2081 logdet = logdet.astype(real_t, copy=False)
2082 return sign, logdet
2085@array_function_dispatch(_unary_dispatcher)
2086def det(a):
2087 """
2088 Compute the determinant of an array.
2090 Parameters
2091 ----------
2092 a : (..., M, M) array_like
2093 Input array to compute determinants for.
2095 Returns
2096 -------
2097 det : (...) array_like
2098 Determinant of `a`.
2100 See Also
2101 --------
2102 slogdet : Another way to represent the determinant, more suitable
2103 for large matrices where underflow/overflow may occur.
2104 scipy.linalg.det : Similar function in SciPy.
2106 Notes
2107 -----
2109 .. versionadded:: 1.8.0
2111 Broadcasting rules apply, see the `numpy.linalg` documentation for
2112 details.
2114 The determinant is computed via LU factorization using the LAPACK
2115 routine ``z/dgetrf``.
2117 Examples
2118 --------
2119 The determinant of a 2-D array [[a, b], [c, d]] is ad - bc:
2121 >>> a = np.array([[1, 2], [3, 4]])
2122 >>> np.linalg.det(a)
2123 -2.0 # may vary
2125 Computing determinants for a stack of matrices:
2127 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2128 >>> a.shape
2129 (3, 2, 2)
2130 >>> np.linalg.det(a)
2131 array([-2., -3., -8.])
2133 """
2134 a = asarray(a)
2135 _assert_stacked_2d(a)
2136 _assert_stacked_square(a)
2137 t, result_t = _commonType(a)
2138 signature = 'D->D' if isComplexType(t) else 'd->d'
2139 r = _umath_linalg.det(a, signature=signature)
2140 r = r.astype(result_t, copy=False)
2141 return r
2144# Linear Least Squares
2146def _lstsq_dispatcher(a, b, rcond=None):
2147 return (a, b)
2150@array_function_dispatch(_lstsq_dispatcher)
2151def lstsq(a, b, rcond="warn"):
2152 r"""
2153 Return the least-squares solution to a linear matrix equation.
2155 Computes the vector `x` that approximately solves the equation
2156 ``a @ x = b``. The equation may be under-, well-, or over-determined
2157 (i.e., the number of linearly independent rows of `a` can be less than,
2158 equal to, or greater than its number of linearly independent columns).
2159 If `a` is square and of full rank, then `x` (but for round-off error)
2160 is the "exact" solution of the equation. Else, `x` minimizes the
2161 Euclidean 2-norm :math:`||b - ax||`. If there are multiple minimizing
2162 solutions, the one with the smallest 2-norm :math:`||x||` is returned.
2164 Parameters
2165 ----------
2166 a : (M, N) array_like
2167 "Coefficient" matrix.
2168 b : {(M,), (M, K)} array_like
2169 Ordinate or "dependent variable" values. If `b` is two-dimensional,
2170 the least-squares solution is calculated for each of the `K` columns
2171 of `b`.
2172 rcond : float, optional
2173 Cut-off ratio for small singular values of `a`.
2174 For the purposes of rank determination, singular values are treated
2175 as zero if they are smaller than `rcond` times the largest singular
2176 value of `a`.
2178 .. versionchanged:: 1.14.0
2179 If not set, a FutureWarning is given. The previous default
2180 of ``-1`` will use the machine precision as `rcond` parameter,
2181 the new default will use the machine precision times `max(M, N)`.
2182 To silence the warning and use the new default, use ``rcond=None``,
2183 to keep using the old behavior, use ``rcond=-1``.
2185 Returns
2186 -------
2187 x : {(N,), (N, K)} ndarray
2188 Least-squares solution. If `b` is two-dimensional,
2189 the solutions are in the `K` columns of `x`.
2190 residuals : {(1,), (K,), (0,)} ndarray
2191 Sums of squared residuals: Squared Euclidean 2-norm for each column in
2192 ``b - a @ x``.
2193 If the rank of `a` is < N or M <= N, this is an empty array.
2194 If `b` is 1-dimensional, this is a (1,) shape array.
2195 Otherwise the shape is (K,).
2196 rank : int
2197 Rank of matrix `a`.
2198 s : (min(M, N),) ndarray
2199 Singular values of `a`.
2201 Raises
2202 ------
2203 LinAlgError
2204 If computation does not converge.
2206 See Also
2207 --------
2208 scipy.linalg.lstsq : Similar function in SciPy.
2210 Notes
2211 -----
2212 If `b` is a matrix, then all array results are returned as matrices.
2214 Examples
2215 --------
2216 Fit a line, ``y = mx + c``, through some noisy data-points:
2218 >>> x = np.array([0, 1, 2, 3])
2219 >>> y = np.array([-1, 0.2, 0.9, 2.1])
2221 By examining the coefficients, we see that the line should have a
2222 gradient of roughly 1 and cut the y-axis at, more or less, -1.
2224 We can rewrite the line equation as ``y = Ap``, where ``A = [[x 1]]``
2225 and ``p = [[m], [c]]``. Now use `lstsq` to solve for `p`:
2227 >>> A = np.vstack([x, np.ones(len(x))]).T
2228 >>> A
2229 array([[ 0., 1.],
2230 [ 1., 1.],
2231 [ 2., 1.],
2232 [ 3., 1.]])
2234 >>> m, c = np.linalg.lstsq(A, y, rcond=None)[0]
2235 >>> m, c
2236 (1.0 -0.95) # may vary
2238 Plot the data along with the fitted line:
2240 >>> import matplotlib.pyplot as plt
2241 >>> _ = plt.plot(x, y, 'o', label='Original data', markersize=10)
2242 >>> _ = plt.plot(x, m*x + c, 'r', label='Fitted line')
2243 >>> _ = plt.legend()
2244 >>> plt.show()
2246 """
2247 a, _ = _makearray(a)
2248 b, wrap = _makearray(b)
2249 is_1d = b.ndim == 1
2250 if is_1d:
2251 b = b[:, newaxis]
2252 _assert_2d(a, b)
2253 m, n = a.shape[-2:]
2254 m2, n_rhs = b.shape[-2:]
2255 if m != m2:
2256 raise LinAlgError('Incompatible dimensions')
2258 t, result_t = _commonType(a, b)
2259 result_real_t = _realType(result_t)
2261 # Determine default rcond value
2262 if rcond == "warn":
2263 # 2017-08-19, 1.14.0
2264 warnings.warn("`rcond` parameter will change to the default of "
2265 "machine precision times ``max(M, N)`` where M and N "
2266 "are the input matrix dimensions.\n"
2267 "To use the future default and silence this warning "
2268 "we advise to pass `rcond=None`, to keep using the old, "
2269 "explicitly pass `rcond=-1`.",
2270 FutureWarning, stacklevel=3)
2271 rcond = -1
2272 if rcond is None:
2273 rcond = finfo(t).eps * max(n, m)
2275 if m <= n:
2276 gufunc = _umath_linalg.lstsq_m
2277 else:
2278 gufunc = _umath_linalg.lstsq_n
2280 signature = 'DDd->Ddid' if isComplexType(t) else 'ddd->ddid'
2281 extobj = get_linalg_error_extobj(_raise_linalgerror_lstsq)
2282 if n_rhs == 0:
2283 # lapack can't handle n_rhs = 0 - so allocate the array one larger in that axis
2284 b = zeros(b.shape[:-2] + (m, n_rhs + 1), dtype=b.dtype)
2285 x, resids, rank, s = gufunc(a, b, rcond, signature=signature, extobj=extobj)
2286 if m == 0:
2287 x[...] = 0
2288 if n_rhs == 0:
2289 # remove the item we added
2290 x = x[..., :n_rhs]
2291 resids = resids[..., :n_rhs]
2293 # remove the axis we added
2294 if is_1d:
2295 x = x.squeeze(axis=-1)
2296 # we probably should squeeze resids too, but we can't
2297 # without breaking compatibility.
2299 # as documented
2300 if rank != n or m <= n:
2301 resids = array([], result_real_t)
2303 # coerce output arrays
2304 s = s.astype(result_real_t, copy=False)
2305 resids = resids.astype(result_real_t, copy=False)
2306 x = x.astype(result_t, copy=True) # Copying lets the memory in r_parts be freed
2307 return wrap(x), wrap(resids), rank, s
2310def _multi_svd_norm(x, row_axis, col_axis, op):
2311 """Compute a function of the singular values of the 2-D matrices in `x`.
2313 This is a private utility function used by `numpy.linalg.norm()`.
2315 Parameters
2316 ----------
2317 x : ndarray
2318 row_axis, col_axis : int
2319 The axes of `x` that hold the 2-D matrices.
2320 op : callable
2321 This should be either numpy.amin or `numpy.amax` or `numpy.sum`.
2323 Returns
2324 -------
2325 result : float or ndarray
2326 If `x` is 2-D, the return values is a float.
2327 Otherwise, it is an array with ``x.ndim - 2`` dimensions.
2328 The return values are either the minimum or maximum or sum of the
2329 singular values of the matrices, depending on whether `op`
2330 is `numpy.amin` or `numpy.amax` or `numpy.sum`.
2332 """
2333 y = moveaxis(x, (row_axis, col_axis), (-2, -1))
2334 result = op(svd(y, compute_uv=False), axis=-1)
2335 return result
2338def _norm_dispatcher(x, ord=None, axis=None, keepdims=None):
2339 return (x,)
2342@array_function_dispatch(_norm_dispatcher)
2343def norm(x, ord=None, axis=None, keepdims=False):
2344 """
2345 Matrix or vector norm.
2347 This function is able to return one of eight different matrix norms,
2348 or one of an infinite number of vector norms (described below), depending
2349 on the value of the ``ord`` parameter.
2351 Parameters
2352 ----------
2353 x : array_like
2354 Input array. If `axis` is None, `x` must be 1-D or 2-D, unless `ord`
2355 is None. If both `axis` and `ord` are None, the 2-norm of
2356 ``x.ravel`` will be returned.
2357 ord : {non-zero int, inf, -inf, 'fro', 'nuc'}, optional
2358 Order of the norm (see table under ``Notes``). inf means numpy's
2359 `inf` object. The default is None.
2360 axis : {None, int, 2-tuple of ints}, optional.
2361 If `axis` is an integer, it specifies the axis of `x` along which to
2362 compute the vector norms. If `axis` is a 2-tuple, it specifies the
2363 axes that hold 2-D matrices, and the matrix norms of these matrices
2364 are computed. If `axis` is None then either a vector norm (when `x`
2365 is 1-D) or a matrix norm (when `x` is 2-D) is returned. The default
2366 is None.
2368 .. versionadded:: 1.8.0
2370 keepdims : bool, optional
2371 If this is set to True, the axes which are normed over are left in the
2372 result as dimensions with size one. With this option the result will
2373 broadcast correctly against the original `x`.
2375 .. versionadded:: 1.10.0
2377 Returns
2378 -------
2379 n : float or ndarray
2380 Norm of the matrix or vector(s).
2382 See Also
2383 --------
2384 scipy.linalg.norm : Similar function in SciPy.
2386 Notes
2387 -----
2388 For values of ``ord < 1``, the result is, strictly speaking, not a
2389 mathematical 'norm', but it may still be useful for various numerical
2390 purposes.
2392 The following norms can be calculated:
2394 ===== ============================ ==========================
2395 ord norm for matrices norm for vectors
2396 ===== ============================ ==========================
2397 None Frobenius norm 2-norm
2398 'fro' Frobenius norm --
2399 'nuc' nuclear norm --
2400 inf max(sum(abs(x), axis=1)) max(abs(x))
2401 -inf min(sum(abs(x), axis=1)) min(abs(x))
2402 0 -- sum(x != 0)
2403 1 max(sum(abs(x), axis=0)) as below
2404 -1 min(sum(abs(x), axis=0)) as below
2405 2 2-norm (largest sing. value) as below
2406 -2 smallest singular value as below
2407 other -- sum(abs(x)**ord)**(1./ord)
2408 ===== ============================ ==========================
2410 The Frobenius norm is given by [1]_:
2412 :math:`||A||_F = [\\sum_{i,j} abs(a_{i,j})^2]^{1/2}`
2414 The nuclear norm is the sum of the singular values.
2416 Both the Frobenius and nuclear norm orders are only defined for
2417 matrices and raise a ValueError when ``x.ndim != 2``.
2419 References
2420 ----------
2421 .. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*,
2422 Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15
2424 Examples
2425 --------
2426 >>> from numpy import linalg as LA
2427 >>> a = np.arange(9) - 4
2428 >>> a
2429 array([-4, -3, -2, ..., 2, 3, 4])
2430 >>> b = a.reshape((3, 3))
2431 >>> b
2432 array([[-4, -3, -2],
2433 [-1, 0, 1],
2434 [ 2, 3, 4]])
2436 >>> LA.norm(a)
2437 7.745966692414834
2438 >>> LA.norm(b)
2439 7.745966692414834
2440 >>> LA.norm(b, 'fro')
2441 7.745966692414834
2442 >>> LA.norm(a, np.inf)
2443 4.0
2444 >>> LA.norm(b, np.inf)
2445 9.0
2446 >>> LA.norm(a, -np.inf)
2447 0.0
2448 >>> LA.norm(b, -np.inf)
2449 2.0
2451 >>> LA.norm(a, 1)
2452 20.0
2453 >>> LA.norm(b, 1)
2454 7.0
2455 >>> LA.norm(a, -1)
2456 -4.6566128774142013e-010
2457 >>> LA.norm(b, -1)
2458 6.0
2459 >>> LA.norm(a, 2)
2460 7.745966692414834
2461 >>> LA.norm(b, 2)
2462 7.3484692283495345
2464 >>> LA.norm(a, -2)
2465 0.0
2466 >>> LA.norm(b, -2)
2467 1.8570331885190563e-016 # may vary
2468 >>> LA.norm(a, 3)
2469 5.8480354764257312 # may vary
2470 >>> LA.norm(a, -3)
2471 0.0
2473 Using the `axis` argument to compute vector norms:
2475 >>> c = np.array([[ 1, 2, 3],
2476 ... [-1, 1, 4]])
2477 >>> LA.norm(c, axis=0)
2478 array([ 1.41421356, 2.23606798, 5. ])
2479 >>> LA.norm(c, axis=1)
2480 array([ 3.74165739, 4.24264069])
2481 >>> LA.norm(c, ord=1, axis=1)
2482 array([ 6., 6.])
2484 Using the `axis` argument to compute matrix norms:
2486 >>> m = np.arange(8).reshape(2,2,2)
2487 >>> LA.norm(m, axis=(1,2))
2488 array([ 3.74165739, 11.22497216])
2489 >>> LA.norm(m[0, :, :]), LA.norm(m[1, :, :])
2490 (3.7416573867739413, 11.224972160321824)
2492 """
2493 x = asarray(x)
2495 if not issubclass(x.dtype.type, (inexact, object_)):
2496 x = x.astype(float)
2498 # Immediately handle some default, simple, fast, and common cases.
2499 if axis is None:
2500 ndim = x.ndim
2501 if ((ord is None) or
2502 (ord in ('f', 'fro') and ndim == 2) or
2503 (ord == 2 and ndim == 1)):
2505 x = x.ravel(order='K')
2506 if isComplexType(x.dtype.type):
2507 x_real = x.real
2508 x_imag = x.imag
2509 sqnorm = x_real.dot(x_real) + x_imag.dot(x_imag)
2510 else:
2511 sqnorm = x.dot(x)
2512 ret = sqrt(sqnorm)
2513 if keepdims:
2514 ret = ret.reshape(ndim*[1])
2515 return ret
2517 # Normalize the `axis` argument to a tuple.
2518 nd = x.ndim
2519 if axis is None:
2520 axis = tuple(range(nd))
2521 elif not isinstance(axis, tuple):
2522 try:
2523 axis = int(axis)
2524 except Exception as e:
2525 raise TypeError("'axis' must be None, an integer or a tuple of integers") from e
2526 axis = (axis,)
2528 if len(axis) == 1:
2529 if ord == Inf:
2530 return abs(x).max(axis=axis, keepdims=keepdims)
2531 elif ord == -Inf:
2532 return abs(x).min(axis=axis, keepdims=keepdims)
2533 elif ord == 0:
2534 # Zero norm
2535 return (x != 0).astype(x.real.dtype).sum(axis=axis, keepdims=keepdims)
2536 elif ord == 1:
2537 # special case for speedup
2538 return add.reduce(abs(x), axis=axis, keepdims=keepdims)
2539 elif ord is None or ord == 2:
2540 # special case for speedup
2541 s = (x.conj() * x).real
2542 return sqrt(add.reduce(s, axis=axis, keepdims=keepdims))
2543 # None of the str-type keywords for ord ('fro', 'nuc')
2544 # are valid for vectors
2545 elif isinstance(ord, str):
2546 raise ValueError(f"Invalid norm order '{ord}' for vectors")
2547 else:
2548 absx = abs(x)
2549 absx **= ord
2550 ret = add.reduce(absx, axis=axis, keepdims=keepdims)
2551 ret **= reciprocal(ord, dtype=ret.dtype)
2552 return ret
2553 elif len(axis) == 2:
2554 row_axis, col_axis = axis
2555 row_axis = normalize_axis_index(row_axis, nd)
2556 col_axis = normalize_axis_index(col_axis, nd)
2557 if row_axis == col_axis:
2558 raise ValueError('Duplicate axes given.')
2559 if ord == 2:
2560 ret = _multi_svd_norm(x, row_axis, col_axis, amax)
2561 elif ord == -2:
2562 ret = _multi_svd_norm(x, row_axis, col_axis, amin)
2563 elif ord == 1:
2564 if col_axis > row_axis:
2565 col_axis -= 1
2566 ret = add.reduce(abs(x), axis=row_axis).max(axis=col_axis)
2567 elif ord == Inf:
2568 if row_axis > col_axis:
2569 row_axis -= 1
2570 ret = add.reduce(abs(x), axis=col_axis).max(axis=row_axis)
2571 elif ord == -1:
2572 if col_axis > row_axis:
2573 col_axis -= 1
2574 ret = add.reduce(abs(x), axis=row_axis).min(axis=col_axis)
2575 elif ord == -Inf:
2576 if row_axis > col_axis:
2577 row_axis -= 1
2578 ret = add.reduce(abs(x), axis=col_axis).min(axis=row_axis)
2579 elif ord in [None, 'fro', 'f']:
2580 ret = sqrt(add.reduce((x.conj() * x).real, axis=axis))
2581 elif ord == 'nuc':
2582 ret = _multi_svd_norm(x, row_axis, col_axis, sum)
2583 else:
2584 raise ValueError("Invalid norm order for matrices.")
2585 if keepdims:
2586 ret_shape = list(x.shape)
2587 ret_shape[axis[0]] = 1
2588 ret_shape[axis[1]] = 1
2589 ret = ret.reshape(ret_shape)
2590 return ret
2591 else:
2592 raise ValueError("Improper number of dimensions to norm.")
2595# multi_dot
2597def _multidot_dispatcher(arrays, *, out=None):
2598 yield from arrays
2599 yield out
2602@array_function_dispatch(_multidot_dispatcher)
2603def multi_dot(arrays, *, out=None):
2604 """
2605 Compute the dot product of two or more arrays in a single function call,
2606 while automatically selecting the fastest evaluation order.
2608 `multi_dot` chains `numpy.dot` and uses optimal parenthesization
2609 of the matrices [1]_ [2]_. Depending on the shapes of the matrices,
2610 this can speed up the multiplication a lot.
2612 If the first argument is 1-D it is treated as a row vector.
2613 If the last argument is 1-D it is treated as a column vector.
2614 The other arguments must be 2-D.
2616 Think of `multi_dot` as::
2618 def multi_dot(arrays): return functools.reduce(np.dot, arrays)
2621 Parameters
2622 ----------
2623 arrays : sequence of array_like
2624 If the first argument is 1-D it is treated as row vector.
2625 If the last argument is 1-D it is treated as column vector.
2626 The other arguments must be 2-D.
2627 out : ndarray, optional
2628 Output argument. This must have the exact kind that would be returned
2629 if it was not used. In particular, it must have the right type, must be
2630 C-contiguous, and its dtype must be the dtype that would be returned
2631 for `dot(a, b)`. This is a performance feature. Therefore, if these
2632 conditions are not met, an exception is raised, instead of attempting
2633 to be flexible.
2635 .. versionadded:: 1.19.0
2637 Returns
2638 -------
2639 output : ndarray
2640 Returns the dot product of the supplied arrays.
2642 See Also
2643 --------
2644 numpy.dot : dot multiplication with two arguments.
2646 References
2647 ----------
2649 .. [1] Cormen, "Introduction to Algorithms", Chapter 15.2, p. 370-378
2650 .. [2] https://en.wikipedia.org/wiki/Matrix_chain_multiplication
2652 Examples
2653 --------
2654 `multi_dot` allows you to write::
2656 >>> from numpy.linalg import multi_dot
2657 >>> # Prepare some data
2658 >>> A = np.random.random((10000, 100))
2659 >>> B = np.random.random((100, 1000))
2660 >>> C = np.random.random((1000, 5))
2661 >>> D = np.random.random((5, 333))
2662 >>> # the actual dot multiplication
2663 >>> _ = multi_dot([A, B, C, D])
2665 instead of::
2667 >>> _ = np.dot(np.dot(np.dot(A, B), C), D)
2668 >>> # or
2669 >>> _ = A.dot(B).dot(C).dot(D)
2671 Notes
2672 -----
2673 The cost for a matrix multiplication can be calculated with the
2674 following function::
2676 def cost(A, B):
2677 return A.shape[0] * A.shape[1] * B.shape[1]
2679 Assume we have three matrices
2680 :math:`A_{10x100}, B_{100x5}, C_{5x50}`.
2682 The costs for the two different parenthesizations are as follows::
2684 cost((AB)C) = 10*100*5 + 10*5*50 = 5000 + 2500 = 7500
2685 cost(A(BC)) = 10*100*50 + 100*5*50 = 50000 + 25000 = 75000
2687 """
2688 n = len(arrays)
2689 # optimization only makes sense for len(arrays) > 2
2690 if n < 2:
2691 raise ValueError("Expecting at least two arrays.")
2692 elif n == 2:
2693 return dot(arrays[0], arrays[1], out=out)
2695 arrays = [asanyarray(a) for a in arrays]
2697 # save original ndim to reshape the result array into the proper form later
2698 ndim_first, ndim_last = arrays[0].ndim, arrays[-1].ndim
2699 # Explicitly convert vectors to 2D arrays to keep the logic of the internal
2700 # _multi_dot_* functions as simple as possible.
2701 if arrays[0].ndim == 1:
2702 arrays[0] = atleast_2d(arrays[0])
2703 if arrays[-1].ndim == 1:
2704 arrays[-1] = atleast_2d(arrays[-1]).T
2705 _assert_2d(*arrays)
2707 # _multi_dot_three is much faster than _multi_dot_matrix_chain_order
2708 if n == 3:
2709 result = _multi_dot_three(arrays[0], arrays[1], arrays[2], out=out)
2710 else:
2711 order = _multi_dot_matrix_chain_order(arrays)
2712 result = _multi_dot(arrays, order, 0, n - 1, out=out)
2714 # return proper shape
2715 if ndim_first == 1 and ndim_last == 1:
2716 return result[0, 0] # scalar
2717 elif ndim_first == 1 or ndim_last == 1:
2718 return result.ravel() # 1-D
2719 else:
2720 return result
2723def _multi_dot_three(A, B, C, out=None):
2724 """
2725 Find the best order for three arrays and do the multiplication.
2727 For three arguments `_multi_dot_three` is approximately 15 times faster
2728 than `_multi_dot_matrix_chain_order`
2730 """
2731 a0, a1b0 = A.shape
2732 b1c0, c1 = C.shape
2733 # cost1 = cost((AB)C) = a0*a1b0*b1c0 + a0*b1c0*c1
2734 cost1 = a0 * b1c0 * (a1b0 + c1)
2735 # cost2 = cost(A(BC)) = a1b0*b1c0*c1 + a0*a1b0*c1
2736 cost2 = a1b0 * c1 * (a0 + b1c0)
2738 if cost1 < cost2:
2739 return dot(dot(A, B), C, out=out)
2740 else:
2741 return dot(A, dot(B, C), out=out)
2744def _multi_dot_matrix_chain_order(arrays, return_costs=False):
2745 """
2746 Return a np.array that encodes the optimal order of mutiplications.
2748 The optimal order array is then used by `_multi_dot()` to do the
2749 multiplication.
2751 Also return the cost matrix if `return_costs` is `True`
2753 The implementation CLOSELY follows Cormen, "Introduction to Algorithms",
2754 Chapter 15.2, p. 370-378. Note that Cormen uses 1-based indices.
2756 cost[i, j] = min([
2757 cost[prefix] + cost[suffix] + cost_mult(prefix, suffix)
2758 for k in range(i, j)])
2760 """
2761 n = len(arrays)
2762 # p stores the dimensions of the matrices
2763 # Example for p: A_{10x100}, B_{100x5}, C_{5x50} --> p = [10, 100, 5, 50]
2764 p = [a.shape[0] for a in arrays] + [arrays[-1].shape[1]]
2765 # m is a matrix of costs of the subproblems
2766 # m[i,j]: min number of scalar multiplications needed to compute A_{i..j}
2767 m = zeros((n, n), dtype=double)
2768 # s is the actual ordering
2769 # s[i, j] is the value of k at which we split the product A_i..A_j
2770 s = empty((n, n), dtype=intp)
2772 for l in range(1, n):
2773 for i in range(n - l):
2774 j = i + l
2775 m[i, j] = Inf
2776 for k in range(i, j):
2777 q = m[i, k] + m[k+1, j] + p[i]*p[k+1]*p[j+1]
2778 if q < m[i, j]:
2779 m[i, j] = q
2780 s[i, j] = k # Note that Cormen uses 1-based index
2782 return (s, m) if return_costs else s
2785def _multi_dot(arrays, order, i, j, out=None):
2786 """Actually do the multiplication with the given order."""
2787 if i == j:
2788 # the initial call with non-None out should never get here
2789 assert out is None
2791 return arrays[i]
2792 else:
2793 return dot(_multi_dot(arrays, order, i, order[i, j]),
2794 _multi_dot(arrays, order, order[i, j] + 1, j),
2795 out=out)