Coverage for /var/srv/projects/api.amasfac.comuna18.com/tmp/venv/lib/python3.9/site-packages/numpy/linalg/linalg.py: 12%
662 statements
« prev ^ index » next coverage.py v6.4.4, created at 2023-07-17 14:22 -0600
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, fastCopyAndTranspose, sum, isfinite,
25 finfo, errstate, geterrobj, moveaxis, amin, amax, product, abs,
26 atleast_2d, intp, asanyarray, object_, matmul,
27 swapaxes, divide, count_nonzero, isnan, sign, argsort, sort,
28 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
161# _fastCopyAndTranpose assumes the input is 2D (as all the calls in here are).
163_fastCT = fastCopyAndTranspose
165def _to_native_byte_order(*arrays):
166 ret = []
167 for arr in arrays:
168 if arr.dtype.byteorder not in ('=', '|'):
169 ret.append(asarray(arr, dtype=arr.dtype.newbyteorder('=')))
170 else:
171 ret.append(arr)
172 if len(ret) == 1:
173 return ret[0]
174 else:
175 return ret
177def _fastCopyAndTranspose(type, *arrays):
178 cast_arrays = ()
179 for a in arrays:
180 if a.dtype.type is not type:
181 a = a.astype(type)
182 cast_arrays = cast_arrays + (_fastCT(a),)
183 if len(cast_arrays) == 1:
184 return cast_arrays[0]
185 else:
186 return cast_arrays
188def _assert_2d(*arrays):
189 for a in arrays:
190 if a.ndim != 2:
191 raise LinAlgError('%d-dimensional array given. Array must be '
192 'two-dimensional' % a.ndim)
194def _assert_stacked_2d(*arrays):
195 for a in arrays:
196 if a.ndim < 2:
197 raise LinAlgError('%d-dimensional array given. Array must be '
198 'at least two-dimensional' % a.ndim)
200def _assert_stacked_square(*arrays):
201 for a in arrays:
202 m, n = a.shape[-2:]
203 if m != n:
204 raise LinAlgError('Last 2 dimensions of the array must be square')
206def _assert_finite(*arrays):
207 for a in arrays:
208 if not isfinite(a).all():
209 raise LinAlgError("Array must not contain infs or NaNs")
211def _is_empty_2d(arr):
212 # check size first for efficiency
213 return arr.size == 0 and product(arr.shape[-2:]) == 0
216def transpose(a):
217 """
218 Transpose each matrix in a stack of matrices.
220 Unlike np.transpose, this only swaps the last two axes, rather than all of
221 them
223 Parameters
224 ----------
225 a : (...,M,N) array_like
227 Returns
228 -------
229 aT : (...,N,M) ndarray
230 """
231 return swapaxes(a, -1, -2)
233# Linear equations
235def _tensorsolve_dispatcher(a, b, axes=None):
236 return (a, b)
239@array_function_dispatch(_tensorsolve_dispatcher)
240def tensorsolve(a, b, axes=None):
241 """
242 Solve the tensor equation ``a x = b`` for x.
244 It is assumed that all indices of `x` are summed over in the product,
245 together with the rightmost indices of `a`, as is done in, for example,
246 ``tensordot(a, x, axes=b.ndim)``.
248 Parameters
249 ----------
250 a : array_like
251 Coefficient tensor, of shape ``b.shape + Q``. `Q`, a tuple, equals
252 the shape of that sub-tensor of `a` consisting of the appropriate
253 number of its rightmost indices, and must be such that
254 ``prod(Q) == prod(b.shape)`` (in which sense `a` is said to be
255 'square').
256 b : array_like
257 Right-hand tensor, which can be of any shape.
258 axes : tuple of ints, optional
259 Axes in `a` to reorder to the right, before inversion.
260 If None (default), no reordering is done.
262 Returns
263 -------
264 x : ndarray, shape Q
266 Raises
267 ------
268 LinAlgError
269 If `a` is singular or not 'square' (in the above sense).
271 See Also
272 --------
273 numpy.tensordot, tensorinv, numpy.einsum
275 Examples
276 --------
277 >>> a = np.eye(2*3*4)
278 >>> a.shape = (2*3, 4, 2, 3, 4)
279 >>> b = np.random.randn(2*3, 4)
280 >>> x = np.linalg.tensorsolve(a, b)
281 >>> x.shape
282 (2, 3, 4)
283 >>> np.allclose(np.tensordot(a, x, axes=3), b)
284 True
286 """
287 a, wrap = _makearray(a)
288 b = asarray(b)
289 an = a.ndim
291 if axes is not None:
292 allaxes = list(range(0, an))
293 for k in axes:
294 allaxes.remove(k)
295 allaxes.insert(an, k)
296 a = a.transpose(allaxes)
298 oldshape = a.shape[-(an-b.ndim):]
299 prod = 1
300 for k in oldshape:
301 prod *= k
303 if a.size != prod ** 2:
304 raise LinAlgError(
305 "Input arrays must satisfy the requirement \
306 prod(a.shape[b.ndim:]) == prod(a.shape[:b.ndim])"
307 )
309 a = a.reshape(prod, prod)
310 b = b.ravel()
311 res = wrap(solve(a, b))
312 res.shape = oldshape
313 return res
316def _solve_dispatcher(a, b):
317 return (a, b)
320@array_function_dispatch(_solve_dispatcher)
321def solve(a, b):
322 """
323 Solve a linear matrix equation, or system of linear scalar equations.
325 Computes the "exact" solution, `x`, of the well-determined, i.e., full
326 rank, linear matrix equation `ax = b`.
328 Parameters
329 ----------
330 a : (..., M, M) array_like
331 Coefficient matrix.
332 b : {(..., M,), (..., M, K)}, array_like
333 Ordinate or "dependent variable" values.
335 Returns
336 -------
337 x : {(..., M,), (..., M, K)} ndarray
338 Solution to the system a x = b. Returned shape is identical to `b`.
340 Raises
341 ------
342 LinAlgError
343 If `a` is singular or not square.
345 See Also
346 --------
347 scipy.linalg.solve : Similar function in SciPy.
349 Notes
350 -----
352 .. versionadded:: 1.8.0
354 Broadcasting rules apply, see the `numpy.linalg` documentation for
355 details.
357 The solutions are computed using LAPACK routine ``_gesv``.
359 `a` must be square and of full-rank, i.e., all rows (or, equivalently,
360 columns) must be linearly independent; if either is not true, use
361 `lstsq` for the least-squares best "solution" of the
362 system/equation.
364 References
365 ----------
366 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
367 FL, Academic Press, Inc., 1980, pg. 22.
369 Examples
370 --------
371 Solve the system of equations ``x0 + 2 * x1 = 1`` and ``3 * x0 + 5 * x1 = 2``:
373 >>> a = np.array([[1, 2], [3, 5]])
374 >>> b = np.array([1, 2])
375 >>> x = np.linalg.solve(a, b)
376 >>> x
377 array([-1., 1.])
379 Check that the solution is correct:
381 >>> np.allclose(np.dot(a, x), b)
382 True
384 """
385 a, _ = _makearray(a)
386 _assert_stacked_2d(a)
387 _assert_stacked_square(a)
388 b, wrap = _makearray(b)
389 t, result_t = _commonType(a, b)
391 # We use the b = (..., M,) logic, only if the number of extra dimensions
392 # match exactly
393 if b.ndim == a.ndim - 1:
394 gufunc = _umath_linalg.solve1
395 else:
396 gufunc = _umath_linalg.solve
398 signature = 'DD->D' if isComplexType(t) else 'dd->d'
399 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
400 r = gufunc(a, b, signature=signature, extobj=extobj)
402 return wrap(r.astype(result_t, copy=False))
405def _tensorinv_dispatcher(a, ind=None):
406 return (a,)
409@array_function_dispatch(_tensorinv_dispatcher)
410def tensorinv(a, ind=2):
411 """
412 Compute the 'inverse' of an N-dimensional array.
414 The result is an inverse for `a` relative to the tensordot operation
415 ``tensordot(a, b, ind)``, i. e., up to floating-point accuracy,
416 ``tensordot(tensorinv(a), a, ind)`` is the "identity" tensor for the
417 tensordot operation.
419 Parameters
420 ----------
421 a : array_like
422 Tensor to 'invert'. Its shape must be 'square', i. e.,
423 ``prod(a.shape[:ind]) == prod(a.shape[ind:])``.
424 ind : int, optional
425 Number of first indices that are involved in the inverse sum.
426 Must be a positive integer, default is 2.
428 Returns
429 -------
430 b : ndarray
431 `a`'s tensordot inverse, shape ``a.shape[ind:] + a.shape[:ind]``.
433 Raises
434 ------
435 LinAlgError
436 If `a` is singular or not 'square' (in the above sense).
438 See Also
439 --------
440 numpy.tensordot, tensorsolve
442 Examples
443 --------
444 >>> a = np.eye(4*6)
445 >>> a.shape = (4, 6, 8, 3)
446 >>> ainv = np.linalg.tensorinv(a, ind=2)
447 >>> ainv.shape
448 (8, 3, 4, 6)
449 >>> b = np.random.randn(4, 6)
450 >>> np.allclose(np.tensordot(ainv, b), np.linalg.tensorsolve(a, b))
451 True
453 >>> a = np.eye(4*6)
454 >>> a.shape = (24, 8, 3)
455 >>> ainv = np.linalg.tensorinv(a, ind=1)
456 >>> ainv.shape
457 (8, 3, 24)
458 >>> b = np.random.randn(24)
459 >>> np.allclose(np.tensordot(ainv, b, 1), np.linalg.tensorsolve(a, b))
460 True
462 """
463 a = asarray(a)
464 oldshape = a.shape
465 prod = 1
466 if ind > 0:
467 invshape = oldshape[ind:] + oldshape[:ind]
468 for k in oldshape[ind:]:
469 prod *= k
470 else:
471 raise ValueError("Invalid ind argument.")
472 a = a.reshape(prod, -1)
473 ia = inv(a)
474 return ia.reshape(*invshape)
477# Matrix inversion
479def _unary_dispatcher(a):
480 return (a,)
483@array_function_dispatch(_unary_dispatcher)
484def inv(a):
485 """
486 Compute the (multiplicative) inverse of a matrix.
488 Given a square matrix `a`, return the matrix `ainv` satisfying
489 ``dot(a, ainv) = dot(ainv, a) = eye(a.shape[0])``.
491 Parameters
492 ----------
493 a : (..., M, M) array_like
494 Matrix to be inverted.
496 Returns
497 -------
498 ainv : (..., M, M) ndarray or matrix
499 (Multiplicative) inverse of the matrix `a`.
501 Raises
502 ------
503 LinAlgError
504 If `a` is not square or inversion fails.
506 See Also
507 --------
508 scipy.linalg.inv : Similar function in SciPy.
510 Notes
511 -----
513 .. versionadded:: 1.8.0
515 Broadcasting rules apply, see the `numpy.linalg` documentation for
516 details.
518 Examples
519 --------
520 >>> from numpy.linalg import inv
521 >>> a = np.array([[1., 2.], [3., 4.]])
522 >>> ainv = inv(a)
523 >>> np.allclose(np.dot(a, ainv), np.eye(2))
524 True
525 >>> np.allclose(np.dot(ainv, a), np.eye(2))
526 True
528 If a is a matrix object, then the return value is a matrix as well:
530 >>> ainv = inv(np.matrix(a))
531 >>> ainv
532 matrix([[-2. , 1. ],
533 [ 1.5, -0.5]])
535 Inverses of several matrices can be computed at once:
537 >>> a = np.array([[[1., 2.], [3., 4.]], [[1, 3], [3, 5]]])
538 >>> inv(a)
539 array([[[-2. , 1. ],
540 [ 1.5 , -0.5 ]],
541 [[-1.25, 0.75],
542 [ 0.75, -0.25]]])
544 """
545 a, wrap = _makearray(a)
546 _assert_stacked_2d(a)
547 _assert_stacked_square(a)
548 t, result_t = _commonType(a)
550 signature = 'D->D' if isComplexType(t) else 'd->d'
551 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
552 ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
553 return wrap(ainv.astype(result_t, copy=False))
556def _matrix_power_dispatcher(a, n):
557 return (a,)
560@array_function_dispatch(_matrix_power_dispatcher)
561def matrix_power(a, n):
562 """
563 Raise a square matrix to the (integer) power `n`.
565 For positive integers `n`, the power is computed by repeated matrix
566 squarings and matrix multiplications. If ``n == 0``, the identity matrix
567 of the same shape as M is returned. If ``n < 0``, the inverse
568 is computed and then raised to the ``abs(n)``.
570 .. note:: Stacks of object matrices are not currently supported.
572 Parameters
573 ----------
574 a : (..., M, M) array_like
575 Matrix to be "powered".
576 n : int
577 The exponent can be any integer or long integer, positive,
578 negative, or zero.
580 Returns
581 -------
582 a**n : (..., M, M) ndarray or matrix object
583 The return value is the same shape and type as `M`;
584 if the exponent is positive or zero then the type of the
585 elements is the same as those of `M`. If the exponent is
586 negative the elements are floating-point.
588 Raises
589 ------
590 LinAlgError
591 For matrices that are not square or that (for negative powers) cannot
592 be inverted numerically.
594 Examples
595 --------
596 >>> from numpy.linalg import matrix_power
597 >>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit
598 >>> matrix_power(i, 3) # should = -i
599 array([[ 0, -1],
600 [ 1, 0]])
601 >>> matrix_power(i, 0)
602 array([[1, 0],
603 [0, 1]])
604 >>> matrix_power(i, -3) # should = 1/(-i) = i, but w/ f.p. elements
605 array([[ 0., 1.],
606 [-1., 0.]])
608 Somewhat more sophisticated example
610 >>> q = np.zeros((4, 4))
611 >>> q[0:2, 0:2] = -i
612 >>> q[2:4, 2:4] = i
613 >>> q # one of the three quaternion units not equal to 1
614 array([[ 0., -1., 0., 0.],
615 [ 1., 0., 0., 0.],
616 [ 0., 0., 0., 1.],
617 [ 0., 0., -1., 0.]])
618 >>> matrix_power(q, 2) # = -np.eye(4)
619 array([[-1., 0., 0., 0.],
620 [ 0., -1., 0., 0.],
621 [ 0., 0., -1., 0.],
622 [ 0., 0., 0., -1.]])
624 """
625 a = asanyarray(a)
626 _assert_stacked_2d(a)
627 _assert_stacked_square(a)
629 try:
630 n = operator.index(n)
631 except TypeError as e:
632 raise TypeError("exponent must be an integer") from e
634 # Fall back on dot for object arrays. Object arrays are not supported by
635 # the current implementation of matmul using einsum
636 if a.dtype != object:
637 fmatmul = matmul
638 elif a.ndim == 2:
639 fmatmul = dot
640 else:
641 raise NotImplementedError(
642 "matrix_power not supported for stacks of object arrays")
644 if n == 0:
645 a = empty_like(a)
646 a[...] = eye(a.shape[-2], dtype=a.dtype)
647 return a
649 elif n < 0:
650 a = inv(a)
651 n = abs(n)
653 # short-cuts.
654 if n == 1:
655 return a
657 elif n == 2:
658 return fmatmul(a, a)
660 elif n == 3:
661 return fmatmul(fmatmul(a, a), a)
663 # Use binary decomposition to reduce the number of matrix multiplications.
664 # Here, we iterate over the bits of n, from LSB to MSB, raise `a` to
665 # increasing powers of 2, and multiply into the result as needed.
666 z = result = None
667 while n > 0:
668 z = a if z is None else fmatmul(z, z)
669 n, bit = divmod(n, 2)
670 if bit:
671 result = z if result is None else fmatmul(result, z)
673 return result
676# Cholesky decomposition
679@array_function_dispatch(_unary_dispatcher)
680def cholesky(a):
681 """
682 Cholesky decomposition.
684 Return the Cholesky decomposition, `L * L.H`, of the square matrix `a`,
685 where `L` is lower-triangular and .H is the conjugate transpose operator
686 (which is the ordinary transpose if `a` is real-valued). `a` must be
687 Hermitian (symmetric if real-valued) and positive-definite. No
688 checking is performed to verify whether `a` is Hermitian or not.
689 In addition, only the lower-triangular and diagonal elements of `a`
690 are used. Only `L` is actually returned.
692 Parameters
693 ----------
694 a : (..., M, M) array_like
695 Hermitian (symmetric if all elements are real), positive-definite
696 input matrix.
698 Returns
699 -------
700 L : (..., M, M) array_like
701 Lower-triangular Cholesky factor of `a`. Returns a matrix object if
702 `a` is a matrix object.
704 Raises
705 ------
706 LinAlgError
707 If the decomposition fails, for example, if `a` is not
708 positive-definite.
710 See Also
711 --------
712 scipy.linalg.cholesky : Similar function in SciPy.
713 scipy.linalg.cholesky_banded : Cholesky decompose a banded Hermitian
714 positive-definite matrix.
715 scipy.linalg.cho_factor : Cholesky decomposition of a matrix, to use in
716 `scipy.linalg.cho_solve`.
718 Notes
719 -----
721 .. versionadded:: 1.8.0
723 Broadcasting rules apply, see the `numpy.linalg` documentation for
724 details.
726 The Cholesky decomposition is often used as a fast way of solving
728 .. math:: A \\mathbf{x} = \\mathbf{b}
730 (when `A` is both Hermitian/symmetric and positive-definite).
732 First, we solve for :math:`\\mathbf{y}` in
734 .. math:: L \\mathbf{y} = \\mathbf{b},
736 and then for :math:`\\mathbf{x}` in
738 .. math:: L.H \\mathbf{x} = \\mathbf{y}.
740 Examples
741 --------
742 >>> A = np.array([[1,-2j],[2j,5]])
743 >>> A
744 array([[ 1.+0.j, -0.-2.j],
745 [ 0.+2.j, 5.+0.j]])
746 >>> L = np.linalg.cholesky(A)
747 >>> L
748 array([[1.+0.j, 0.+0.j],
749 [0.+2.j, 1.+0.j]])
750 >>> np.dot(L, L.T.conj()) # verify that L * L.H = A
751 array([[1.+0.j, 0.-2.j],
752 [0.+2.j, 5.+0.j]])
753 >>> A = [[1,-2j],[2j,5]] # what happens if A is only array_like?
754 >>> np.linalg.cholesky(A) # an ndarray object is returned
755 array([[1.+0.j, 0.+0.j],
756 [0.+2.j, 1.+0.j]])
757 >>> # But a matrix object is returned if A is a matrix object
758 >>> np.linalg.cholesky(np.matrix(A))
759 matrix([[ 1.+0.j, 0.+0.j],
760 [ 0.+2.j, 1.+0.j]])
762 """
763 extobj = get_linalg_error_extobj(_raise_linalgerror_nonposdef)
764 gufunc = _umath_linalg.cholesky_lo
765 a, wrap = _makearray(a)
766 _assert_stacked_2d(a)
767 _assert_stacked_square(a)
768 t, result_t = _commonType(a)
769 signature = 'D->D' if isComplexType(t) else 'd->d'
770 r = gufunc(a, signature=signature, extobj=extobj)
771 return wrap(r.astype(result_t, copy=False))
774# QR decomposition
776def _qr_dispatcher(a, mode=None):
777 return (a,)
780@array_function_dispatch(_qr_dispatcher)
781def qr(a, mode='reduced'):
782 """
783 Compute the qr factorization of a matrix.
785 Factor the matrix `a` as *qr*, where `q` is orthonormal and `r` is
786 upper-triangular.
788 Parameters
789 ----------
790 a : array_like, shape (..., M, N)
791 An array-like object with the dimensionality of at least 2.
792 mode : {'reduced', 'complete', 'r', 'raw'}, optional
793 If K = min(M, N), then
795 * 'reduced' : returns q, r with dimensions
796 (..., M, K), (..., K, N) (default)
797 * 'complete' : returns q, r with dimensions (..., M, M), (..., M, N)
798 * 'r' : returns r only with dimensions (..., K, N)
799 * 'raw' : returns h, tau with dimensions (..., N, M), (..., K,)
801 The options 'reduced', 'complete, and 'raw' are new in numpy 1.8,
802 see the notes for more information. The default is 'reduced', and to
803 maintain backward compatibility with earlier versions of numpy both
804 it and the old default 'full' can be omitted. Note that array h
805 returned in 'raw' mode is transposed for calling Fortran. The
806 'economic' mode is deprecated. The modes 'full' and 'economic' may
807 be passed using only the first letter for backwards compatibility,
808 but all others must be spelled out. See the Notes for more
809 explanation.
812 Returns
813 -------
814 q : ndarray of float or complex, optional
815 A matrix with orthonormal columns. When mode = 'complete' the
816 result is an orthogonal/unitary matrix depending on whether or not
817 a is real/complex. The determinant may be either +/- 1 in that
818 case. In case the number of dimensions in the input array is
819 greater than 2 then a stack of the matrices with above properties
820 is returned.
821 r : ndarray of float or complex, optional
822 The upper-triangular matrix or a stack of upper-triangular
823 matrices if the number of dimensions in the input array is greater
824 than 2.
825 (h, tau) : ndarrays of np.double or np.cdouble, optional
826 The array h contains the Householder reflectors that generate q
827 along with r. The tau array contains scaling factors for the
828 reflectors. In the deprecated 'economic' mode only h is returned.
830 Raises
831 ------
832 LinAlgError
833 If factoring fails.
835 See Also
836 --------
837 scipy.linalg.qr : Similar function in SciPy.
838 scipy.linalg.rq : Compute RQ decomposition of a matrix.
840 Notes
841 -----
842 This is an interface to the LAPACK routines ``dgeqrf``, ``zgeqrf``,
843 ``dorgqr``, and ``zungqr``.
845 For more information on the qr factorization, see for example:
846 https://en.wikipedia.org/wiki/QR_factorization
848 Subclasses of `ndarray` are preserved except for the 'raw' mode. So if
849 `a` is of type `matrix`, all the return values will be matrices too.
851 New 'reduced', 'complete', and 'raw' options for mode were added in
852 NumPy 1.8.0 and the old option 'full' was made an alias of 'reduced'. In
853 addition the options 'full' and 'economic' were deprecated. Because
854 'full' was the previous default and 'reduced' is the new default,
855 backward compatibility can be maintained by letting `mode` default.
856 The 'raw' option was added so that LAPACK routines that can multiply
857 arrays by q using the Householder reflectors can be used. Note that in
858 this case the returned arrays are of type np.double or np.cdouble and
859 the h array is transposed to be FORTRAN compatible. No routines using
860 the 'raw' return are currently exposed by numpy, but some are available
861 in lapack_lite and just await the necessary work.
863 Examples
864 --------
865 >>> a = np.random.randn(9, 6)
866 >>> q, r = np.linalg.qr(a)
867 >>> np.allclose(a, np.dot(q, r)) # a does equal qr
868 True
869 >>> r2 = np.linalg.qr(a, mode='r')
870 >>> np.allclose(r, r2) # mode='r' returns the same r as mode='full'
871 True
872 >>> a = np.random.normal(size=(3, 2, 2)) # Stack of 2 x 2 matrices as input
873 >>> q, r = np.linalg.qr(a)
874 >>> q.shape
875 (3, 2, 2)
876 >>> r.shape
877 (3, 2, 2)
878 >>> np.allclose(a, np.matmul(q, r))
879 True
881 Example illustrating a common use of `qr`: solving of least squares
882 problems
884 What are the least-squares-best `m` and `y0` in ``y = y0 + mx`` for
885 the following data: {(0,1), (1,0), (1,2), (2,1)}. (Graph the points
886 and you'll see that it should be y0 = 0, m = 1.) The answer is provided
887 by solving the over-determined matrix equation ``Ax = b``, where::
889 A = array([[0, 1], [1, 1], [1, 1], [2, 1]])
890 x = array([[y0], [m]])
891 b = array([[1], [0], [2], [1]])
893 If A = qr such that q is orthonormal (which is always possible via
894 Gram-Schmidt), then ``x = inv(r) * (q.T) * b``. (In numpy practice,
895 however, we simply use `lstsq`.)
897 >>> A = np.array([[0, 1], [1, 1], [1, 1], [2, 1]])
898 >>> A
899 array([[0, 1],
900 [1, 1],
901 [1, 1],
902 [2, 1]])
903 >>> b = np.array([1, 2, 2, 3])
904 >>> q, r = np.linalg.qr(A)
905 >>> p = np.dot(q.T, b)
906 >>> np.dot(np.linalg.inv(r), p)
907 array([ 1., 1.])
909 """
910 if mode not in ('reduced', 'complete', 'r', 'raw'):
911 if mode in ('f', 'full'):
912 # 2013-04-01, 1.8
913 msg = "".join((
914 "The 'full' option is deprecated in favor of 'reduced'.\n",
915 "For backward compatibility let mode default."))
916 warnings.warn(msg, DeprecationWarning, stacklevel=3)
917 mode = 'reduced'
918 elif mode in ('e', 'economic'):
919 # 2013-04-01, 1.8
920 msg = "The 'economic' option is deprecated."
921 warnings.warn(msg, DeprecationWarning, stacklevel=3)
922 mode = 'economic'
923 else:
924 raise ValueError(f"Unrecognized mode '{mode}'")
926 a, wrap = _makearray(a)
927 _assert_stacked_2d(a)
928 m, n = a.shape[-2:]
929 t, result_t = _commonType(a)
930 a = a.astype(t, copy=True)
931 a = _to_native_byte_order(a)
932 mn = min(m, n)
934 if m <= n:
935 gufunc = _umath_linalg.qr_r_raw_m
936 else:
937 gufunc = _umath_linalg.qr_r_raw_n
939 signature = 'D->D' if isComplexType(t) else 'd->d'
940 extobj = get_linalg_error_extobj(_raise_linalgerror_qr)
941 tau = gufunc(a, signature=signature, extobj=extobj)
943 # handle modes that don't return q
944 if mode == 'r':
945 r = triu(a[..., :mn, :])
946 r = r.astype(result_t, copy=False)
947 return wrap(r)
949 if mode == 'raw':
950 q = transpose(a)
951 q = q.astype(result_t, copy=False)
952 tau = tau.astype(result_t, copy=False)
953 return wrap(q), tau
955 if mode == 'economic':
956 a = a.astype(result_t, copy=False)
957 return wrap(a)
959 # mc is the number of columns in the resulting q
960 # matrix. If the mode is complete then it is
961 # same as number of rows, and if the mode is reduced,
962 # then it is the minimum of number of rows and columns.
963 if mode == 'complete' and m > n:
964 mc = m
965 gufunc = _umath_linalg.qr_complete
966 else:
967 mc = mn
968 gufunc = _umath_linalg.qr_reduced
970 signature = 'DD->D' if isComplexType(t) else 'dd->d'
971 extobj = get_linalg_error_extobj(_raise_linalgerror_qr)
972 q = gufunc(a, tau, signature=signature, extobj=extobj)
973 r = triu(a[..., :mc, :])
975 q = q.astype(result_t, copy=False)
976 r = r.astype(result_t, copy=False)
978 return wrap(q), wrap(r)
980# Eigenvalues
983@array_function_dispatch(_unary_dispatcher)
984def eigvals(a):
985 """
986 Compute the eigenvalues of a general matrix.
988 Main difference between `eigvals` and `eig`: the eigenvectors aren't
989 returned.
991 Parameters
992 ----------
993 a : (..., M, M) array_like
994 A complex- or real-valued matrix whose eigenvalues will be computed.
996 Returns
997 -------
998 w : (..., M,) ndarray
999 The eigenvalues, each repeated according to its multiplicity.
1000 They are not necessarily ordered, nor are they necessarily
1001 real for real matrices.
1003 Raises
1004 ------
1005 LinAlgError
1006 If the eigenvalue computation does not converge.
1008 See Also
1009 --------
1010 eig : eigenvalues and right eigenvectors of general arrays
1011 eigvalsh : eigenvalues of real symmetric or complex Hermitian
1012 (conjugate symmetric) arrays.
1013 eigh : eigenvalues and eigenvectors of real symmetric or complex
1014 Hermitian (conjugate symmetric) arrays.
1015 scipy.linalg.eigvals : Similar function in SciPy.
1017 Notes
1018 -----
1020 .. versionadded:: 1.8.0
1022 Broadcasting rules apply, see the `numpy.linalg` documentation for
1023 details.
1025 This is implemented using the ``_geev`` LAPACK routines which compute
1026 the eigenvalues and eigenvectors of general square arrays.
1028 Examples
1029 --------
1030 Illustration, using the fact that the eigenvalues of a diagonal matrix
1031 are its diagonal elements, that multiplying a matrix on the left
1032 by an orthogonal matrix, `Q`, and on the right by `Q.T` (the transpose
1033 of `Q`), preserves the eigenvalues of the "middle" matrix. In other words,
1034 if `Q` is orthogonal, then ``Q * A * Q.T`` has the same eigenvalues as
1035 ``A``:
1037 >>> from numpy import linalg as LA
1038 >>> x = np.random.random()
1039 >>> Q = np.array([[np.cos(x), -np.sin(x)], [np.sin(x), np.cos(x)]])
1040 >>> LA.norm(Q[0, :]), LA.norm(Q[1, :]), np.dot(Q[0, :],Q[1, :])
1041 (1.0, 1.0, 0.0)
1043 Now multiply a diagonal matrix by ``Q`` on one side and by ``Q.T`` on the other:
1045 >>> D = np.diag((-1,1))
1046 >>> LA.eigvals(D)
1047 array([-1., 1.])
1048 >>> A = np.dot(Q, D)
1049 >>> A = np.dot(A, Q.T)
1050 >>> LA.eigvals(A)
1051 array([ 1., -1.]) # random
1053 """
1054 a, wrap = _makearray(a)
1055 _assert_stacked_2d(a)
1056 _assert_stacked_square(a)
1057 _assert_finite(a)
1058 t, result_t = _commonType(a)
1060 extobj = get_linalg_error_extobj(
1061 _raise_linalgerror_eigenvalues_nonconvergence)
1062 signature = 'D->D' if isComplexType(t) else 'd->D'
1063 w = _umath_linalg.eigvals(a, signature=signature, extobj=extobj)
1065 if not isComplexType(t):
1066 if all(w.imag == 0):
1067 w = w.real
1068 result_t = _realType(result_t)
1069 else:
1070 result_t = _complexType(result_t)
1072 return w.astype(result_t, copy=False)
1075def _eigvalsh_dispatcher(a, UPLO=None):
1076 return (a,)
1079@array_function_dispatch(_eigvalsh_dispatcher)
1080def eigvalsh(a, UPLO='L'):
1081 """
1082 Compute the eigenvalues of a complex Hermitian or real symmetric matrix.
1084 Main difference from eigh: the eigenvectors are not computed.
1086 Parameters
1087 ----------
1088 a : (..., M, M) array_like
1089 A complex- or real-valued matrix whose eigenvalues are to be
1090 computed.
1091 UPLO : {'L', 'U'}, optional
1092 Specifies whether the calculation is done with the lower triangular
1093 part of `a` ('L', default) or the upper triangular part ('U').
1094 Irrespective of this value only the real parts of the diagonal will
1095 be considered in the computation to preserve the notion of a Hermitian
1096 matrix. It therefore follows that the imaginary part of the diagonal
1097 will always be treated as zero.
1099 Returns
1100 -------
1101 w : (..., M,) ndarray
1102 The eigenvalues in ascending order, each repeated according to
1103 its multiplicity.
1105 Raises
1106 ------
1107 LinAlgError
1108 If the eigenvalue computation does not converge.
1110 See Also
1111 --------
1112 eigh : eigenvalues and eigenvectors of real symmetric or complex Hermitian
1113 (conjugate symmetric) arrays.
1114 eigvals : eigenvalues of general real or complex arrays.
1115 eig : eigenvalues and right eigenvectors of general real or complex
1116 arrays.
1117 scipy.linalg.eigvalsh : Similar function in SciPy.
1119 Notes
1120 -----
1122 .. versionadded:: 1.8.0
1124 Broadcasting rules apply, see the `numpy.linalg` documentation for
1125 details.
1127 The eigenvalues are computed using LAPACK routines ``_syevd``, ``_heevd``.
1129 Examples
1130 --------
1131 >>> from numpy import linalg as LA
1132 >>> a = np.array([[1, -2j], [2j, 5]])
1133 >>> LA.eigvalsh(a)
1134 array([ 0.17157288, 5.82842712]) # may vary
1136 >>> # demonstrate the treatment of the imaginary part of the diagonal
1137 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1138 >>> a
1139 array([[5.+2.j, 9.-2.j],
1140 [0.+2.j, 2.-1.j]])
1141 >>> # with UPLO='L' this is numerically equivalent to using LA.eigvals()
1142 >>> # with:
1143 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1144 >>> b
1145 array([[5.+0.j, 0.-2.j],
1146 [0.+2.j, 2.+0.j]])
1147 >>> wa = LA.eigvalsh(a)
1148 >>> wb = LA.eigvals(b)
1149 >>> wa; wb
1150 array([1., 6.])
1151 array([6.+0.j, 1.+0.j])
1153 """
1154 UPLO = UPLO.upper()
1155 if UPLO not in ('L', 'U'):
1156 raise ValueError("UPLO argument must be 'L' or 'U'")
1158 extobj = get_linalg_error_extobj(
1159 _raise_linalgerror_eigenvalues_nonconvergence)
1160 if UPLO == 'L':
1161 gufunc = _umath_linalg.eigvalsh_lo
1162 else:
1163 gufunc = _umath_linalg.eigvalsh_up
1165 a, wrap = _makearray(a)
1166 _assert_stacked_2d(a)
1167 _assert_stacked_square(a)
1168 t, result_t = _commonType(a)
1169 signature = 'D->d' if isComplexType(t) else 'd->d'
1170 w = gufunc(a, signature=signature, extobj=extobj)
1171 return w.astype(_realType(result_t), copy=False)
1173def _convertarray(a):
1174 t, result_t = _commonType(a)
1175 a = _fastCT(a.astype(t))
1176 return a, t, result_t
1179# Eigenvectors
1182@array_function_dispatch(_unary_dispatcher)
1183def eig(a):
1184 """
1185 Compute the eigenvalues and right eigenvectors of a square array.
1187 Parameters
1188 ----------
1189 a : (..., M, M) array
1190 Matrices for which the eigenvalues and right eigenvectors will
1191 be computed
1193 Returns
1194 -------
1195 w : (..., M) array
1196 The eigenvalues, each repeated according to its multiplicity.
1197 The eigenvalues are not necessarily ordered. The resulting
1198 array will be of complex type, unless the imaginary part is
1199 zero in which case it will be cast to a real type. When `a`
1200 is real the resulting eigenvalues will be real (0 imaginary
1201 part) or occur in conjugate pairs
1203 v : (..., M, M) array
1204 The normalized (unit "length") eigenvectors, such that the
1205 column ``v[:,i]`` is the eigenvector corresponding to the
1206 eigenvalue ``w[i]``.
1208 Raises
1209 ------
1210 LinAlgError
1211 If the eigenvalue computation does not converge.
1213 See Also
1214 --------
1215 eigvals : eigenvalues of a non-symmetric array.
1216 eigh : eigenvalues and eigenvectors of a real symmetric or complex
1217 Hermitian (conjugate symmetric) array.
1218 eigvalsh : eigenvalues of a real symmetric or complex Hermitian
1219 (conjugate symmetric) array.
1220 scipy.linalg.eig : Similar function in SciPy that also solves the
1221 generalized eigenvalue problem.
1222 scipy.linalg.schur : Best choice for unitary and other non-Hermitian
1223 normal matrices.
1225 Notes
1226 -----
1228 .. versionadded:: 1.8.0
1230 Broadcasting rules apply, see the `numpy.linalg` documentation for
1231 details.
1233 This is implemented using the ``_geev`` LAPACK routines which compute
1234 the eigenvalues and eigenvectors of general square arrays.
1236 The number `w` is an eigenvalue of `a` if there exists a vector
1237 `v` such that ``a @ v = w * v``. Thus, the arrays `a`, `w`, and
1238 `v` satisfy the equations ``a @ v[:,i] = w[i] * v[:,i]``
1239 for :math:`i \\in \\{0,...,M-1\\}`.
1241 The array `v` of eigenvectors may not be of maximum rank, that is, some
1242 of the columns may be linearly dependent, although round-off error may
1243 obscure that fact. If the eigenvalues are all different, then theoretically
1244 the eigenvectors are linearly independent and `a` can be diagonalized by
1245 a similarity transformation using `v`, i.e, ``inv(v) @ a @ v`` is diagonal.
1247 For non-Hermitian normal matrices the SciPy function `scipy.linalg.schur`
1248 is preferred because the matrix `v` is guaranteed to be unitary, which is
1249 not the case when using `eig`. The Schur factorization produces an
1250 upper triangular matrix rather than a diagonal matrix, but for normal
1251 matrices only the diagonal of the upper triangular matrix is needed, the
1252 rest is roundoff error.
1254 Finally, it is emphasized that `v` consists of the *right* (as in
1255 right-hand side) eigenvectors of `a`. A vector `y` satisfying
1256 ``y.T @ a = z * y.T`` for some number `z` is called a *left*
1257 eigenvector of `a`, and, in general, the left and right eigenvectors
1258 of a matrix are not necessarily the (perhaps conjugate) transposes
1259 of each other.
1261 References
1262 ----------
1263 G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, FL,
1264 Academic Press, Inc., 1980, Various pp.
1266 Examples
1267 --------
1268 >>> from numpy import linalg as LA
1270 (Almost) trivial example with real e-values and e-vectors.
1272 >>> w, v = LA.eig(np.diag((1, 2, 3)))
1273 >>> w; v
1274 array([1., 2., 3.])
1275 array([[1., 0., 0.],
1276 [0., 1., 0.],
1277 [0., 0., 1.]])
1279 Real matrix possessing complex e-values and e-vectors; note that the
1280 e-values are complex conjugates of each other.
1282 >>> w, v = LA.eig(np.array([[1, -1], [1, 1]]))
1283 >>> w; v
1284 array([1.+1.j, 1.-1.j])
1285 array([[0.70710678+0.j , 0.70710678-0.j ],
1286 [0. -0.70710678j, 0. +0.70710678j]])
1288 Complex-valued matrix with real e-values (but complex-valued e-vectors);
1289 note that ``a.conj().T == a``, i.e., `a` is Hermitian.
1291 >>> a = np.array([[1, 1j], [-1j, 1]])
1292 >>> w, v = LA.eig(a)
1293 >>> w; v
1294 array([2.+0.j, 0.+0.j])
1295 array([[ 0. +0.70710678j, 0.70710678+0.j ], # may vary
1296 [ 0.70710678+0.j , -0. +0.70710678j]])
1298 Be careful about round-off error!
1300 >>> a = np.array([[1 + 1e-9, 0], [0, 1 - 1e-9]])
1301 >>> # Theor. e-values are 1 +/- 1e-9
1302 >>> w, v = LA.eig(a)
1303 >>> w; v
1304 array([1., 1.])
1305 array([[1., 0.],
1306 [0., 1.]])
1308 """
1309 a, wrap = _makearray(a)
1310 _assert_stacked_2d(a)
1311 _assert_stacked_square(a)
1312 _assert_finite(a)
1313 t, result_t = _commonType(a)
1315 extobj = get_linalg_error_extobj(
1316 _raise_linalgerror_eigenvalues_nonconvergence)
1317 signature = 'D->DD' if isComplexType(t) else 'd->DD'
1318 w, vt = _umath_linalg.eig(a, signature=signature, extobj=extobj)
1320 if not isComplexType(t) and all(w.imag == 0.0):
1321 w = w.real
1322 vt = vt.real
1323 result_t = _realType(result_t)
1324 else:
1325 result_t = _complexType(result_t)
1327 vt = vt.astype(result_t, copy=False)
1328 return w.astype(result_t, copy=False), wrap(vt)
1331@array_function_dispatch(_eigvalsh_dispatcher)
1332def eigh(a, UPLO='L'):
1333 """
1334 Return the eigenvalues and eigenvectors of a complex Hermitian
1335 (conjugate symmetric) or a real symmetric matrix.
1337 Returns two objects, a 1-D array containing the eigenvalues of `a`, and
1338 a 2-D square array or matrix (depending on the input type) of the
1339 corresponding eigenvectors (in columns).
1341 Parameters
1342 ----------
1343 a : (..., M, M) array
1344 Hermitian or real symmetric matrices whose eigenvalues and
1345 eigenvectors are to be computed.
1346 UPLO : {'L', 'U'}, optional
1347 Specifies whether the calculation is done with the lower triangular
1348 part of `a` ('L', default) or the upper triangular part ('U').
1349 Irrespective of this value only the real parts of the diagonal will
1350 be considered in the computation to preserve the notion of a Hermitian
1351 matrix. It therefore follows that the imaginary part of the diagonal
1352 will always be treated as zero.
1354 Returns
1355 -------
1356 w : (..., M) ndarray
1357 The eigenvalues in ascending order, each repeated according to
1358 its multiplicity.
1359 v : {(..., M, M) ndarray, (..., M, M) matrix}
1360 The column ``v[:, i]`` is the normalized eigenvector corresponding
1361 to the eigenvalue ``w[i]``. Will return a matrix object if `a` is
1362 a matrix object.
1364 Raises
1365 ------
1366 LinAlgError
1367 If the eigenvalue computation does not converge.
1369 See Also
1370 --------
1371 eigvalsh : eigenvalues of real symmetric or complex Hermitian
1372 (conjugate symmetric) arrays.
1373 eig : eigenvalues and right eigenvectors for non-symmetric arrays.
1374 eigvals : eigenvalues of non-symmetric arrays.
1375 scipy.linalg.eigh : Similar function in SciPy (but also solves the
1376 generalized eigenvalue problem).
1378 Notes
1379 -----
1381 .. versionadded:: 1.8.0
1383 Broadcasting rules apply, see the `numpy.linalg` documentation for
1384 details.
1386 The eigenvalues/eigenvectors are computed using LAPACK routines ``_syevd``,
1387 ``_heevd``.
1389 The eigenvalues of real symmetric or complex Hermitian matrices are
1390 always real. [1]_ The array `v` of (column) eigenvectors is unitary
1391 and `a`, `w`, and `v` satisfy the equations
1392 ``dot(a, v[:, i]) = w[i] * v[:, i]``.
1394 References
1395 ----------
1396 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
1397 FL, Academic Press, Inc., 1980, pg. 222.
1399 Examples
1400 --------
1401 >>> from numpy import linalg as LA
1402 >>> a = np.array([[1, -2j], [2j, 5]])
1403 >>> a
1404 array([[ 1.+0.j, -0.-2.j],
1405 [ 0.+2.j, 5.+0.j]])
1406 >>> w, v = LA.eigh(a)
1407 >>> w; v
1408 array([0.17157288, 5.82842712])
1409 array([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1410 [ 0. +0.38268343j, 0. -0.92387953j]])
1412 >>> np.dot(a, v[:, 0]) - w[0] * v[:, 0] # verify 1st e-val/vec pair
1413 array([5.55111512e-17+0.0000000e+00j, 0.00000000e+00+1.2490009e-16j])
1414 >>> np.dot(a, v[:, 1]) - w[1] * v[:, 1] # verify 2nd e-val/vec pair
1415 array([0.+0.j, 0.+0.j])
1417 >>> A = np.matrix(a) # what happens if input is a matrix object
1418 >>> A
1419 matrix([[ 1.+0.j, -0.-2.j],
1420 [ 0.+2.j, 5.+0.j]])
1421 >>> w, v = LA.eigh(A)
1422 >>> w; v
1423 array([0.17157288, 5.82842712])
1424 matrix([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1425 [ 0. +0.38268343j, 0. -0.92387953j]])
1427 >>> # demonstrate the treatment of the imaginary part of the diagonal
1428 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1429 >>> a
1430 array([[5.+2.j, 9.-2.j],
1431 [0.+2.j, 2.-1.j]])
1432 >>> # with UPLO='L' this is numerically equivalent to using LA.eig() with:
1433 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1434 >>> b
1435 array([[5.+0.j, 0.-2.j],
1436 [0.+2.j, 2.+0.j]])
1437 >>> wa, va = LA.eigh(a)
1438 >>> wb, vb = LA.eig(b)
1439 >>> wa; wb
1440 array([1., 6.])
1441 array([6.+0.j, 1.+0.j])
1442 >>> va; vb
1443 array([[-0.4472136 +0.j , -0.89442719+0.j ], # may vary
1444 [ 0. +0.89442719j, 0. -0.4472136j ]])
1445 array([[ 0.89442719+0.j , -0. +0.4472136j],
1446 [-0. +0.4472136j, 0.89442719+0.j ]])
1447 """
1448 UPLO = UPLO.upper()
1449 if UPLO not in ('L', 'U'):
1450 raise ValueError("UPLO argument must be 'L' or 'U'")
1452 a, wrap = _makearray(a)
1453 _assert_stacked_2d(a)
1454 _assert_stacked_square(a)
1455 t, result_t = _commonType(a)
1457 extobj = get_linalg_error_extobj(
1458 _raise_linalgerror_eigenvalues_nonconvergence)
1459 if UPLO == 'L':
1460 gufunc = _umath_linalg.eigh_lo
1461 else:
1462 gufunc = _umath_linalg.eigh_up
1464 signature = 'D->dD' if isComplexType(t) else 'd->dd'
1465 w, vt = gufunc(a, signature=signature, extobj=extobj)
1466 w = w.astype(_realType(result_t), copy=False)
1467 vt = vt.astype(result_t, copy=False)
1468 return w, wrap(vt)
1471# Singular value decomposition
1473def _svd_dispatcher(a, full_matrices=None, compute_uv=None, hermitian=None):
1474 return (a,)
1477@array_function_dispatch(_svd_dispatcher)
1478def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
1479 """
1480 Singular Value Decomposition.
1482 When `a` is a 2D array, and ``full_matrices=False``, then it is
1483 factorized as ``u @ np.diag(s) @ vh = (u * s) @ vh``, where
1484 `u` and the Hermitian transpose of `vh` are 2D arrays with
1485 orthonormal columns and `s` is a 1D array of `a`'s singular
1486 values. When `a` is higher-dimensional, SVD is applied in
1487 stacked mode as explained below.
1489 Parameters
1490 ----------
1491 a : (..., M, N) array_like
1492 A real or complex array with ``a.ndim >= 2``.
1493 full_matrices : bool, optional
1494 If True (default), `u` and `vh` have the shapes ``(..., M, M)`` and
1495 ``(..., N, N)``, respectively. Otherwise, the shapes are
1496 ``(..., M, K)`` and ``(..., K, N)``, respectively, where
1497 ``K = min(M, N)``.
1498 compute_uv : bool, optional
1499 Whether or not to compute `u` and `vh` in addition to `s`. True
1500 by default.
1501 hermitian : bool, optional
1502 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
1503 enabling a more efficient method for finding singular values.
1504 Defaults to False.
1506 .. versionadded:: 1.17.0
1508 Returns
1509 -------
1510 u : { (..., M, M), (..., M, K) } array
1511 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1512 size as those of the input `a`. The size of the last two dimensions
1513 depends on the value of `full_matrices`. Only returned when
1514 `compute_uv` is True.
1515 s : (..., K) array
1516 Vector(s) with the singular values, within each vector sorted in
1517 descending order. The first ``a.ndim - 2`` dimensions have the same
1518 size as those of the input `a`.
1519 vh : { (..., N, N), (..., K, N) } array
1520 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1521 size as those of the input `a`. The size of the last two dimensions
1522 depends on the value of `full_matrices`. Only returned when
1523 `compute_uv` is True.
1525 Raises
1526 ------
1527 LinAlgError
1528 If SVD computation does not converge.
1530 See Also
1531 --------
1532 scipy.linalg.svd : Similar function in SciPy.
1533 scipy.linalg.svdvals : Compute singular values of a matrix.
1535 Notes
1536 -----
1538 .. versionchanged:: 1.8.0
1539 Broadcasting rules apply, see the `numpy.linalg` documentation for
1540 details.
1542 The decomposition is performed using LAPACK routine ``_gesdd``.
1544 SVD is usually described for the factorization of a 2D matrix :math:`A`.
1545 The higher-dimensional case will be discussed below. In the 2D case, SVD is
1546 written as :math:`A = U S V^H`, where :math:`A = a`, :math:`U= u`,
1547 :math:`S= \\mathtt{np.diag}(s)` and :math:`V^H = vh`. The 1D array `s`
1548 contains the singular values of `a` and `u` and `vh` are unitary. The rows
1549 of `vh` are the eigenvectors of :math:`A^H A` and the columns of `u` are
1550 the eigenvectors of :math:`A A^H`. In both cases the corresponding
1551 (possibly non-zero) eigenvalues are given by ``s**2``.
1553 If `a` has more than two dimensions, then broadcasting rules apply, as
1554 explained in :ref:`routines.linalg-broadcasting`. This means that SVD is
1555 working in "stacked" mode: it iterates over all indices of the first
1556 ``a.ndim - 2`` dimensions and for each combination SVD is applied to the
1557 last two indices. The matrix `a` can be reconstructed from the
1558 decomposition with either ``(u * s[..., None, :]) @ vh`` or
1559 ``u @ (s[..., None] * vh)``. (The ``@`` operator can be replaced by the
1560 function ``np.matmul`` for python versions below 3.5.)
1562 If `a` is a ``matrix`` object (as opposed to an ``ndarray``), then so are
1563 all the return values.
1565 Examples
1566 --------
1567 >>> a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6)
1568 >>> b = np.random.randn(2, 7, 8, 3) + 1j*np.random.randn(2, 7, 8, 3)
1570 Reconstruction based on full SVD, 2D case:
1572 >>> u, s, vh = np.linalg.svd(a, full_matrices=True)
1573 >>> u.shape, s.shape, vh.shape
1574 ((9, 9), (6,), (6, 6))
1575 >>> np.allclose(a, np.dot(u[:, :6] * s, vh))
1576 True
1577 >>> smat = np.zeros((9, 6), dtype=complex)
1578 >>> smat[:6, :6] = np.diag(s)
1579 >>> np.allclose(a, np.dot(u, np.dot(smat, vh)))
1580 True
1582 Reconstruction based on reduced SVD, 2D case:
1584 >>> u, s, vh = np.linalg.svd(a, full_matrices=False)
1585 >>> u.shape, s.shape, vh.shape
1586 ((9, 6), (6,), (6, 6))
1587 >>> np.allclose(a, np.dot(u * s, vh))
1588 True
1589 >>> smat = np.diag(s)
1590 >>> np.allclose(a, np.dot(u, np.dot(smat, vh)))
1591 True
1593 Reconstruction based on full SVD, 4D case:
1595 >>> u, s, vh = np.linalg.svd(b, full_matrices=True)
1596 >>> u.shape, s.shape, vh.shape
1597 ((2, 7, 8, 8), (2, 7, 3), (2, 7, 3, 3))
1598 >>> np.allclose(b, np.matmul(u[..., :3] * s[..., None, :], vh))
1599 True
1600 >>> np.allclose(b, np.matmul(u[..., :3], s[..., None] * vh))
1601 True
1603 Reconstruction based on reduced SVD, 4D case:
1605 >>> u, s, vh = np.linalg.svd(b, full_matrices=False)
1606 >>> u.shape, s.shape, vh.shape
1607 ((2, 7, 8, 3), (2, 7, 3), (2, 7, 3, 3))
1608 >>> np.allclose(b, np.matmul(u * s[..., None, :], vh))
1609 True
1610 >>> np.allclose(b, np.matmul(u, s[..., None] * vh))
1611 True
1613 """
1614 import numpy as _nx
1615 a, wrap = _makearray(a)
1617 if hermitian:
1618 # note: lapack svd returns eigenvalues with s ** 2 sorted descending,
1619 # but eig returns s sorted ascending, so we re-order the eigenvalues
1620 # and related arrays to have the correct order
1621 if compute_uv:
1622 s, u = eigh(a)
1623 sgn = sign(s)
1624 s = abs(s)
1625 sidx = argsort(s)[..., ::-1]
1626 sgn = _nx.take_along_axis(sgn, sidx, axis=-1)
1627 s = _nx.take_along_axis(s, sidx, axis=-1)
1628 u = _nx.take_along_axis(u, sidx[..., None, :], axis=-1)
1629 # singular values are unsigned, move the sign into v
1630 vt = transpose(u * sgn[..., None, :]).conjugate()
1631 return wrap(u), s, wrap(vt)
1632 else:
1633 s = eigvalsh(a)
1634 s = s[..., ::-1]
1635 s = abs(s)
1636 return sort(s)[..., ::-1]
1638 _assert_stacked_2d(a)
1639 t, result_t = _commonType(a)
1641 extobj = get_linalg_error_extobj(_raise_linalgerror_svd_nonconvergence)
1643 m, n = a.shape[-2:]
1644 if compute_uv:
1645 if full_matrices:
1646 if m < n:
1647 gufunc = _umath_linalg.svd_m_f
1648 else:
1649 gufunc = _umath_linalg.svd_n_f
1650 else:
1651 if m < n:
1652 gufunc = _umath_linalg.svd_m_s
1653 else:
1654 gufunc = _umath_linalg.svd_n_s
1656 signature = 'D->DdD' if isComplexType(t) else 'd->ddd'
1657 u, s, vh = gufunc(a, signature=signature, extobj=extobj)
1658 u = u.astype(result_t, copy=False)
1659 s = s.astype(_realType(result_t), copy=False)
1660 vh = vh.astype(result_t, copy=False)
1661 return wrap(u), s, wrap(vh)
1662 else:
1663 if m < n:
1664 gufunc = _umath_linalg.svd_m
1665 else:
1666 gufunc = _umath_linalg.svd_n
1668 signature = 'D->d' if isComplexType(t) else 'd->d'
1669 s = gufunc(a, signature=signature, extobj=extobj)
1670 s = s.astype(_realType(result_t), copy=False)
1671 return s
1674def _cond_dispatcher(x, p=None):
1675 return (x,)
1678@array_function_dispatch(_cond_dispatcher)
1679def cond(x, p=None):
1680 """
1681 Compute the condition number of a matrix.
1683 This function is capable of returning the condition number using
1684 one of seven different norms, depending on the value of `p` (see
1685 Parameters below).
1687 Parameters
1688 ----------
1689 x : (..., M, N) array_like
1690 The matrix whose condition number is sought.
1691 p : {None, 1, -1, 2, -2, inf, -inf, 'fro'}, optional
1692 Order of the norm used in the condition number computation:
1694 ===== ============================
1695 p norm for matrices
1696 ===== ============================
1697 None 2-norm, computed directly using the ``SVD``
1698 'fro' Frobenius norm
1699 inf max(sum(abs(x), axis=1))
1700 -inf min(sum(abs(x), axis=1))
1701 1 max(sum(abs(x), axis=0))
1702 -1 min(sum(abs(x), axis=0))
1703 2 2-norm (largest sing. value)
1704 -2 smallest singular value
1705 ===== ============================
1707 inf means the `numpy.inf` object, and the Frobenius norm is
1708 the root-of-sum-of-squares norm.
1710 Returns
1711 -------
1712 c : {float, inf}
1713 The condition number of the matrix. May be infinite.
1715 See Also
1716 --------
1717 numpy.linalg.norm
1719 Notes
1720 -----
1721 The condition number of `x` is defined as the norm of `x` times the
1722 norm of the inverse of `x` [1]_; the norm can be the usual L2-norm
1723 (root-of-sum-of-squares) or one of a number of other matrix norms.
1725 References
1726 ----------
1727 .. [1] G. Strang, *Linear Algebra and Its Applications*, Orlando, FL,
1728 Academic Press, Inc., 1980, pg. 285.
1730 Examples
1731 --------
1732 >>> from numpy import linalg as LA
1733 >>> a = np.array([[1, 0, -1], [0, 1, 0], [1, 0, 1]])
1734 >>> a
1735 array([[ 1, 0, -1],
1736 [ 0, 1, 0],
1737 [ 1, 0, 1]])
1738 >>> LA.cond(a)
1739 1.4142135623730951
1740 >>> LA.cond(a, 'fro')
1741 3.1622776601683795
1742 >>> LA.cond(a, np.inf)
1743 2.0
1744 >>> LA.cond(a, -np.inf)
1745 1.0
1746 >>> LA.cond(a, 1)
1747 2.0
1748 >>> LA.cond(a, -1)
1749 1.0
1750 >>> LA.cond(a, 2)
1751 1.4142135623730951
1752 >>> LA.cond(a, -2)
1753 0.70710678118654746 # may vary
1754 >>> min(LA.svd(a, compute_uv=False))*min(LA.svd(LA.inv(a), compute_uv=False))
1755 0.70710678118654746 # may vary
1757 """
1758 x = asarray(x) # in case we have a matrix
1759 if _is_empty_2d(x):
1760 raise LinAlgError("cond is not defined on empty arrays")
1761 if p is None or p == 2 or p == -2:
1762 s = svd(x, compute_uv=False)
1763 with errstate(all='ignore'):
1764 if p == -2:
1765 r = s[..., -1] / s[..., 0]
1766 else:
1767 r = s[..., 0] / s[..., -1]
1768 else:
1769 # Call inv(x) ignoring errors. The result array will
1770 # contain nans in the entries where inversion failed.
1771 _assert_stacked_2d(x)
1772 _assert_stacked_square(x)
1773 t, result_t = _commonType(x)
1774 signature = 'D->D' if isComplexType(t) else 'd->d'
1775 with errstate(all='ignore'):
1776 invx = _umath_linalg.inv(x, signature=signature)
1777 r = norm(x, p, axis=(-2, -1)) * norm(invx, p, axis=(-2, -1))
1778 r = r.astype(result_t, copy=False)
1780 # Convert nans to infs unless the original array had nan entries
1781 r = asarray(r)
1782 nan_mask = isnan(r)
1783 if nan_mask.any():
1784 nan_mask &= ~isnan(x).any(axis=(-2, -1))
1785 if r.ndim > 0:
1786 r[nan_mask] = Inf
1787 elif nan_mask:
1788 r[()] = Inf
1790 # Convention is to return scalars instead of 0d arrays
1791 if r.ndim == 0:
1792 r = r[()]
1794 return r
1797def _matrix_rank_dispatcher(A, tol=None, hermitian=None):
1798 return (A,)
1801@array_function_dispatch(_matrix_rank_dispatcher)
1802def matrix_rank(A, tol=None, hermitian=False):
1803 """
1804 Return matrix rank of array using SVD method
1806 Rank of the array is the number of singular values of the array that are
1807 greater than `tol`.
1809 .. versionchanged:: 1.14
1810 Can now operate on stacks of matrices
1812 Parameters
1813 ----------
1814 A : {(M,), (..., M, N)} array_like
1815 Input vector or stack of matrices.
1816 tol : (...) array_like, float, optional
1817 Threshold below which SVD values are considered zero. If `tol` is
1818 None, and ``S`` is an array with singular values for `M`, and
1819 ``eps`` is the epsilon value for datatype of ``S``, then `tol` is
1820 set to ``S.max() * max(M, N) * eps``.
1822 .. versionchanged:: 1.14
1823 Broadcasted against the stack of matrices
1824 hermitian : bool, optional
1825 If True, `A` is assumed to be Hermitian (symmetric if real-valued),
1826 enabling a more efficient method for finding singular values.
1827 Defaults to False.
1829 .. versionadded:: 1.14
1831 Returns
1832 -------
1833 rank : (...) array_like
1834 Rank of A.
1836 Notes
1837 -----
1838 The default threshold to detect rank deficiency is a test on the magnitude
1839 of the singular values of `A`. By default, we identify singular values less
1840 than ``S.max() * max(M, N) * eps`` as indicating rank deficiency (with
1841 the symbols defined above). This is the algorithm MATLAB uses [1]. It also
1842 appears in *Numerical recipes* in the discussion of SVD solutions for linear
1843 least squares [2].
1845 This default threshold is designed to detect rank deficiency accounting for
1846 the numerical errors of the SVD computation. Imagine that there is a column
1847 in `A` that is an exact (in floating point) linear combination of other
1848 columns in `A`. Computing the SVD on `A` will not produce a singular value
1849 exactly equal to 0 in general: any difference of the smallest SVD value from
1850 0 will be caused by numerical imprecision in the calculation of the SVD.
1851 Our threshold for small SVD values takes this numerical imprecision into
1852 account, and the default threshold will detect such numerical rank
1853 deficiency. The threshold may declare a matrix `A` rank deficient even if
1854 the linear combination of some columns of `A` is not exactly equal to
1855 another column of `A` but only numerically very close to another column of
1856 `A`.
1858 We chose our default threshold because it is in wide use. Other thresholds
1859 are possible. For example, elsewhere in the 2007 edition of *Numerical
1860 recipes* there is an alternative threshold of ``S.max() *
1861 np.finfo(A.dtype).eps / 2. * np.sqrt(m + n + 1.)``. The authors describe
1862 this threshold as being based on "expected roundoff error" (p 71).
1864 The thresholds above deal with floating point roundoff error in the
1865 calculation of the SVD. However, you may have more information about the
1866 sources of error in `A` that would make you consider other tolerance values
1867 to detect *effective* rank deficiency. The most useful measure of the
1868 tolerance depends on the operations you intend to use on your matrix. For
1869 example, if your data come from uncertain measurements with uncertainties
1870 greater than floating point epsilon, choosing a tolerance near that
1871 uncertainty may be preferable. The tolerance may be absolute if the
1872 uncertainties are absolute rather than relative.
1874 References
1875 ----------
1876 .. [1] MATLAB reference documentation, "Rank"
1877 https://www.mathworks.com/help/techdoc/ref/rank.html
1878 .. [2] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery,
1879 "Numerical Recipes (3rd edition)", Cambridge University Press, 2007,
1880 page 795.
1882 Examples
1883 --------
1884 >>> from numpy.linalg import matrix_rank
1885 >>> matrix_rank(np.eye(4)) # Full rank matrix
1886 4
1887 >>> I=np.eye(4); I[-1,-1] = 0. # rank deficient matrix
1888 >>> matrix_rank(I)
1889 3
1890 >>> matrix_rank(np.ones((4,))) # 1 dimension - rank 1 unless all 0
1891 1
1892 >>> matrix_rank(np.zeros((4,)))
1893 0
1894 """
1895 A = asarray(A)
1896 if A.ndim < 2:
1897 return int(not all(A==0))
1898 S = svd(A, compute_uv=False, hermitian=hermitian)
1899 if tol is None:
1900 tol = S.max(axis=-1, keepdims=True) * max(A.shape[-2:]) * finfo(S.dtype).eps
1901 else:
1902 tol = asarray(tol)[..., newaxis]
1903 return count_nonzero(S > tol, axis=-1)
1906# Generalized inverse
1908def _pinv_dispatcher(a, rcond=None, hermitian=None):
1909 return (a,)
1912@array_function_dispatch(_pinv_dispatcher)
1913def pinv(a, rcond=1e-15, hermitian=False):
1914 """
1915 Compute the (Moore-Penrose) pseudo-inverse of a matrix.
1917 Calculate the generalized inverse of a matrix using its
1918 singular-value decomposition (SVD) and including all
1919 *large* singular values.
1921 .. versionchanged:: 1.14
1922 Can now operate on stacks of matrices
1924 Parameters
1925 ----------
1926 a : (..., M, N) array_like
1927 Matrix or stack of matrices to be pseudo-inverted.
1928 rcond : (...) array_like of float
1929 Cutoff for small singular values.
1930 Singular values less than or equal to
1931 ``rcond * largest_singular_value`` are set to zero.
1932 Broadcasts against the stack of matrices.
1933 hermitian : bool, optional
1934 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
1935 enabling a more efficient method for finding singular values.
1936 Defaults to False.
1938 .. versionadded:: 1.17.0
1940 Returns
1941 -------
1942 B : (..., N, M) ndarray
1943 The pseudo-inverse of `a`. If `a` is a `matrix` instance, then so
1944 is `B`.
1946 Raises
1947 ------
1948 LinAlgError
1949 If the SVD computation does not converge.
1951 See Also
1952 --------
1953 scipy.linalg.pinv : Similar function in SciPy.
1954 scipy.linalg.pinvh : Compute the (Moore-Penrose) pseudo-inverse of a
1955 Hermitian matrix.
1957 Notes
1958 -----
1959 The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
1960 defined as: "the matrix that 'solves' [the least-squares problem]
1961 :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then
1962 :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`.
1964 It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular
1965 value decomposition of A, then
1966 :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are
1967 orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting
1968 of A's so-called singular values, (followed, typically, by
1969 zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix
1970 consisting of the reciprocals of A's singular values
1971 (again, followed by zeros). [1]_
1973 References
1974 ----------
1975 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
1976 FL, Academic Press, Inc., 1980, pp. 139-142.
1978 Examples
1979 --------
1980 The following example checks that ``a * a+ * a == a`` and
1981 ``a+ * a * a+ == a+``:
1983 >>> a = np.random.randn(9, 6)
1984 >>> B = np.linalg.pinv(a)
1985 >>> np.allclose(a, np.dot(a, np.dot(B, a)))
1986 True
1987 >>> np.allclose(B, np.dot(B, np.dot(a, B)))
1988 True
1990 """
1991 a, wrap = _makearray(a)
1992 rcond = asarray(rcond)
1993 if _is_empty_2d(a):
1994 m, n = a.shape[-2:]
1995 res = empty(a.shape[:-2] + (n, m), dtype=a.dtype)
1996 return wrap(res)
1997 a = a.conjugate()
1998 u, s, vt = svd(a, full_matrices=False, hermitian=hermitian)
2000 # discard small singular values
2001 cutoff = rcond[..., newaxis] * amax(s, axis=-1, keepdims=True)
2002 large = s > cutoff
2003 s = divide(1, s, where=large, out=s)
2004 s[~large] = 0
2006 res = matmul(transpose(vt), multiply(s[..., newaxis], transpose(u)))
2007 return wrap(res)
2010# Determinant
2013@array_function_dispatch(_unary_dispatcher)
2014def slogdet(a):
2015 """
2016 Compute the sign and (natural) logarithm of the determinant of an array.
2018 If an array has a very small or very large determinant, then a call to
2019 `det` may overflow or underflow. This routine is more robust against such
2020 issues, because it computes the logarithm of the determinant rather than
2021 the determinant itself.
2023 Parameters
2024 ----------
2025 a : (..., M, M) array_like
2026 Input array, has to be a square 2-D array.
2028 Returns
2029 -------
2030 sign : (...) array_like
2031 A number representing the sign of the determinant. For a real matrix,
2032 this is 1, 0, or -1. For a complex matrix, this is a complex number
2033 with absolute value 1 (i.e., it is on the unit circle), or else 0.
2034 logdet : (...) array_like
2035 The natural log of the absolute value of the determinant.
2037 If the determinant is zero, then `sign` will be 0 and `logdet` will be
2038 -Inf. In all cases, the determinant is equal to ``sign * np.exp(logdet)``.
2040 See Also
2041 --------
2042 det
2044 Notes
2045 -----
2047 .. versionadded:: 1.8.0
2049 Broadcasting rules apply, see the `numpy.linalg` documentation for
2050 details.
2052 .. versionadded:: 1.6.0
2054 The determinant is computed via LU factorization using the LAPACK
2055 routine ``z/dgetrf``.
2058 Examples
2059 --------
2060 The determinant of a 2-D array ``[[a, b], [c, d]]`` is ``ad - bc``:
2062 >>> a = np.array([[1, 2], [3, 4]])
2063 >>> (sign, logdet) = np.linalg.slogdet(a)
2064 >>> (sign, logdet)
2065 (-1, 0.69314718055994529) # may vary
2066 >>> sign * np.exp(logdet)
2067 -2.0
2069 Computing log-determinants for a stack of matrices:
2071 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2072 >>> a.shape
2073 (3, 2, 2)
2074 >>> sign, logdet = np.linalg.slogdet(a)
2075 >>> (sign, logdet)
2076 (array([-1., -1., -1.]), array([ 0.69314718, 1.09861229, 2.07944154]))
2077 >>> sign * np.exp(logdet)
2078 array([-2., -3., -8.])
2080 This routine succeeds where ordinary `det` does not:
2082 >>> np.linalg.det(np.eye(500) * 0.1)
2083 0.0
2084 >>> np.linalg.slogdet(np.eye(500) * 0.1)
2085 (1, -1151.2925464970228)
2087 """
2088 a = asarray(a)
2089 _assert_stacked_2d(a)
2090 _assert_stacked_square(a)
2091 t, result_t = _commonType(a)
2092 real_t = _realType(result_t)
2093 signature = 'D->Dd' if isComplexType(t) else 'd->dd'
2094 sign, logdet = _umath_linalg.slogdet(a, signature=signature)
2095 sign = sign.astype(result_t, copy=False)
2096 logdet = logdet.astype(real_t, copy=False)
2097 return sign, logdet
2100@array_function_dispatch(_unary_dispatcher)
2101def det(a):
2102 """
2103 Compute the determinant of an array.
2105 Parameters
2106 ----------
2107 a : (..., M, M) array_like
2108 Input array to compute determinants for.
2110 Returns
2111 -------
2112 det : (...) array_like
2113 Determinant of `a`.
2115 See Also
2116 --------
2117 slogdet : Another way to represent the determinant, more suitable
2118 for large matrices where underflow/overflow may occur.
2119 scipy.linalg.det : Similar function in SciPy.
2121 Notes
2122 -----
2124 .. versionadded:: 1.8.0
2126 Broadcasting rules apply, see the `numpy.linalg` documentation for
2127 details.
2129 The determinant is computed via LU factorization using the LAPACK
2130 routine ``z/dgetrf``.
2132 Examples
2133 --------
2134 The determinant of a 2-D array [[a, b], [c, d]] is ad - bc:
2136 >>> a = np.array([[1, 2], [3, 4]])
2137 >>> np.linalg.det(a)
2138 -2.0 # may vary
2140 Computing determinants for a stack of matrices:
2142 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2143 >>> a.shape
2144 (3, 2, 2)
2145 >>> np.linalg.det(a)
2146 array([-2., -3., -8.])
2148 """
2149 a = asarray(a)
2150 _assert_stacked_2d(a)
2151 _assert_stacked_square(a)
2152 t, result_t = _commonType(a)
2153 signature = 'D->D' if isComplexType(t) else 'd->d'
2154 r = _umath_linalg.det(a, signature=signature)
2155 r = r.astype(result_t, copy=False)
2156 return r
2159# Linear Least Squares
2161def _lstsq_dispatcher(a, b, rcond=None):
2162 return (a, b)
2165@array_function_dispatch(_lstsq_dispatcher)
2166def lstsq(a, b, rcond="warn"):
2167 r"""
2168 Return the least-squares solution to a linear matrix equation.
2170 Computes the vector `x` that approximately solves the equation
2171 ``a @ x = b``. The equation may be under-, well-, or over-determined
2172 (i.e., the number of linearly independent rows of `a` can be less than,
2173 equal to, or greater than its number of linearly independent columns).
2174 If `a` is square and of full rank, then `x` (but for round-off error)
2175 is the "exact" solution of the equation. Else, `x` minimizes the
2176 Euclidean 2-norm :math:`||b - ax||`. If there are multiple minimizing
2177 solutions, the one with the smallest 2-norm :math:`||x||` is returned.
2179 Parameters
2180 ----------
2181 a : (M, N) array_like
2182 "Coefficient" matrix.
2183 b : {(M,), (M, K)} array_like
2184 Ordinate or "dependent variable" values. If `b` is two-dimensional,
2185 the least-squares solution is calculated for each of the `K` columns
2186 of `b`.
2187 rcond : float, optional
2188 Cut-off ratio for small singular values of `a`.
2189 For the purposes of rank determination, singular values are treated
2190 as zero if they are smaller than `rcond` times the largest singular
2191 value of `a`.
2193 .. versionchanged:: 1.14.0
2194 If not set, a FutureWarning is given. The previous default
2195 of ``-1`` will use the machine precision as `rcond` parameter,
2196 the new default will use the machine precision times `max(M, N)`.
2197 To silence the warning and use the new default, use ``rcond=None``,
2198 to keep using the old behavior, use ``rcond=-1``.
2200 Returns
2201 -------
2202 x : {(N,), (N, K)} ndarray
2203 Least-squares solution. If `b` is two-dimensional,
2204 the solutions are in the `K` columns of `x`.
2205 residuals : {(1,), (K,), (0,)} ndarray
2206 Sums of squared residuals: Squared Euclidean 2-norm for each column in
2207 ``b - a @ x``.
2208 If the rank of `a` is < N or M <= N, this is an empty array.
2209 If `b` is 1-dimensional, this is a (1,) shape array.
2210 Otherwise the shape is (K,).
2211 rank : int
2212 Rank of matrix `a`.
2213 s : (min(M, N),) ndarray
2214 Singular values of `a`.
2216 Raises
2217 ------
2218 LinAlgError
2219 If computation does not converge.
2221 See Also
2222 --------
2223 scipy.linalg.lstsq : Similar function in SciPy.
2225 Notes
2226 -----
2227 If `b` is a matrix, then all array results are returned as matrices.
2229 Examples
2230 --------
2231 Fit a line, ``y = mx + c``, through some noisy data-points:
2233 >>> x = np.array([0, 1, 2, 3])
2234 >>> y = np.array([-1, 0.2, 0.9, 2.1])
2236 By examining the coefficients, we see that the line should have a
2237 gradient of roughly 1 and cut the y-axis at, more or less, -1.
2239 We can rewrite the line equation as ``y = Ap``, where ``A = [[x 1]]``
2240 and ``p = [[m], [c]]``. Now use `lstsq` to solve for `p`:
2242 >>> A = np.vstack([x, np.ones(len(x))]).T
2243 >>> A
2244 array([[ 0., 1.],
2245 [ 1., 1.],
2246 [ 2., 1.],
2247 [ 3., 1.]])
2249 >>> m, c = np.linalg.lstsq(A, y, rcond=None)[0]
2250 >>> m, c
2251 (1.0 -0.95) # may vary
2253 Plot the data along with the fitted line:
2255 >>> import matplotlib.pyplot as plt
2256 >>> _ = plt.plot(x, y, 'o', label='Original data', markersize=10)
2257 >>> _ = plt.plot(x, m*x + c, 'r', label='Fitted line')
2258 >>> _ = plt.legend()
2259 >>> plt.show()
2261 """
2262 a, _ = _makearray(a)
2263 b, wrap = _makearray(b)
2264 is_1d = b.ndim == 1
2265 if is_1d:
2266 b = b[:, newaxis]
2267 _assert_2d(a, b)
2268 m, n = a.shape[-2:]
2269 m2, n_rhs = b.shape[-2:]
2270 if m != m2:
2271 raise LinAlgError('Incompatible dimensions')
2273 t, result_t = _commonType(a, b)
2274 result_real_t = _realType(result_t)
2276 # Determine default rcond value
2277 if rcond == "warn":
2278 # 2017-08-19, 1.14.0
2279 warnings.warn("`rcond` parameter will change to the default of "
2280 "machine precision times ``max(M, N)`` where M and N "
2281 "are the input matrix dimensions.\n"
2282 "To use the future default and silence this warning "
2283 "we advise to pass `rcond=None`, to keep using the old, "
2284 "explicitly pass `rcond=-1`.",
2285 FutureWarning, stacklevel=3)
2286 rcond = -1
2287 if rcond is None:
2288 rcond = finfo(t).eps * max(n, m)
2290 if m <= n:
2291 gufunc = _umath_linalg.lstsq_m
2292 else:
2293 gufunc = _umath_linalg.lstsq_n
2295 signature = 'DDd->Ddid' if isComplexType(t) else 'ddd->ddid'
2296 extobj = get_linalg_error_extobj(_raise_linalgerror_lstsq)
2297 if n_rhs == 0:
2298 # lapack can't handle n_rhs = 0 - so allocate the array one larger in that axis
2299 b = zeros(b.shape[:-2] + (m, n_rhs + 1), dtype=b.dtype)
2300 x, resids, rank, s = gufunc(a, b, rcond, signature=signature, extobj=extobj)
2301 if m == 0:
2302 x[...] = 0
2303 if n_rhs == 0:
2304 # remove the item we added
2305 x = x[..., :n_rhs]
2306 resids = resids[..., :n_rhs]
2308 # remove the axis we added
2309 if is_1d:
2310 x = x.squeeze(axis=-1)
2311 # we probably should squeeze resids too, but we can't
2312 # without breaking compatibility.
2314 # as documented
2315 if rank != n or m <= n:
2316 resids = array([], result_real_t)
2318 # coerce output arrays
2319 s = s.astype(result_real_t, copy=False)
2320 resids = resids.astype(result_real_t, copy=False)
2321 x = x.astype(result_t, copy=True) # Copying lets the memory in r_parts be freed
2322 return wrap(x), wrap(resids), rank, s
2325def _multi_svd_norm(x, row_axis, col_axis, op):
2326 """Compute a function of the singular values of the 2-D matrices in `x`.
2328 This is a private utility function used by `numpy.linalg.norm()`.
2330 Parameters
2331 ----------
2332 x : ndarray
2333 row_axis, col_axis : int
2334 The axes of `x` that hold the 2-D matrices.
2335 op : callable
2336 This should be either numpy.amin or `numpy.amax` or `numpy.sum`.
2338 Returns
2339 -------
2340 result : float or ndarray
2341 If `x` is 2-D, the return values is a float.
2342 Otherwise, it is an array with ``x.ndim - 2`` dimensions.
2343 The return values are either the minimum or maximum or sum of the
2344 singular values of the matrices, depending on whether `op`
2345 is `numpy.amin` or `numpy.amax` or `numpy.sum`.
2347 """
2348 y = moveaxis(x, (row_axis, col_axis), (-2, -1))
2349 result = op(svd(y, compute_uv=False), axis=-1)
2350 return result
2353def _norm_dispatcher(x, ord=None, axis=None, keepdims=None):
2354 return (x,)
2357@array_function_dispatch(_norm_dispatcher)
2358def norm(x, ord=None, axis=None, keepdims=False):
2359 """
2360 Matrix or vector norm.
2362 This function is able to return one of eight different matrix norms,
2363 or one of an infinite number of vector norms (described below), depending
2364 on the value of the ``ord`` parameter.
2366 Parameters
2367 ----------
2368 x : array_like
2369 Input array. If `axis` is None, `x` must be 1-D or 2-D, unless `ord`
2370 is None. If both `axis` and `ord` are None, the 2-norm of
2371 ``x.ravel`` will be returned.
2372 ord : {non-zero int, inf, -inf, 'fro', 'nuc'}, optional
2373 Order of the norm (see table under ``Notes``). inf means numpy's
2374 `inf` object. The default is None.
2375 axis : {None, int, 2-tuple of ints}, optional.
2376 If `axis` is an integer, it specifies the axis of `x` along which to
2377 compute the vector norms. If `axis` is a 2-tuple, it specifies the
2378 axes that hold 2-D matrices, and the matrix norms of these matrices
2379 are computed. If `axis` is None then either a vector norm (when `x`
2380 is 1-D) or a matrix norm (when `x` is 2-D) is returned. The default
2381 is None.
2383 .. versionadded:: 1.8.0
2385 keepdims : bool, optional
2386 If this is set to True, the axes which are normed over are left in the
2387 result as dimensions with size one. With this option the result will
2388 broadcast correctly against the original `x`.
2390 .. versionadded:: 1.10.0
2392 Returns
2393 -------
2394 n : float or ndarray
2395 Norm of the matrix or vector(s).
2397 See Also
2398 --------
2399 scipy.linalg.norm : Similar function in SciPy.
2401 Notes
2402 -----
2403 For values of ``ord < 1``, the result is, strictly speaking, not a
2404 mathematical 'norm', but it may still be useful for various numerical
2405 purposes.
2407 The following norms can be calculated:
2409 ===== ============================ ==========================
2410 ord norm for matrices norm for vectors
2411 ===== ============================ ==========================
2412 None Frobenius norm 2-norm
2413 'fro' Frobenius norm --
2414 'nuc' nuclear norm --
2415 inf max(sum(abs(x), axis=1)) max(abs(x))
2416 -inf min(sum(abs(x), axis=1)) min(abs(x))
2417 0 -- sum(x != 0)
2418 1 max(sum(abs(x), axis=0)) as below
2419 -1 min(sum(abs(x), axis=0)) as below
2420 2 2-norm (largest sing. value) as below
2421 -2 smallest singular value as below
2422 other -- sum(abs(x)**ord)**(1./ord)
2423 ===== ============================ ==========================
2425 The Frobenius norm is given by [1]_:
2427 :math:`||A||_F = [\\sum_{i,j} abs(a_{i,j})^2]^{1/2}`
2429 The nuclear norm is the sum of the singular values.
2431 Both the Frobenius and nuclear norm orders are only defined for
2432 matrices and raise a ValueError when ``x.ndim != 2``.
2434 References
2435 ----------
2436 .. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*,
2437 Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15
2439 Examples
2440 --------
2441 >>> from numpy import linalg as LA
2442 >>> a = np.arange(9) - 4
2443 >>> a
2444 array([-4, -3, -2, ..., 2, 3, 4])
2445 >>> b = a.reshape((3, 3))
2446 >>> b
2447 array([[-4, -3, -2],
2448 [-1, 0, 1],
2449 [ 2, 3, 4]])
2451 >>> LA.norm(a)
2452 7.745966692414834
2453 >>> LA.norm(b)
2454 7.745966692414834
2455 >>> LA.norm(b, 'fro')
2456 7.745966692414834
2457 >>> LA.norm(a, np.inf)
2458 4.0
2459 >>> LA.norm(b, np.inf)
2460 9.0
2461 >>> LA.norm(a, -np.inf)
2462 0.0
2463 >>> LA.norm(b, -np.inf)
2464 2.0
2466 >>> LA.norm(a, 1)
2467 20.0
2468 >>> LA.norm(b, 1)
2469 7.0
2470 >>> LA.norm(a, -1)
2471 -4.6566128774142013e-010
2472 >>> LA.norm(b, -1)
2473 6.0
2474 >>> LA.norm(a, 2)
2475 7.745966692414834
2476 >>> LA.norm(b, 2)
2477 7.3484692283495345
2479 >>> LA.norm(a, -2)
2480 0.0
2481 >>> LA.norm(b, -2)
2482 1.8570331885190563e-016 # may vary
2483 >>> LA.norm(a, 3)
2484 5.8480354764257312 # may vary
2485 >>> LA.norm(a, -3)
2486 0.0
2488 Using the `axis` argument to compute vector norms:
2490 >>> c = np.array([[ 1, 2, 3],
2491 ... [-1, 1, 4]])
2492 >>> LA.norm(c, axis=0)
2493 array([ 1.41421356, 2.23606798, 5. ])
2494 >>> LA.norm(c, axis=1)
2495 array([ 3.74165739, 4.24264069])
2496 >>> LA.norm(c, ord=1, axis=1)
2497 array([ 6., 6.])
2499 Using the `axis` argument to compute matrix norms:
2501 >>> m = np.arange(8).reshape(2,2,2)
2502 >>> LA.norm(m, axis=(1,2))
2503 array([ 3.74165739, 11.22497216])
2504 >>> LA.norm(m[0, :, :]), LA.norm(m[1, :, :])
2505 (3.7416573867739413, 11.224972160321824)
2507 """
2508 x = asarray(x)
2510 if not issubclass(x.dtype.type, (inexact, object_)):
2511 x = x.astype(float)
2513 # Immediately handle some default, simple, fast, and common cases.
2514 if axis is None:
2515 ndim = x.ndim
2516 if ((ord is None) or
2517 (ord in ('f', 'fro') and ndim == 2) or
2518 (ord == 2 and ndim == 1)):
2520 x = x.ravel(order='K')
2521 if isComplexType(x.dtype.type):
2522 x_real = x.real
2523 x_imag = x.imag
2524 sqnorm = x_real.dot(x_real) + x_imag.dot(x_imag)
2525 else:
2526 sqnorm = x.dot(x)
2527 ret = sqrt(sqnorm)
2528 if keepdims:
2529 ret = ret.reshape(ndim*[1])
2530 return ret
2532 # Normalize the `axis` argument to a tuple.
2533 nd = x.ndim
2534 if axis is None:
2535 axis = tuple(range(nd))
2536 elif not isinstance(axis, tuple):
2537 try:
2538 axis = int(axis)
2539 except Exception as e:
2540 raise TypeError("'axis' must be None, an integer or a tuple of integers") from e
2541 axis = (axis,)
2543 if len(axis) == 1:
2544 if ord == Inf:
2545 return abs(x).max(axis=axis, keepdims=keepdims)
2546 elif ord == -Inf:
2547 return abs(x).min(axis=axis, keepdims=keepdims)
2548 elif ord == 0:
2549 # Zero norm
2550 return (x != 0).astype(x.real.dtype).sum(axis=axis, keepdims=keepdims)
2551 elif ord == 1:
2552 # special case for speedup
2553 return add.reduce(abs(x), axis=axis, keepdims=keepdims)
2554 elif ord is None or ord == 2:
2555 # special case for speedup
2556 s = (x.conj() * x).real
2557 return sqrt(add.reduce(s, axis=axis, keepdims=keepdims))
2558 # None of the str-type keywords for ord ('fro', 'nuc')
2559 # are valid for vectors
2560 elif isinstance(ord, str):
2561 raise ValueError(f"Invalid norm order '{ord}' for vectors")
2562 else:
2563 absx = abs(x)
2564 absx **= ord
2565 ret = add.reduce(absx, axis=axis, keepdims=keepdims)
2566 ret **= reciprocal(ord, dtype=ret.dtype)
2567 return ret
2568 elif len(axis) == 2:
2569 row_axis, col_axis = axis
2570 row_axis = normalize_axis_index(row_axis, nd)
2571 col_axis = normalize_axis_index(col_axis, nd)
2572 if row_axis == col_axis:
2573 raise ValueError('Duplicate axes given.')
2574 if ord == 2:
2575 ret = _multi_svd_norm(x, row_axis, col_axis, amax)
2576 elif ord == -2:
2577 ret = _multi_svd_norm(x, row_axis, col_axis, amin)
2578 elif ord == 1:
2579 if col_axis > row_axis:
2580 col_axis -= 1
2581 ret = add.reduce(abs(x), axis=row_axis).max(axis=col_axis)
2582 elif ord == Inf:
2583 if row_axis > col_axis:
2584 row_axis -= 1
2585 ret = add.reduce(abs(x), axis=col_axis).max(axis=row_axis)
2586 elif ord == -1:
2587 if col_axis > row_axis:
2588 col_axis -= 1
2589 ret = add.reduce(abs(x), axis=row_axis).min(axis=col_axis)
2590 elif ord == -Inf:
2591 if row_axis > col_axis:
2592 row_axis -= 1
2593 ret = add.reduce(abs(x), axis=col_axis).min(axis=row_axis)
2594 elif ord in [None, 'fro', 'f']:
2595 ret = sqrt(add.reduce((x.conj() * x).real, axis=axis))
2596 elif ord == 'nuc':
2597 ret = _multi_svd_norm(x, row_axis, col_axis, sum)
2598 else:
2599 raise ValueError("Invalid norm order for matrices.")
2600 if keepdims:
2601 ret_shape = list(x.shape)
2602 ret_shape[axis[0]] = 1
2603 ret_shape[axis[1]] = 1
2604 ret = ret.reshape(ret_shape)
2605 return ret
2606 else:
2607 raise ValueError("Improper number of dimensions to norm.")
2610# multi_dot
2612def _multidot_dispatcher(arrays, *, out=None):
2613 yield from arrays
2614 yield out
2617@array_function_dispatch(_multidot_dispatcher)
2618def multi_dot(arrays, *, out=None):
2619 """
2620 Compute the dot product of two or more arrays in a single function call,
2621 while automatically selecting the fastest evaluation order.
2623 `multi_dot` chains `numpy.dot` and uses optimal parenthesization
2624 of the matrices [1]_ [2]_. Depending on the shapes of the matrices,
2625 this can speed up the multiplication a lot.
2627 If the first argument is 1-D it is treated as a row vector.
2628 If the last argument is 1-D it is treated as a column vector.
2629 The other arguments must be 2-D.
2631 Think of `multi_dot` as::
2633 def multi_dot(arrays): return functools.reduce(np.dot, arrays)
2636 Parameters
2637 ----------
2638 arrays : sequence of array_like
2639 If the first argument is 1-D it is treated as row vector.
2640 If the last argument is 1-D it is treated as column vector.
2641 The other arguments must be 2-D.
2642 out : ndarray, optional
2643 Output argument. This must have the exact kind that would be returned
2644 if it was not used. In particular, it must have the right type, must be
2645 C-contiguous, and its dtype must be the dtype that would be returned
2646 for `dot(a, b)`. This is a performance feature. Therefore, if these
2647 conditions are not met, an exception is raised, instead of attempting
2648 to be flexible.
2650 .. versionadded:: 1.19.0
2652 Returns
2653 -------
2654 output : ndarray
2655 Returns the dot product of the supplied arrays.
2657 See Also
2658 --------
2659 numpy.dot : dot multiplication with two arguments.
2661 References
2662 ----------
2664 .. [1] Cormen, "Introduction to Algorithms", Chapter 15.2, p. 370-378
2665 .. [2] https://en.wikipedia.org/wiki/Matrix_chain_multiplication
2667 Examples
2668 --------
2669 `multi_dot` allows you to write::
2671 >>> from numpy.linalg import multi_dot
2672 >>> # Prepare some data
2673 >>> A = np.random.random((10000, 100))
2674 >>> B = np.random.random((100, 1000))
2675 >>> C = np.random.random((1000, 5))
2676 >>> D = np.random.random((5, 333))
2677 >>> # the actual dot multiplication
2678 >>> _ = multi_dot([A, B, C, D])
2680 instead of::
2682 >>> _ = np.dot(np.dot(np.dot(A, B), C), D)
2683 >>> # or
2684 >>> _ = A.dot(B).dot(C).dot(D)
2686 Notes
2687 -----
2688 The cost for a matrix multiplication can be calculated with the
2689 following function::
2691 def cost(A, B):
2692 return A.shape[0] * A.shape[1] * B.shape[1]
2694 Assume we have three matrices
2695 :math:`A_{10x100}, B_{100x5}, C_{5x50}`.
2697 The costs for the two different parenthesizations are as follows::
2699 cost((AB)C) = 10*100*5 + 10*5*50 = 5000 + 2500 = 7500
2700 cost(A(BC)) = 10*100*50 + 100*5*50 = 50000 + 25000 = 75000
2702 """
2703 n = len(arrays)
2704 # optimization only makes sense for len(arrays) > 2
2705 if n < 2:
2706 raise ValueError("Expecting at least two arrays.")
2707 elif n == 2:
2708 return dot(arrays[0], arrays[1], out=out)
2710 arrays = [asanyarray(a) for a in arrays]
2712 # save original ndim to reshape the result array into the proper form later
2713 ndim_first, ndim_last = arrays[0].ndim, arrays[-1].ndim
2714 # Explicitly convert vectors to 2D arrays to keep the logic of the internal
2715 # _multi_dot_* functions as simple as possible.
2716 if arrays[0].ndim == 1:
2717 arrays[0] = atleast_2d(arrays[0])
2718 if arrays[-1].ndim == 1:
2719 arrays[-1] = atleast_2d(arrays[-1]).T
2720 _assert_2d(*arrays)
2722 # _multi_dot_three is much faster than _multi_dot_matrix_chain_order
2723 if n == 3:
2724 result = _multi_dot_three(arrays[0], arrays[1], arrays[2], out=out)
2725 else:
2726 order = _multi_dot_matrix_chain_order(arrays)
2727 result = _multi_dot(arrays, order, 0, n - 1, out=out)
2729 # return proper shape
2730 if ndim_first == 1 and ndim_last == 1:
2731 return result[0, 0] # scalar
2732 elif ndim_first == 1 or ndim_last == 1:
2733 return result.ravel() # 1-D
2734 else:
2735 return result
2738def _multi_dot_three(A, B, C, out=None):
2739 """
2740 Find the best order for three arrays and do the multiplication.
2742 For three arguments `_multi_dot_three` is approximately 15 times faster
2743 than `_multi_dot_matrix_chain_order`
2745 """
2746 a0, a1b0 = A.shape
2747 b1c0, c1 = C.shape
2748 # cost1 = cost((AB)C) = a0*a1b0*b1c0 + a0*b1c0*c1
2749 cost1 = a0 * b1c0 * (a1b0 + c1)
2750 # cost2 = cost(A(BC)) = a1b0*b1c0*c1 + a0*a1b0*c1
2751 cost2 = a1b0 * c1 * (a0 + b1c0)
2753 if cost1 < cost2:
2754 return dot(dot(A, B), C, out=out)
2755 else:
2756 return dot(A, dot(B, C), out=out)
2759def _multi_dot_matrix_chain_order(arrays, return_costs=False):
2760 """
2761 Return a np.array that encodes the optimal order of mutiplications.
2763 The optimal order array is then used by `_multi_dot()` to do the
2764 multiplication.
2766 Also return the cost matrix if `return_costs` is `True`
2768 The implementation CLOSELY follows Cormen, "Introduction to Algorithms",
2769 Chapter 15.2, p. 370-378. Note that Cormen uses 1-based indices.
2771 cost[i, j] = min([
2772 cost[prefix] + cost[suffix] + cost_mult(prefix, suffix)
2773 for k in range(i, j)])
2775 """
2776 n = len(arrays)
2777 # p stores the dimensions of the matrices
2778 # Example for p: A_{10x100}, B_{100x5}, C_{5x50} --> p = [10, 100, 5, 50]
2779 p = [a.shape[0] for a in arrays] + [arrays[-1].shape[1]]
2780 # m is a matrix of costs of the subproblems
2781 # m[i,j]: min number of scalar multiplications needed to compute A_{i..j}
2782 m = zeros((n, n), dtype=double)
2783 # s is the actual ordering
2784 # s[i, j] is the value of k at which we split the product A_i..A_j
2785 s = empty((n, n), dtype=intp)
2787 for l in range(1, n):
2788 for i in range(n - l):
2789 j = i + l
2790 m[i, j] = Inf
2791 for k in range(i, j):
2792 q = m[i, k] + m[k+1, j] + p[i]*p[k+1]*p[j+1]
2793 if q < m[i, j]:
2794 m[i, j] = q
2795 s[i, j] = k # Note that Cormen uses 1-based index
2797 return (s, m) if return_costs else s
2800def _multi_dot(arrays, order, i, j, out=None):
2801 """Actually do the multiplication with the given order."""
2802 if i == j:
2803 # the initial call with non-None out should never get here
2804 assert out is None
2806 return arrays[i]
2807 else:
2808 return dot(_multi_dot(arrays, order, i, order[i, j]),
2809 _multi_dot(arrays, order, order[i, j] + 1, j),
2810 out=out)