Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/numpy/linalg/linalg.py : 14%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""Lite version of scipy.linalg.
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)
29from numpy.core.multiarray import normalize_axis_index
30from numpy.core.overrides import set_module
31from numpy.core import overrides
32from numpy.lib.twodim_base import triu, eye
33from numpy.linalg import lapack_lite, _umath_linalg
36array_function_dispatch = functools.partial(
37 overrides.array_function_dispatch, module='numpy.linalg')
40fortran_int = intc
43@set_module('numpy.linalg')
44class LinAlgError(Exception):
45 """
46 Generic Python-exception-derived object raised by linalg functions.
48 General purpose exception class, derived from Python's exception.Exception
49 class, programmatically raised in linalg functions when a Linear
50 Algebra-related condition would prevent further correct execution of the
51 function.
53 Parameters
54 ----------
55 None
57 Examples
58 --------
59 >>> from numpy import linalg as LA
60 >>> LA.inv(np.zeros((2,2)))
61 Traceback (most recent call last):
62 File "<stdin>", line 1, in <module>
63 File "...linalg.py", line 350,
64 in inv return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))
65 File "...linalg.py", line 249,
66 in solve
67 raise LinAlgError('Singular matrix')
68 numpy.linalg.LinAlgError: Singular matrix
70 """
73def _determine_error_states():
74 errobj = geterrobj()
75 bufsize = errobj[0]
77 with errstate(invalid='call', over='ignore',
78 divide='ignore', under='ignore'):
79 invalid_call_errmask = geterrobj()[1]
81 return [bufsize, invalid_call_errmask, None]
83# Dealing with errors in _umath_linalg
84_linalg_error_extobj = _determine_error_states()
85del _determine_error_states
87def _raise_linalgerror_singular(err, flag):
88 raise LinAlgError("Singular matrix")
90def _raise_linalgerror_nonposdef(err, flag):
91 raise LinAlgError("Matrix is not positive definite")
93def _raise_linalgerror_eigenvalues_nonconvergence(err, flag):
94 raise LinAlgError("Eigenvalues did not converge")
96def _raise_linalgerror_svd_nonconvergence(err, flag):
97 raise LinAlgError("SVD did not converge")
99def _raise_linalgerror_lstsq(err, flag):
100 raise LinAlgError("SVD did not converge in Linear Least Squares")
102def get_linalg_error_extobj(callback):
103 extobj = list(_linalg_error_extobj) # make a copy
104 extobj[2] = callback
105 return extobj
107def _makearray(a):
108 new = asarray(a)
109 wrap = getattr(a, "__array_prepare__", new.__array_wrap__)
110 return new, wrap
112def isComplexType(t):
113 return issubclass(t, complexfloating)
115_real_types_map = {single : single,
116 double : double,
117 csingle : single,
118 cdouble : double}
120_complex_types_map = {single : csingle,
121 double : cdouble,
122 csingle : csingle,
123 cdouble : cdouble}
125def _realType(t, default=double):
126 return _real_types_map.get(t, default)
128def _complexType(t, default=cdouble):
129 return _complex_types_map.get(t, default)
131def _linalgRealType(t):
132 """Cast the type t to either double or cdouble."""
133 return double
135def _commonType(*arrays):
136 # in lite version, use higher precision (always double or cdouble)
137 result_type = single
138 is_complex = False
139 for a in arrays:
140 if issubclass(a.dtype.type, inexact):
141 if isComplexType(a.dtype.type):
142 is_complex = True
143 rt = _realType(a.dtype.type, default=None)
144 if rt is None:
145 # unsupported inexact scalar
146 raise TypeError("array type %s is unsupported in linalg" %
147 (a.dtype.name,))
148 else:
149 rt = double
150 if rt is double:
151 result_type = double
152 if is_complex:
153 t = cdouble
154 result_type = _complex_types_map[result_type]
155 else:
156 t = double
157 return t, result_type
160# _fastCopyAndTranpose assumes the input is 2D (as all the calls in here are).
162_fastCT = fastCopyAndTranspose
164def _to_native_byte_order(*arrays):
165 ret = []
166 for arr in arrays:
167 if arr.dtype.byteorder not in ('=', '|'):
168 ret.append(asarray(arr, dtype=arr.dtype.newbyteorder('=')))
169 else:
170 ret.append(arr)
171 if len(ret) == 1:
172 return ret[0]
173 else:
174 return ret
176def _fastCopyAndTranspose(type, *arrays):
177 cast_arrays = ()
178 for a in arrays:
179 if a.dtype.type is type:
180 cast_arrays = cast_arrays + (_fastCT(a),)
181 else:
182 cast_arrays = cast_arrays + (_fastCT(a.astype(type)),)
183 if len(cast_arrays) == 1:
184 return cast_arrays[0]
185 else:
186 return cast_arrays
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 a = a.reshape(-1, prod)
304 b = b.ravel()
305 res = wrap(solve(a, b))
306 res.shape = oldshape
307 return res
310def _solve_dispatcher(a, b):
311 return (a, b)
314@array_function_dispatch(_solve_dispatcher)
315def solve(a, b):
316 """
317 Solve a linear matrix equation, or system of linear scalar equations.
319 Computes the "exact" solution, `x`, of the well-determined, i.e., full
320 rank, linear matrix equation `ax = b`.
322 Parameters
323 ----------
324 a : (..., M, M) array_like
325 Coefficient matrix.
326 b : {(..., M,), (..., M, K)}, array_like
327 Ordinate or "dependent variable" values.
329 Returns
330 -------
331 x : {(..., M,), (..., M, K)} ndarray
332 Solution to the system a x = b. Returned shape is identical to `b`.
334 Raises
335 ------
336 LinAlgError
337 If `a` is singular or not square.
339 See Also
340 --------
341 scipy.linalg.solve : Similar function in SciPy.
343 Notes
344 -----
346 .. versionadded:: 1.8.0
348 Broadcasting rules apply, see the `numpy.linalg` documentation for
349 details.
351 The solutions are computed using LAPACK routine ``_gesv``.
353 `a` must be square and of full-rank, i.e., all rows (or, equivalently,
354 columns) must be linearly independent; if either is not true, use
355 `lstsq` for the least-squares best "solution" of the
356 system/equation.
358 References
359 ----------
360 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
361 FL, Academic Press, Inc., 1980, pg. 22.
363 Examples
364 --------
365 Solve the system of equations ``3 * x0 + x1 = 9`` and ``x0 + 2 * x1 = 8``:
367 >>> a = np.array([[3,1], [1,2]])
368 >>> b = np.array([9,8])
369 >>> x = np.linalg.solve(a, b)
370 >>> x
371 array([2., 3.])
373 Check that the solution is correct:
375 >>> np.allclose(np.dot(a, x), b)
376 True
378 """
379 a, _ = _makearray(a)
380 _assert_stacked_2d(a)
381 _assert_stacked_square(a)
382 b, wrap = _makearray(b)
383 t, result_t = _commonType(a, b)
385 # We use the b = (..., M,) logic, only if the number of extra dimensions
386 # match exactly
387 if b.ndim == a.ndim - 1:
388 gufunc = _umath_linalg.solve1
389 else:
390 gufunc = _umath_linalg.solve
392 signature = 'DD->D' if isComplexType(t) else 'dd->d'
393 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
394 r = gufunc(a, b, signature=signature, extobj=extobj)
396 return wrap(r.astype(result_t, copy=False))
399def _tensorinv_dispatcher(a, ind=None):
400 return (a,)
403@array_function_dispatch(_tensorinv_dispatcher)
404def tensorinv(a, ind=2):
405 """
406 Compute the 'inverse' of an N-dimensional array.
408 The result is an inverse for `a` relative to the tensordot operation
409 ``tensordot(a, b, ind)``, i. e., up to floating-point accuracy,
410 ``tensordot(tensorinv(a), a, ind)`` is the "identity" tensor for the
411 tensordot operation.
413 Parameters
414 ----------
415 a : array_like
416 Tensor to 'invert'. Its shape must be 'square', i. e.,
417 ``prod(a.shape[:ind]) == prod(a.shape[ind:])``.
418 ind : int, optional
419 Number of first indices that are involved in the inverse sum.
420 Must be a positive integer, default is 2.
422 Returns
423 -------
424 b : ndarray
425 `a`'s tensordot inverse, shape ``a.shape[ind:] + a.shape[:ind]``.
427 Raises
428 ------
429 LinAlgError
430 If `a` is singular or not 'square' (in the above sense).
432 See Also
433 --------
434 numpy.tensordot, tensorsolve
436 Examples
437 --------
438 >>> a = np.eye(4*6)
439 >>> a.shape = (4, 6, 8, 3)
440 >>> ainv = np.linalg.tensorinv(a, ind=2)
441 >>> ainv.shape
442 (8, 3, 4, 6)
443 >>> b = np.random.randn(4, 6)
444 >>> np.allclose(np.tensordot(ainv, b), np.linalg.tensorsolve(a, b))
445 True
447 >>> a = np.eye(4*6)
448 >>> a.shape = (24, 8, 3)
449 >>> ainv = np.linalg.tensorinv(a, ind=1)
450 >>> ainv.shape
451 (8, 3, 24)
452 >>> b = np.random.randn(24)
453 >>> np.allclose(np.tensordot(ainv, b, 1), np.linalg.tensorsolve(a, b))
454 True
456 """
457 a = asarray(a)
458 oldshape = a.shape
459 prod = 1
460 if ind > 0:
461 invshape = oldshape[ind:] + oldshape[:ind]
462 for k in oldshape[ind:]:
463 prod *= k
464 else:
465 raise ValueError("Invalid ind argument.")
466 a = a.reshape(prod, -1)
467 ia = inv(a)
468 return ia.reshape(*invshape)
471# Matrix inversion
473def _unary_dispatcher(a):
474 return (a,)
477@array_function_dispatch(_unary_dispatcher)
478def inv(a):
479 """
480 Compute the (multiplicative) inverse of a matrix.
482 Given a square matrix `a`, return the matrix `ainv` satisfying
483 ``dot(a, ainv) = dot(ainv, a) = eye(a.shape[0])``.
485 Parameters
486 ----------
487 a : (..., M, M) array_like
488 Matrix to be inverted.
490 Returns
491 -------
492 ainv : (..., M, M) ndarray or matrix
493 (Multiplicative) inverse of the matrix `a`.
495 Raises
496 ------
497 LinAlgError
498 If `a` is not square or inversion fails.
500 See Also
501 --------
502 scipy.linalg.inv : Similar function in SciPy.
504 Notes
505 -----
507 .. versionadded:: 1.8.0
509 Broadcasting rules apply, see the `numpy.linalg` documentation for
510 details.
512 Examples
513 --------
514 >>> from numpy.linalg import inv
515 >>> a = np.array([[1., 2.], [3., 4.]])
516 >>> ainv = inv(a)
517 >>> np.allclose(np.dot(a, ainv), np.eye(2))
518 True
519 >>> np.allclose(np.dot(ainv, a), np.eye(2))
520 True
522 If a is a matrix object, then the return value is a matrix as well:
524 >>> ainv = inv(np.matrix(a))
525 >>> ainv
526 matrix([[-2. , 1. ],
527 [ 1.5, -0.5]])
529 Inverses of several matrices can be computed at once:
531 >>> a = np.array([[[1., 2.], [3., 4.]], [[1, 3], [3, 5]]])
532 >>> inv(a)
533 array([[[-2. , 1. ],
534 [ 1.5 , -0.5 ]],
535 [[-1.25, 0.75],
536 [ 0.75, -0.25]]])
538 """
539 a, wrap = _makearray(a)
540 _assert_stacked_2d(a)
541 _assert_stacked_square(a)
542 t, result_t = _commonType(a)
544 signature = 'D->D' if isComplexType(t) else 'd->d'
545 extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
546 ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
547 return wrap(ainv.astype(result_t, copy=False))
550def _matrix_power_dispatcher(a, n):
551 return (a,)
554@array_function_dispatch(_matrix_power_dispatcher)
555def matrix_power(a, n):
556 """
557 Raise a square matrix to the (integer) power `n`.
559 For positive integers `n`, the power is computed by repeated matrix
560 squarings and matrix multiplications. If ``n == 0``, the identity matrix
561 of the same shape as M is returned. If ``n < 0``, the inverse
562 is computed and then raised to the ``abs(n)``.
564 .. note:: Stacks of object matrices are not currently supported.
566 Parameters
567 ----------
568 a : (..., M, M) array_like
569 Matrix to be "powered".
570 n : int
571 The exponent can be any integer or long integer, positive,
572 negative, or zero.
574 Returns
575 -------
576 a**n : (..., M, M) ndarray or matrix object
577 The return value is the same shape and type as `M`;
578 if the exponent is positive or zero then the type of the
579 elements is the same as those of `M`. If the exponent is
580 negative the elements are floating-point.
582 Raises
583 ------
584 LinAlgError
585 For matrices that are not square or that (for negative powers) cannot
586 be inverted numerically.
588 Examples
589 --------
590 >>> from numpy.linalg import matrix_power
591 >>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit
592 >>> matrix_power(i, 3) # should = -i
593 array([[ 0, -1],
594 [ 1, 0]])
595 >>> matrix_power(i, 0)
596 array([[1, 0],
597 [0, 1]])
598 >>> matrix_power(i, -3) # should = 1/(-i) = i, but w/ f.p. elements
599 array([[ 0., 1.],
600 [-1., 0.]])
602 Somewhat more sophisticated example
604 >>> q = np.zeros((4, 4))
605 >>> q[0:2, 0:2] = -i
606 >>> q[2:4, 2:4] = i
607 >>> q # one of the three quaternion units not equal to 1
608 array([[ 0., -1., 0., 0.],
609 [ 1., 0., 0., 0.],
610 [ 0., 0., 0., 1.],
611 [ 0., 0., -1., 0.]])
612 >>> matrix_power(q, 2) # = -np.eye(4)
613 array([[-1., 0., 0., 0.],
614 [ 0., -1., 0., 0.],
615 [ 0., 0., -1., 0.],
616 [ 0., 0., 0., -1.]])
618 """
619 a = asanyarray(a)
620 _assert_stacked_2d(a)
621 _assert_stacked_square(a)
623 try:
624 n = operator.index(n)
625 except TypeError as e:
626 raise TypeError("exponent must be an integer") from e
628 # Fall back on dot for object arrays. Object arrays are not supported by
629 # the current implementation of matmul using einsum
630 if a.dtype != object:
631 fmatmul = matmul
632 elif a.ndim == 2:
633 fmatmul = dot
634 else:
635 raise NotImplementedError(
636 "matrix_power not supported for stacks of object arrays")
638 if n == 0:
639 a = empty_like(a)
640 a[...] = eye(a.shape[-2], dtype=a.dtype)
641 return a
643 elif n < 0:
644 a = inv(a)
645 n = abs(n)
647 # short-cuts.
648 if n == 1:
649 return a
651 elif n == 2:
652 return fmatmul(a, a)
654 elif n == 3:
655 return fmatmul(fmatmul(a, a), a)
657 # Use binary decomposition to reduce the number of matrix multiplications.
658 # Here, we iterate over the bits of n, from LSB to MSB, raise `a` to
659 # increasing powers of 2, and multiply into the result as needed.
660 z = result = None
661 while n > 0:
662 z = a if z is None else fmatmul(z, z)
663 n, bit = divmod(n, 2)
664 if bit:
665 result = z if result is None else fmatmul(result, z)
667 return result
670# Cholesky decomposition
673@array_function_dispatch(_unary_dispatcher)
674def cholesky(a):
675 """
676 Cholesky decomposition.
678 Return the Cholesky decomposition, `L * L.H`, of the square matrix `a`,
679 where `L` is lower-triangular and .H is the conjugate transpose operator
680 (which is the ordinary transpose if `a` is real-valued). `a` must be
681 Hermitian (symmetric if real-valued) and positive-definite. No
682 checking is performed to verify whether `a` is Hermitian or not.
683 In addition, only the lower-triangular and diagonal elements of `a`
684 are used. Only `L` is actually returned.
686 Parameters
687 ----------
688 a : (..., M, M) array_like
689 Hermitian (symmetric if all elements are real), positive-definite
690 input matrix.
692 Returns
693 -------
694 L : (..., M, M) array_like
695 Upper or lower-triangular Cholesky factor of `a`. Returns a
696 matrix object if `a` is a matrix object.
698 Raises
699 ------
700 LinAlgError
701 If the decomposition fails, for example, if `a` is not
702 positive-definite.
704 See Also
705 --------
706 scipy.linalg.cholesky : Similar function in SciPy.
707 scipy.linalg.cholesky_banded : Cholesky decompose a banded Hermitian
708 positive-definite matrix.
709 scipy.linalg.cho_factor : Cholesky decomposition of a matrix, to use in
710 `scipy.linalg.cho_solve`.
712 Notes
713 -----
715 .. versionadded:: 1.8.0
717 Broadcasting rules apply, see the `numpy.linalg` documentation for
718 details.
720 The Cholesky decomposition is often used as a fast way of solving
722 .. math:: A \\mathbf{x} = \\mathbf{b}
724 (when `A` is both Hermitian/symmetric and positive-definite).
726 First, we solve for :math:`\\mathbf{y}` in
728 .. math:: L \\mathbf{y} = \\mathbf{b},
730 and then for :math:`\\mathbf{x}` in
732 .. math:: L.H \\mathbf{x} = \\mathbf{y}.
734 Examples
735 --------
736 >>> A = np.array([[1,-2j],[2j,5]])
737 >>> A
738 array([[ 1.+0.j, -0.-2.j],
739 [ 0.+2.j, 5.+0.j]])
740 >>> L = np.linalg.cholesky(A)
741 >>> L
742 array([[1.+0.j, 0.+0.j],
743 [0.+2.j, 1.+0.j]])
744 >>> np.dot(L, L.T.conj()) # verify that L * L.H = A
745 array([[1.+0.j, 0.-2.j],
746 [0.+2.j, 5.+0.j]])
747 >>> A = [[1,-2j],[2j,5]] # what happens if A is only array_like?
748 >>> np.linalg.cholesky(A) # an ndarray object is returned
749 array([[1.+0.j, 0.+0.j],
750 [0.+2.j, 1.+0.j]])
751 >>> # But a matrix object is returned if A is a matrix object
752 >>> np.linalg.cholesky(np.matrix(A))
753 matrix([[ 1.+0.j, 0.+0.j],
754 [ 0.+2.j, 1.+0.j]])
756 """
757 extobj = get_linalg_error_extobj(_raise_linalgerror_nonposdef)
758 gufunc = _umath_linalg.cholesky_lo
759 a, wrap = _makearray(a)
760 _assert_stacked_2d(a)
761 _assert_stacked_square(a)
762 t, result_t = _commonType(a)
763 signature = 'D->D' if isComplexType(t) else 'd->d'
764 r = gufunc(a, signature=signature, extobj=extobj)
765 return wrap(r.astype(result_t, copy=False))
768# QR decomposition
770def _qr_dispatcher(a, mode=None):
771 return (a,)
774@array_function_dispatch(_qr_dispatcher)
775def qr(a, mode='reduced'):
776 """
777 Compute the qr factorization of a matrix.
779 Factor the matrix `a` as *qr*, where `q` is orthonormal and `r` is
780 upper-triangular.
782 Parameters
783 ----------
784 a : array_like, shape (M, N)
785 Matrix to be factored.
786 mode : {'reduced', 'complete', 'r', 'raw'}, optional
787 If K = min(M, N), then
789 * 'reduced' : returns q, r with dimensions (M, K), (K, N) (default)
790 * 'complete' : returns q, r with dimensions (M, M), (M, N)
791 * 'r' : returns r only with dimensions (K, N)
792 * 'raw' : returns h, tau with dimensions (N, M), (K,)
794 The options 'reduced', 'complete, and 'raw' are new in numpy 1.8,
795 see the notes for more information. The default is 'reduced', and to
796 maintain backward compatibility with earlier versions of numpy both
797 it and the old default 'full' can be omitted. Note that array h
798 returned in 'raw' mode is transposed for calling Fortran. The
799 'economic' mode is deprecated. The modes 'full' and 'economic' may
800 be passed using only the first letter for backwards compatibility,
801 but all others must be spelled out. See the Notes for more
802 explanation.
805 Returns
806 -------
807 q : ndarray of float or complex, optional
808 A matrix with orthonormal columns. When mode = 'complete' the
809 result is an orthogonal/unitary matrix depending on whether or not
810 a is real/complex. The determinant may be either +/- 1 in that
811 case.
812 r : ndarray of float or complex, optional
813 The upper-triangular matrix.
814 (h, tau) : ndarrays of np.double or np.cdouble, optional
815 The array h contains the Householder reflectors that generate q
816 along with r. The tau array contains scaling factors for the
817 reflectors. In the deprecated 'economic' mode only h is returned.
819 Raises
820 ------
821 LinAlgError
822 If factoring fails.
824 See Also
825 --------
826 scipy.linalg.qr : Similar function in SciPy.
827 scipy.linalg.rq : Compute RQ decomposition of a matrix.
829 Notes
830 -----
831 This is an interface to the LAPACK routines ``dgeqrf``, ``zgeqrf``,
832 ``dorgqr``, and ``zungqr``.
834 For more information on the qr factorization, see for example:
835 https://en.wikipedia.org/wiki/QR_factorization
837 Subclasses of `ndarray` are preserved except for the 'raw' mode. So if
838 `a` is of type `matrix`, all the return values will be matrices too.
840 New 'reduced', 'complete', and 'raw' options for mode were added in
841 NumPy 1.8.0 and the old option 'full' was made an alias of 'reduced'. In
842 addition the options 'full' and 'economic' were deprecated. Because
843 'full' was the previous default and 'reduced' is the new default,
844 backward compatibility can be maintained by letting `mode` default.
845 The 'raw' option was added so that LAPACK routines that can multiply
846 arrays by q using the Householder reflectors can be used. Note that in
847 this case the returned arrays are of type np.double or np.cdouble and
848 the h array is transposed to be FORTRAN compatible. No routines using
849 the 'raw' return are currently exposed by numpy, but some are available
850 in lapack_lite and just await the necessary work.
852 Examples
853 --------
854 >>> a = np.random.randn(9, 6)
855 >>> q, r = np.linalg.qr(a)
856 >>> np.allclose(a, np.dot(q, r)) # a does equal qr
857 True
858 >>> r2 = np.linalg.qr(a, mode='r')
859 >>> np.allclose(r, r2) # mode='r' returns the same r as mode='full'
860 True
862 Example illustrating a common use of `qr`: solving of least squares
863 problems
865 What are the least-squares-best `m` and `y0` in ``y = y0 + mx`` for
866 the following data: {(0,1), (1,0), (1,2), (2,1)}. (Graph the points
867 and you'll see that it should be y0 = 0, m = 1.) The answer is provided
868 by solving the over-determined matrix equation ``Ax = b``, where::
870 A = array([[0, 1], [1, 1], [1, 1], [2, 1]])
871 x = array([[y0], [m]])
872 b = array([[1], [0], [2], [1]])
874 If A = qr such that q is orthonormal (which is always possible via
875 Gram-Schmidt), then ``x = inv(r) * (q.T) * b``. (In numpy practice,
876 however, we simply use `lstsq`.)
878 >>> A = np.array([[0, 1], [1, 1], [1, 1], [2, 1]])
879 >>> A
880 array([[0, 1],
881 [1, 1],
882 [1, 1],
883 [2, 1]])
884 >>> b = np.array([1, 0, 2, 1])
885 >>> q, r = np.linalg.qr(A)
886 >>> p = np.dot(q.T, b)
887 >>> np.dot(np.linalg.inv(r), p)
888 array([ 1.1e-16, 1.0e+00])
890 """
891 if mode not in ('reduced', 'complete', 'r', 'raw'):
892 if mode in ('f', 'full'):
893 # 2013-04-01, 1.8
894 msg = "".join((
895 "The 'full' option is deprecated in favor of 'reduced'.\n",
896 "For backward compatibility let mode default."))
897 warnings.warn(msg, DeprecationWarning, stacklevel=3)
898 mode = 'reduced'
899 elif mode in ('e', 'economic'):
900 # 2013-04-01, 1.8
901 msg = "The 'economic' option is deprecated."
902 warnings.warn(msg, DeprecationWarning, stacklevel=3)
903 mode = 'economic'
904 else:
905 raise ValueError("Unrecognized mode '%s'" % mode)
907 a, wrap = _makearray(a)
908 _assert_2d(a)
909 m, n = a.shape
910 t, result_t = _commonType(a)
911 a = _fastCopyAndTranspose(t, a)
912 a = _to_native_byte_order(a)
913 mn = min(m, n)
914 tau = zeros((mn,), t)
916 if isComplexType(t):
917 lapack_routine = lapack_lite.zgeqrf
918 routine_name = 'zgeqrf'
919 else:
920 lapack_routine = lapack_lite.dgeqrf
921 routine_name = 'dgeqrf'
923 # calculate optimal size of work data 'work'
924 lwork = 1
925 work = zeros((lwork,), t)
926 results = lapack_routine(m, n, a, max(1, m), tau, work, -1, 0)
927 if results['info'] != 0:
928 raise LinAlgError('%s returns %d' % (routine_name, results['info']))
930 # do qr decomposition
931 lwork = max(1, n, int(abs(work[0])))
932 work = zeros((lwork,), t)
933 results = lapack_routine(m, n, a, max(1, m), tau, work, lwork, 0)
934 if results['info'] != 0:
935 raise LinAlgError('%s returns %d' % (routine_name, results['info']))
937 # handle modes that don't return q
938 if mode == 'r':
939 r = _fastCopyAndTranspose(result_t, a[:, :mn])
940 return wrap(triu(r))
942 if mode == 'raw':
943 return a, tau
945 if mode == 'economic':
946 if t != result_t :
947 a = a.astype(result_t, copy=False)
948 return wrap(a.T)
950 # generate q from a
951 if mode == 'complete' and m > n:
952 mc = m
953 q = empty((m, m), t)
954 else:
955 mc = mn
956 q = empty((n, m), t)
957 q[:n] = a
959 if isComplexType(t):
960 lapack_routine = lapack_lite.zungqr
961 routine_name = 'zungqr'
962 else:
963 lapack_routine = lapack_lite.dorgqr
964 routine_name = 'dorgqr'
966 # determine optimal lwork
967 lwork = 1
968 work = zeros((lwork,), t)
969 results = lapack_routine(m, mc, mn, q, max(1, m), tau, work, -1, 0)
970 if results['info'] != 0:
971 raise LinAlgError('%s returns %d' % (routine_name, results['info']))
973 # compute q
974 lwork = max(1, n, int(abs(work[0])))
975 work = zeros((lwork,), t)
976 results = lapack_routine(m, mc, mn, q, max(1, m), tau, work, lwork, 0)
977 if results['info'] != 0:
978 raise LinAlgError('%s returns %d' % (routine_name, results['info']))
980 q = _fastCopyAndTranspose(result_t, q[:mc])
981 r = _fastCopyAndTranspose(result_t, a[:, :mc])
983 return wrap(q), wrap(triu(r))
986# Eigenvalues
989@array_function_dispatch(_unary_dispatcher)
990def eigvals(a):
991 """
992 Compute the eigenvalues of a general matrix.
994 Main difference between `eigvals` and `eig`: the eigenvectors aren't
995 returned.
997 Parameters
998 ----------
999 a : (..., M, M) array_like
1000 A complex- or real-valued matrix whose eigenvalues will be computed.
1002 Returns
1003 -------
1004 w : (..., M,) ndarray
1005 The eigenvalues, each repeated according to its multiplicity.
1006 They are not necessarily ordered, nor are they necessarily
1007 real for real matrices.
1009 Raises
1010 ------
1011 LinAlgError
1012 If the eigenvalue computation does not converge.
1014 See Also
1015 --------
1016 eig : eigenvalues and right eigenvectors of general arrays
1017 eigvalsh : eigenvalues of real symmetric or complex Hermitian
1018 (conjugate symmetric) arrays.
1019 eigh : eigenvalues and eigenvectors of real symmetric or complex
1020 Hermitian (conjugate symmetric) arrays.
1021 scipy.linalg.eigvals : Similar function in SciPy.
1023 Notes
1024 -----
1026 .. versionadded:: 1.8.0
1028 Broadcasting rules apply, see the `numpy.linalg` documentation for
1029 details.
1031 This is implemented using the ``_geev`` LAPACK routines which compute
1032 the eigenvalues and eigenvectors of general square arrays.
1034 Examples
1035 --------
1036 Illustration, using the fact that the eigenvalues of a diagonal matrix
1037 are its diagonal elements, that multiplying a matrix on the left
1038 by an orthogonal matrix, `Q`, and on the right by `Q.T` (the transpose
1039 of `Q`), preserves the eigenvalues of the "middle" matrix. In other words,
1040 if `Q` is orthogonal, then ``Q * A * Q.T`` has the same eigenvalues as
1041 ``A``:
1043 >>> from numpy import linalg as LA
1044 >>> x = np.random.random()
1045 >>> Q = np.array([[np.cos(x), -np.sin(x)], [np.sin(x), np.cos(x)]])
1046 >>> LA.norm(Q[0, :]), LA.norm(Q[1, :]), np.dot(Q[0, :],Q[1, :])
1047 (1.0, 1.0, 0.0)
1049 Now multiply a diagonal matrix by ``Q`` on one side and by ``Q.T`` on the other:
1051 >>> D = np.diag((-1,1))
1052 >>> LA.eigvals(D)
1053 array([-1., 1.])
1054 >>> A = np.dot(Q, D)
1055 >>> A = np.dot(A, Q.T)
1056 >>> LA.eigvals(A)
1057 array([ 1., -1.]) # random
1059 """
1060 a, wrap = _makearray(a)
1061 _assert_stacked_2d(a)
1062 _assert_stacked_square(a)
1063 _assert_finite(a)
1064 t, result_t = _commonType(a)
1066 extobj = get_linalg_error_extobj(
1067 _raise_linalgerror_eigenvalues_nonconvergence)
1068 signature = 'D->D' if isComplexType(t) else 'd->D'
1069 w = _umath_linalg.eigvals(a, signature=signature, extobj=extobj)
1071 if not isComplexType(t):
1072 if all(w.imag == 0):
1073 w = w.real
1074 result_t = _realType(result_t)
1075 else:
1076 result_t = _complexType(result_t)
1078 return w.astype(result_t, copy=False)
1081def _eigvalsh_dispatcher(a, UPLO=None):
1082 return (a,)
1085@array_function_dispatch(_eigvalsh_dispatcher)
1086def eigvalsh(a, UPLO='L'):
1087 """
1088 Compute the eigenvalues of a complex Hermitian or real symmetric matrix.
1090 Main difference from eigh: the eigenvectors are not computed.
1092 Parameters
1093 ----------
1094 a : (..., M, M) array_like
1095 A complex- or real-valued matrix whose eigenvalues are to be
1096 computed.
1097 UPLO : {'L', 'U'}, optional
1098 Specifies whether the calculation is done with the lower triangular
1099 part of `a` ('L', default) or the upper triangular part ('U').
1100 Irrespective of this value only the real parts of the diagonal will
1101 be considered in the computation to preserve the notion of a Hermitian
1102 matrix. It therefore follows that the imaginary part of the diagonal
1103 will always be treated as zero.
1105 Returns
1106 -------
1107 w : (..., M,) ndarray
1108 The eigenvalues in ascending order, each repeated according to
1109 its multiplicity.
1111 Raises
1112 ------
1113 LinAlgError
1114 If the eigenvalue computation does not converge.
1116 See Also
1117 --------
1118 eigh : eigenvalues and eigenvectors of real symmetric or complex Hermitian
1119 (conjugate symmetric) arrays.
1120 eigvals : eigenvalues of general real or complex arrays.
1121 eig : eigenvalues and right eigenvectors of general real or complex
1122 arrays.
1123 scipy.linalg.eigvalsh : Similar function in SciPy.
1125 Notes
1126 -----
1128 .. versionadded:: 1.8.0
1130 Broadcasting rules apply, see the `numpy.linalg` documentation for
1131 details.
1133 The eigenvalues are computed using LAPACK routines ``_syevd``, ``_heevd``.
1135 Examples
1136 --------
1137 >>> from numpy import linalg as LA
1138 >>> a = np.array([[1, -2j], [2j, 5]])
1139 >>> LA.eigvalsh(a)
1140 array([ 0.17157288, 5.82842712]) # may vary
1142 >>> # demonstrate the treatment of the imaginary part of the diagonal
1143 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1144 >>> a
1145 array([[5.+2.j, 9.-2.j],
1146 [0.+2.j, 2.-1.j]])
1147 >>> # with UPLO='L' this is numerically equivalent to using LA.eigvals()
1148 >>> # with:
1149 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1150 >>> b
1151 array([[5.+0.j, 0.-2.j],
1152 [0.+2.j, 2.+0.j]])
1153 >>> wa = LA.eigvalsh(a)
1154 >>> wb = LA.eigvals(b)
1155 >>> wa; wb
1156 array([1., 6.])
1157 array([6.+0.j, 1.+0.j])
1159 """
1160 UPLO = UPLO.upper()
1161 if UPLO not in ('L', 'U'):
1162 raise ValueError("UPLO argument must be 'L' or 'U'")
1164 extobj = get_linalg_error_extobj(
1165 _raise_linalgerror_eigenvalues_nonconvergence)
1166 if UPLO == 'L':
1167 gufunc = _umath_linalg.eigvalsh_lo
1168 else:
1169 gufunc = _umath_linalg.eigvalsh_up
1171 a, wrap = _makearray(a)
1172 _assert_stacked_2d(a)
1173 _assert_stacked_square(a)
1174 t, result_t = _commonType(a)
1175 signature = 'D->d' if isComplexType(t) else 'd->d'
1176 w = gufunc(a, signature=signature, extobj=extobj)
1177 return w.astype(_realType(result_t), copy=False)
1179def _convertarray(a):
1180 t, result_t = _commonType(a)
1181 a = _fastCT(a.astype(t))
1182 return a, t, result_t
1185# Eigenvectors
1188@array_function_dispatch(_unary_dispatcher)
1189def eig(a):
1190 """
1191 Compute the eigenvalues and right eigenvectors of a square array.
1193 Parameters
1194 ----------
1195 a : (..., M, M) array
1196 Matrices for which the eigenvalues and right eigenvectors will
1197 be computed
1199 Returns
1200 -------
1201 w : (..., M) array
1202 The eigenvalues, each repeated according to its multiplicity.
1203 The eigenvalues are not necessarily ordered. The resulting
1204 array will be of complex type, unless the imaginary part is
1205 zero in which case it will be cast to a real type. When `a`
1206 is real the resulting eigenvalues will be real (0 imaginary
1207 part) or occur in conjugate pairs
1209 v : (..., M, M) array
1210 The normalized (unit "length") eigenvectors, such that the
1211 column ``v[:,i]`` is the eigenvector corresponding to the
1212 eigenvalue ``w[i]``.
1214 Raises
1215 ------
1216 LinAlgError
1217 If the eigenvalue computation does not converge.
1219 See Also
1220 --------
1221 eigvals : eigenvalues of a non-symmetric array.
1222 eigh : eigenvalues and eigenvectors of a real symmetric or complex
1223 Hermitian (conjugate symmetric) array.
1224 eigvalsh : eigenvalues of a real symmetric or complex Hermitian
1225 (conjugate symmetric) array.
1226 scipy.linalg.eig : Similar function in SciPy that also solves the
1227 generalized eigenvalue problem.
1228 scipy.linalg.schur : Best choice for unitary and other non-Hermitian
1229 normal matrices.
1231 Notes
1232 -----
1234 .. versionadded:: 1.8.0
1236 Broadcasting rules apply, see the `numpy.linalg` documentation for
1237 details.
1239 This is implemented using the ``_geev`` LAPACK routines which compute
1240 the eigenvalues and eigenvectors of general square arrays.
1242 The number `w` is an eigenvalue of `a` if there exists a vector
1243 `v` such that ``a @ v = w * v``. Thus, the arrays `a`, `w`, and
1244 `v` satisfy the equations ``a @ v[:,i] = w[i] * v[:,i]``
1245 for :math:`i \\in \\{0,...,M-1\\}`.
1247 The array `v` of eigenvectors may not be of maximum rank, that is, some
1248 of the columns may be linearly dependent, although round-off error may
1249 obscure that fact. If the eigenvalues are all different, then theoretically
1250 the eigenvectors are linearly independent and `a` can be diagonalized by
1251 a similarity transformation using `v`, i.e, ``inv(v) @ a @ v`` is diagonal.
1253 For non-Hermitian normal matrices the SciPy function `scipy.linalg.schur`
1254 is preferred because the matrix `v` is guaranteed to be unitary, which is
1255 not the case when using `eig`. The Schur factorization produces an
1256 upper triangular matrix rather than a diagonal matrix, but for normal
1257 matrices only the diagonal of the upper triangular matrix is needed, the
1258 rest is roundoff error.
1260 Finally, it is emphasized that `v` consists of the *right* (as in
1261 right-hand side) eigenvectors of `a`. A vector `y` satisfying
1262 ``y.T @ a = z * y.T`` for some number `z` is called a *left*
1263 eigenvector of `a`, and, in general, the left and right eigenvectors
1264 of a matrix are not necessarily the (perhaps conjugate) transposes
1265 of each other.
1267 References
1268 ----------
1269 G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, FL,
1270 Academic Press, Inc., 1980, Various pp.
1272 Examples
1273 --------
1274 >>> from numpy import linalg as LA
1276 (Almost) trivial example with real e-values and e-vectors.
1278 >>> w, v = LA.eig(np.diag((1, 2, 3)))
1279 >>> w; v
1280 array([1., 2., 3.])
1281 array([[1., 0., 0.],
1282 [0., 1., 0.],
1283 [0., 0., 1.]])
1285 Real matrix possessing complex e-values and e-vectors; note that the
1286 e-values are complex conjugates of each other.
1288 >>> w, v = LA.eig(np.array([[1, -1], [1, 1]]))
1289 >>> w; v
1290 array([1.+1.j, 1.-1.j])
1291 array([[0.70710678+0.j , 0.70710678-0.j ],
1292 [0. -0.70710678j, 0. +0.70710678j]])
1294 Complex-valued matrix with real e-values (but complex-valued e-vectors);
1295 note that ``a.conj().T == a``, i.e., `a` is Hermitian.
1297 >>> a = np.array([[1, 1j], [-1j, 1]])
1298 >>> w, v = LA.eig(a)
1299 >>> w; v
1300 array([2.+0.j, 0.+0.j])
1301 array([[ 0. +0.70710678j, 0.70710678+0.j ], # may vary
1302 [ 0.70710678+0.j , -0. +0.70710678j]])
1304 Be careful about round-off error!
1306 >>> a = np.array([[1 + 1e-9, 0], [0, 1 - 1e-9]])
1307 >>> # Theor. e-values are 1 +/- 1e-9
1308 >>> w, v = LA.eig(a)
1309 >>> w; v
1310 array([1., 1.])
1311 array([[1., 0.],
1312 [0., 1.]])
1314 """
1315 a, wrap = _makearray(a)
1316 _assert_stacked_2d(a)
1317 _assert_stacked_square(a)
1318 _assert_finite(a)
1319 t, result_t = _commonType(a)
1321 extobj = get_linalg_error_extobj(
1322 _raise_linalgerror_eigenvalues_nonconvergence)
1323 signature = 'D->DD' if isComplexType(t) else 'd->DD'
1324 w, vt = _umath_linalg.eig(a, signature=signature, extobj=extobj)
1326 if not isComplexType(t) and all(w.imag == 0.0):
1327 w = w.real
1328 vt = vt.real
1329 result_t = _realType(result_t)
1330 else:
1331 result_t = _complexType(result_t)
1333 vt = vt.astype(result_t, copy=False)
1334 return w.astype(result_t, copy=False), wrap(vt)
1337@array_function_dispatch(_eigvalsh_dispatcher)
1338def eigh(a, UPLO='L'):
1339 """
1340 Return the eigenvalues and eigenvectors of a complex Hermitian
1341 (conjugate symmetric) or a real symmetric matrix.
1343 Returns two objects, a 1-D array containing the eigenvalues of `a`, and
1344 a 2-D square array or matrix (depending on the input type) of the
1345 corresponding eigenvectors (in columns).
1347 Parameters
1348 ----------
1349 a : (..., M, M) array
1350 Hermitian or real symmetric matrices whose eigenvalues and
1351 eigenvectors are to be computed.
1352 UPLO : {'L', 'U'}, optional
1353 Specifies whether the calculation is done with the lower triangular
1354 part of `a` ('L', default) or the upper triangular part ('U').
1355 Irrespective of this value only the real parts of the diagonal will
1356 be considered in the computation to preserve the notion of a Hermitian
1357 matrix. It therefore follows that the imaginary part of the diagonal
1358 will always be treated as zero.
1360 Returns
1361 -------
1362 w : (..., M) ndarray
1363 The eigenvalues in ascending order, each repeated according to
1364 its multiplicity.
1365 v : {(..., M, M) ndarray, (..., M, M) matrix}
1366 The column ``v[:, i]`` is the normalized eigenvector corresponding
1367 to the eigenvalue ``w[i]``. Will return a matrix object if `a` is
1368 a matrix object.
1370 Raises
1371 ------
1372 LinAlgError
1373 If the eigenvalue computation does not converge.
1375 See Also
1376 --------
1377 eigvalsh : eigenvalues of real symmetric or complex Hermitian
1378 (conjugate symmetric) arrays.
1379 eig : eigenvalues and right eigenvectors for non-symmetric arrays.
1380 eigvals : eigenvalues of non-symmetric arrays.
1381 scipy.linalg.eigh : Similar function in SciPy (but also solves the
1382 generalized eigenvalue problem).
1384 Notes
1385 -----
1387 .. versionadded:: 1.8.0
1389 Broadcasting rules apply, see the `numpy.linalg` documentation for
1390 details.
1392 The eigenvalues/eigenvectors are computed using LAPACK routines ``_syevd``,
1393 ``_heevd``.
1395 The eigenvalues of real symmetric or complex Hermitian matrices are
1396 always real. [1]_ The array `v` of (column) eigenvectors is unitary
1397 and `a`, `w`, and `v` satisfy the equations
1398 ``dot(a, v[:, i]) = w[i] * v[:, i]``.
1400 References
1401 ----------
1402 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
1403 FL, Academic Press, Inc., 1980, pg. 222.
1405 Examples
1406 --------
1407 >>> from numpy import linalg as LA
1408 >>> a = np.array([[1, -2j], [2j, 5]])
1409 >>> a
1410 array([[ 1.+0.j, -0.-2.j],
1411 [ 0.+2.j, 5.+0.j]])
1412 >>> w, v = LA.eigh(a)
1413 >>> w; v
1414 array([0.17157288, 5.82842712])
1415 array([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1416 [ 0. +0.38268343j, 0. -0.92387953j]])
1418 >>> np.dot(a, v[:, 0]) - w[0] * v[:, 0] # verify 1st e-val/vec pair
1419 array([5.55111512e-17+0.0000000e+00j, 0.00000000e+00+1.2490009e-16j])
1420 >>> np.dot(a, v[:, 1]) - w[1] * v[:, 1] # verify 2nd e-val/vec pair
1421 array([0.+0.j, 0.+0.j])
1423 >>> A = np.matrix(a) # what happens if input is a matrix object
1424 >>> A
1425 matrix([[ 1.+0.j, -0.-2.j],
1426 [ 0.+2.j, 5.+0.j]])
1427 >>> w, v = LA.eigh(A)
1428 >>> w; v
1429 array([0.17157288, 5.82842712])
1430 matrix([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1431 [ 0. +0.38268343j, 0. -0.92387953j]])
1433 >>> # demonstrate the treatment of the imaginary part of the diagonal
1434 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1435 >>> a
1436 array([[5.+2.j, 9.-2.j],
1437 [0.+2.j, 2.-1.j]])
1438 >>> # with UPLO='L' this is numerically equivalent to using LA.eig() with:
1439 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1440 >>> b
1441 array([[5.+0.j, 0.-2.j],
1442 [0.+2.j, 2.+0.j]])
1443 >>> wa, va = LA.eigh(a)
1444 >>> wb, vb = LA.eig(b)
1445 >>> wa; wb
1446 array([1., 6.])
1447 array([6.+0.j, 1.+0.j])
1448 >>> va; vb
1449 array([[-0.4472136 +0.j , -0.89442719+0.j ], # may vary
1450 [ 0. +0.89442719j, 0. -0.4472136j ]])
1451 array([[ 0.89442719+0.j , -0. +0.4472136j],
1452 [-0. +0.4472136j, 0.89442719+0.j ]])
1453 """
1454 UPLO = UPLO.upper()
1455 if UPLO not in ('L', 'U'):
1456 raise ValueError("UPLO argument must be 'L' or 'U'")
1458 a, wrap = _makearray(a)
1459 _assert_stacked_2d(a)
1460 _assert_stacked_square(a)
1461 t, result_t = _commonType(a)
1463 extobj = get_linalg_error_extobj(
1464 _raise_linalgerror_eigenvalues_nonconvergence)
1465 if UPLO == 'L':
1466 gufunc = _umath_linalg.eigh_lo
1467 else:
1468 gufunc = _umath_linalg.eigh_up
1470 signature = 'D->dD' if isComplexType(t) else 'd->dd'
1471 w, vt = gufunc(a, signature=signature, extobj=extobj)
1472 w = w.astype(_realType(result_t), copy=False)
1473 vt = vt.astype(result_t, copy=False)
1474 return w, wrap(vt)
1477# Singular value decomposition
1479def _svd_dispatcher(a, full_matrices=None, compute_uv=None, hermitian=None):
1480 return (a,)
1483@array_function_dispatch(_svd_dispatcher)
1484def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
1485 """
1486 Singular Value Decomposition.
1488 When `a` is a 2D array, it is factorized as ``u @ np.diag(s) @ vh
1489 = (u * s) @ vh``, where `u` and `vh` are 2D unitary arrays and `s` is a 1D
1490 array of `a`'s singular values. When `a` is higher-dimensional, SVD is
1491 applied in stacked mode as explained below.
1493 Parameters
1494 ----------
1495 a : (..., M, N) array_like
1496 A real or complex array with ``a.ndim >= 2``.
1497 full_matrices : bool, optional
1498 If True (default), `u` and `vh` have the shapes ``(..., M, M)`` and
1499 ``(..., N, N)``, respectively. Otherwise, the shapes are
1500 ``(..., M, K)`` and ``(..., K, N)``, respectively, where
1501 ``K = min(M, N)``.
1502 compute_uv : bool, optional
1503 Whether or not to compute `u` and `vh` in addition to `s`. True
1504 by default.
1505 hermitian : bool, optional
1506 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
1507 enabling a more efficient method for finding singular values.
1508 Defaults to False.
1510 .. versionadded:: 1.17.0
1512 Returns
1513 -------
1514 u : { (..., M, M), (..., M, K) } array
1515 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1516 size as those of the input `a`. The size of the last two dimensions
1517 depends on the value of `full_matrices`. Only returned when
1518 `compute_uv` is True.
1519 s : (..., K) array
1520 Vector(s) with the singular values, within each vector sorted in
1521 descending order. The first ``a.ndim - 2`` dimensions have the same
1522 size as those of the input `a`.
1523 vh : { (..., N, N), (..., K, N) } array
1524 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1525 size as those of the input `a`. The size of the last two dimensions
1526 depends on the value of `full_matrices`. Only returned when
1527 `compute_uv` is True.
1529 Raises
1530 ------
1531 LinAlgError
1532 If SVD computation does not converge.
1534 See Also
1535 --------
1536 scipy.linalg.svd : Similar function in SciPy.
1537 scipy.linalg.svdvals : Compute singular values of a matrix.
1539 Notes
1540 -----
1542 .. versionchanged:: 1.8.0
1543 Broadcasting rules apply, see the `numpy.linalg` documentation for
1544 details.
1546 The decomposition is performed using LAPACK routine ``_gesdd``.
1548 SVD is usually described for the factorization of a 2D matrix :math:`A`.
1549 The higher-dimensional case will be discussed below. In the 2D case, SVD is
1550 written as :math:`A = U S V^H`, where :math:`A = a`, :math:`U= u`,
1551 :math:`S= \\mathtt{np.diag}(s)` and :math:`V^H = vh`. The 1D array `s`
1552 contains the singular values of `a` and `u` and `vh` are unitary. The rows
1553 of `vh` are the eigenvectors of :math:`A^H A` and the columns of `u` are
1554 the eigenvectors of :math:`A A^H`. In both cases the corresponding
1555 (possibly non-zero) eigenvalues are given by ``s**2``.
1557 If `a` has more than two dimensions, then broadcasting rules apply, as
1558 explained in :ref:`routines.linalg-broadcasting`. This means that SVD is
1559 working in "stacked" mode: it iterates over all indices of the first
1560 ``a.ndim - 2`` dimensions and for each combination SVD is applied to the
1561 last two indices. The matrix `a` can be reconstructed from the
1562 decomposition with either ``(u * s[..., None, :]) @ vh`` or
1563 ``u @ (s[..., None] * vh)``. (The ``@`` operator can be replaced by the
1564 function ``np.matmul`` for python versions below 3.5.)
1566 If `a` is a ``matrix`` object (as opposed to an ``ndarray``), then so are
1567 all the return values.
1569 Examples
1570 --------
1571 >>> a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6)
1572 >>> b = np.random.randn(2, 7, 8, 3) + 1j*np.random.randn(2, 7, 8, 3)
1574 Reconstruction based on full SVD, 2D case:
1576 >>> u, s, vh = np.linalg.svd(a, full_matrices=True)
1577 >>> u.shape, s.shape, vh.shape
1578 ((9, 9), (6,), (6, 6))
1579 >>> np.allclose(a, np.dot(u[:, :6] * s, vh))
1580 True
1581 >>> smat = np.zeros((9, 6), dtype=complex)
1582 >>> smat[:6, :6] = np.diag(s)
1583 >>> np.allclose(a, np.dot(u, np.dot(smat, vh)))
1584 True
1586 Reconstruction based on reduced SVD, 2D case:
1588 >>> u, s, vh = np.linalg.svd(a, full_matrices=False)
1589 >>> u.shape, s.shape, vh.shape
1590 ((9, 6), (6,), (6, 6))
1591 >>> np.allclose(a, np.dot(u * s, vh))
1592 True
1593 >>> smat = np.diag(s)
1594 >>> np.allclose(a, np.dot(u, np.dot(smat, vh)))
1595 True
1597 Reconstruction based on full SVD, 4D case:
1599 >>> u, s, vh = np.linalg.svd(b, full_matrices=True)
1600 >>> u.shape, s.shape, vh.shape
1601 ((2, 7, 8, 8), (2, 7, 3), (2, 7, 3, 3))
1602 >>> np.allclose(b, np.matmul(u[..., :3] * s[..., None, :], vh))
1603 True
1604 >>> np.allclose(b, np.matmul(u[..., :3], s[..., None] * vh))
1605 True
1607 Reconstruction based on reduced SVD, 4D case:
1609 >>> u, s, vh = np.linalg.svd(b, full_matrices=False)
1610 >>> u.shape, s.shape, vh.shape
1611 ((2, 7, 8, 3), (2, 7, 3), (2, 7, 3, 3))
1612 >>> np.allclose(b, np.matmul(u * s[..., None, :], vh))
1613 True
1614 >>> np.allclose(b, np.matmul(u, s[..., None] * vh))
1615 True
1617 """
1618 import numpy as _nx
1619 a, wrap = _makearray(a)
1621 if hermitian:
1622 # note: lapack svd returns eigenvalues with s ** 2 sorted descending,
1623 # but eig returns s sorted ascending, so we re-order the eigenvalues
1624 # and related arrays to have the correct order
1625 if compute_uv:
1626 s, u = eigh(a)
1627 sgn = sign(s)
1628 s = abs(s)
1629 sidx = argsort(s)[..., ::-1]
1630 sgn = _nx.take_along_axis(sgn, sidx, axis=-1)
1631 s = _nx.take_along_axis(s, sidx, axis=-1)
1632 u = _nx.take_along_axis(u, sidx[..., None, :], axis=-1)
1633 # singular values are unsigned, move the sign into v
1634 vt = transpose(u * sgn[..., None, :]).conjugate()
1635 return wrap(u), s, wrap(vt)
1636 else:
1637 s = eigvalsh(a)
1638 s = s[..., ::-1]
1639 s = abs(s)
1640 return sort(s)[..., ::-1]
1642 _assert_stacked_2d(a)
1643 t, result_t = _commonType(a)
1645 extobj = get_linalg_error_extobj(_raise_linalgerror_svd_nonconvergence)
1647 m, n = a.shape[-2:]
1648 if compute_uv:
1649 if full_matrices:
1650 if m < n:
1651 gufunc = _umath_linalg.svd_m_f
1652 else:
1653 gufunc = _umath_linalg.svd_n_f
1654 else:
1655 if m < n:
1656 gufunc = _umath_linalg.svd_m_s
1657 else:
1658 gufunc = _umath_linalg.svd_n_s
1660 signature = 'D->DdD' if isComplexType(t) else 'd->ddd'
1661 u, s, vh = gufunc(a, signature=signature, extobj=extobj)
1662 u = u.astype(result_t, copy=False)
1663 s = s.astype(_realType(result_t), copy=False)
1664 vh = vh.astype(result_t, copy=False)
1665 return wrap(u), s, wrap(vh)
1666 else:
1667 if m < n:
1668 gufunc = _umath_linalg.svd_m
1669 else:
1670 gufunc = _umath_linalg.svd_n
1672 signature = 'D->d' if isComplexType(t) else 'd->d'
1673 s = gufunc(a, signature=signature, extobj=extobj)
1674 s = s.astype(_realType(result_t), copy=False)
1675 return s
1678def _cond_dispatcher(x, p=None):
1679 return (x,)
1682@array_function_dispatch(_cond_dispatcher)
1683def cond(x, p=None):
1684 """
1685 Compute the condition number of a matrix.
1687 This function is capable of returning the condition number using
1688 one of seven different norms, depending on the value of `p` (see
1689 Parameters below).
1691 Parameters
1692 ----------
1693 x : (..., M, N) array_like
1694 The matrix whose condition number is sought.
1695 p : {None, 1, -1, 2, -2, inf, -inf, 'fro'}, optional
1696 Order of the norm:
1698 ===== ============================
1699 p norm for matrices
1700 ===== ============================
1701 None 2-norm, computed directly using the ``SVD``
1702 'fro' Frobenius norm
1703 inf max(sum(abs(x), axis=1))
1704 -inf min(sum(abs(x), axis=1))
1705 1 max(sum(abs(x), axis=0))
1706 -1 min(sum(abs(x), axis=0))
1707 2 2-norm (largest sing. value)
1708 -2 smallest singular value
1709 ===== ============================
1711 inf means the numpy.inf object, and the Frobenius norm is
1712 the root-of-sum-of-squares norm.
1714 Returns
1715 -------
1716 c : {float, inf}
1717 The condition number of the matrix. May be infinite.
1719 See Also
1720 --------
1721 numpy.linalg.norm
1723 Notes
1724 -----
1725 The condition number of `x` is defined as the norm of `x` times the
1726 norm of the inverse of `x` [1]_; the norm can be the usual L2-norm
1727 (root-of-sum-of-squares) or one of a number of other matrix norms.
1729 References
1730 ----------
1731 .. [1] G. Strang, *Linear Algebra and Its Applications*, Orlando, FL,
1732 Academic Press, Inc., 1980, pg. 285.
1734 Examples
1735 --------
1736 >>> from numpy import linalg as LA
1737 >>> a = np.array([[1, 0, -1], [0, 1, 0], [1, 0, 1]])
1738 >>> a
1739 array([[ 1, 0, -1],
1740 [ 0, 1, 0],
1741 [ 1, 0, 1]])
1742 >>> LA.cond(a)
1743 1.4142135623730951
1744 >>> LA.cond(a, 'fro')
1745 3.1622776601683795
1746 >>> LA.cond(a, np.inf)
1747 2.0
1748 >>> LA.cond(a, -np.inf)
1749 1.0
1750 >>> LA.cond(a, 1)
1751 2.0
1752 >>> LA.cond(a, -1)
1753 1.0
1754 >>> LA.cond(a, 2)
1755 1.4142135623730951
1756 >>> LA.cond(a, -2)
1757 0.70710678118654746 # may vary
1758 >>> min(LA.svd(a, compute_uv=False))*min(LA.svd(LA.inv(a), compute_uv=False))
1759 0.70710678118654746 # may vary
1761 """
1762 x = asarray(x) # in case we have a matrix
1763 if _is_empty_2d(x):
1764 raise LinAlgError("cond is not defined on empty arrays")
1765 if p is None or p == 2 or p == -2:
1766 s = svd(x, compute_uv=False)
1767 with errstate(all='ignore'):
1768 if p == -2:
1769 r = s[..., -1] / s[..., 0]
1770 else:
1771 r = s[..., 0] / s[..., -1]
1772 else:
1773 # Call inv(x) ignoring errors. The result array will
1774 # contain nans in the entries where inversion failed.
1775 _assert_stacked_2d(x)
1776 _assert_stacked_square(x)
1777 t, result_t = _commonType(x)
1778 signature = 'D->D' if isComplexType(t) else 'd->d'
1779 with errstate(all='ignore'):
1780 invx = _umath_linalg.inv(x, signature=signature)
1781 r = norm(x, p, axis=(-2, -1)) * norm(invx, p, axis=(-2, -1))
1782 r = r.astype(result_t, copy=False)
1784 # Convert nans to infs unless the original array had nan entries
1785 r = asarray(r)
1786 nan_mask = isnan(r)
1787 if nan_mask.any():
1788 nan_mask &= ~isnan(x).any(axis=(-2, -1))
1789 if r.ndim > 0:
1790 r[nan_mask] = Inf
1791 elif nan_mask:
1792 r[()] = Inf
1794 # Convention is to return scalars instead of 0d arrays
1795 if r.ndim == 0:
1796 r = r[()]
1798 return r
1801def _matrix_rank_dispatcher(M, tol=None, hermitian=None):
1802 return (M,)
1805@array_function_dispatch(_matrix_rank_dispatcher)
1806def matrix_rank(M, tol=None, hermitian=False):
1807 """
1808 Return matrix rank of array using SVD method
1810 Rank of the array is the number of singular values of the array that are
1811 greater than `tol`.
1813 .. versionchanged:: 1.14
1814 Can now operate on stacks of matrices
1816 Parameters
1817 ----------
1818 M : {(M,), (..., M, N)} array_like
1819 Input vector or stack of matrices.
1820 tol : (...) array_like, float, optional
1821 Threshold below which SVD values are considered zero. If `tol` is
1822 None, and ``S`` is an array with singular values for `M`, and
1823 ``eps`` is the epsilon value for datatype of ``S``, then `tol` is
1824 set to ``S.max() * max(M.shape) * eps``.
1826 .. versionchanged:: 1.14
1827 Broadcasted against the stack of matrices
1828 hermitian : bool, optional
1829 If True, `M` is assumed to be Hermitian (symmetric if real-valued),
1830 enabling a more efficient method for finding singular values.
1831 Defaults to False.
1833 .. versionadded:: 1.14
1835 Returns
1836 -------
1837 rank : (...) array_like
1838 Rank of M.
1840 Notes
1841 -----
1842 The default threshold to detect rank deficiency is a test on the magnitude
1843 of the singular values of `M`. By default, we identify singular values less
1844 than ``S.max() * max(M.shape) * eps`` as indicating rank deficiency (with
1845 the symbols defined above). This is the algorithm MATLAB uses [1]. It also
1846 appears in *Numerical recipes* in the discussion of SVD solutions for linear
1847 least squares [2].
1849 This default threshold is designed to detect rank deficiency accounting for
1850 the numerical errors of the SVD computation. Imagine that there is a column
1851 in `M` that is an exact (in floating point) linear combination of other
1852 columns in `M`. Computing the SVD on `M` will not produce a singular value
1853 exactly equal to 0 in general: any difference of the smallest SVD value from
1854 0 will be caused by numerical imprecision in the calculation of the SVD.
1855 Our threshold for small SVD values takes this numerical imprecision into
1856 account, and the default threshold will detect such numerical rank
1857 deficiency. The threshold may declare a matrix `M` rank deficient even if
1858 the linear combination of some columns of `M` is not exactly equal to
1859 another column of `M` but only numerically very close to another column of
1860 `M`.
1862 We chose our default threshold because it is in wide use. Other thresholds
1863 are possible. For example, elsewhere in the 2007 edition of *Numerical
1864 recipes* there is an alternative threshold of ``S.max() *
1865 np.finfo(M.dtype).eps / 2. * np.sqrt(m + n + 1.)``. The authors describe
1866 this threshold as being based on "expected roundoff error" (p 71).
1868 The thresholds above deal with floating point roundoff error in the
1869 calculation of the SVD. However, you may have more information about the
1870 sources of error in `M` that would make you consider other tolerance values
1871 to detect *effective* rank deficiency. The most useful measure of the
1872 tolerance depends on the operations you intend to use on your matrix. For
1873 example, if your data come from uncertain measurements with uncertainties
1874 greater than floating point epsilon, choosing a tolerance near that
1875 uncertainty may be preferable. The tolerance may be absolute if the
1876 uncertainties are absolute rather than relative.
1878 References
1879 ----------
1880 .. [1] MATLAB reference documention, "Rank"
1881 https://www.mathworks.com/help/techdoc/ref/rank.html
1882 .. [2] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery,
1883 "Numerical Recipes (3rd edition)", Cambridge University Press, 2007,
1884 page 795.
1886 Examples
1887 --------
1888 >>> from numpy.linalg import matrix_rank
1889 >>> matrix_rank(np.eye(4)) # Full rank matrix
1890 4
1891 >>> I=np.eye(4); I[-1,-1] = 0. # rank deficient matrix
1892 >>> matrix_rank(I)
1893 3
1894 >>> matrix_rank(np.ones((4,))) # 1 dimension - rank 1 unless all 0
1895 1
1896 >>> matrix_rank(np.zeros((4,)))
1897 0
1898 """
1899 M = asarray(M)
1900 if M.ndim < 2:
1901 return int(not all(M==0))
1902 S = svd(M, compute_uv=False, hermitian=hermitian)
1903 if tol is None:
1904 tol = S.max(axis=-1, keepdims=True) * max(M.shape[-2:]) * finfo(S.dtype).eps
1905 else:
1906 tol = asarray(tol)[..., newaxis]
1907 return count_nonzero(S > tol, axis=-1)
1910# Generalized inverse
1912def _pinv_dispatcher(a, rcond=None, hermitian=None):
1913 return (a,)
1916@array_function_dispatch(_pinv_dispatcher)
1917def pinv(a, rcond=1e-15, hermitian=False):
1918 """
1919 Compute the (Moore-Penrose) pseudo-inverse of a matrix.
1921 Calculate the generalized inverse of a matrix using its
1922 singular-value decomposition (SVD) and including all
1923 *large* singular values.
1925 .. versionchanged:: 1.14
1926 Can now operate on stacks of matrices
1928 Parameters
1929 ----------
1930 a : (..., M, N) array_like
1931 Matrix or stack of matrices to be pseudo-inverted.
1932 rcond : (...) array_like of float
1933 Cutoff for small singular values.
1934 Singular values less than or equal to
1935 ``rcond * largest_singular_value`` are set to zero.
1936 Broadcasts against the stack of matrices.
1937 hermitian : bool, optional
1938 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
1939 enabling a more efficient method for finding singular values.
1940 Defaults to False.
1942 .. versionadded:: 1.17.0
1944 Returns
1945 -------
1946 B : (..., N, M) ndarray
1947 The pseudo-inverse of `a`. If `a` is a `matrix` instance, then so
1948 is `B`.
1950 Raises
1951 ------
1952 LinAlgError
1953 If the SVD computation does not converge.
1955 See Also
1956 --------
1957 scipy.linalg.pinv : Similar function in SciPy.
1958 scipy.linalg.pinv2 : Similar function in SciPy (SVD-based).
1959 scipy.linalg.pinvh : Compute the (Moore-Penrose) pseudo-inverse of a
1960 Hermitian matrix.
1962 Notes
1963 -----
1964 The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
1965 defined as: "the matrix that 'solves' [the least-squares problem]
1966 :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then
1967 :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`.
1969 It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular
1970 value decomposition of A, then
1971 :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are
1972 orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting
1973 of A's so-called singular values, (followed, typically, by
1974 zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix
1975 consisting of the reciprocals of A's singular values
1976 (again, followed by zeros). [1]_
1978 References
1979 ----------
1980 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
1981 FL, Academic Press, Inc., 1980, pp. 139-142.
1983 Examples
1984 --------
1985 The following example checks that ``a * a+ * a == a`` and
1986 ``a+ * a * a+ == a+``:
1988 >>> a = np.random.randn(9, 6)
1989 >>> B = np.linalg.pinv(a)
1990 >>> np.allclose(a, np.dot(a, np.dot(B, a)))
1991 True
1992 >>> np.allclose(B, np.dot(B, np.dot(a, B)))
1993 True
1995 """
1996 a, wrap = _makearray(a)
1997 rcond = asarray(rcond)
1998 if _is_empty_2d(a):
1999 m, n = a.shape[-2:]
2000 res = empty(a.shape[:-2] + (n, m), dtype=a.dtype)
2001 return wrap(res)
2002 a = a.conjugate()
2003 u, s, vt = svd(a, full_matrices=False, hermitian=hermitian)
2005 # discard small singular values
2006 cutoff = rcond[..., newaxis] * amax(s, axis=-1, keepdims=True)
2007 large = s > cutoff
2008 s = divide(1, s, where=large, out=s)
2009 s[~large] = 0
2011 res = matmul(transpose(vt), multiply(s[..., newaxis], transpose(u)))
2012 return wrap(res)
2015# Determinant
2018@array_function_dispatch(_unary_dispatcher)
2019def slogdet(a):
2020 """
2021 Compute the sign and (natural) logarithm of the determinant of an array.
2023 If an array has a very small or very large determinant, then a call to
2024 `det` may overflow or underflow. This routine is more robust against such
2025 issues, because it computes the logarithm of the determinant rather than
2026 the determinant itself.
2028 Parameters
2029 ----------
2030 a : (..., M, M) array_like
2031 Input array, has to be a square 2-D array.
2033 Returns
2034 -------
2035 sign : (...) array_like
2036 A number representing the sign of the determinant. For a real matrix,
2037 this is 1, 0, or -1. For a complex matrix, this is a complex number
2038 with absolute value 1 (i.e., it is on the unit circle), or else 0.
2039 logdet : (...) array_like
2040 The natural log of the absolute value of the determinant.
2042 If the determinant is zero, then `sign` will be 0 and `logdet` will be
2043 -Inf. In all cases, the determinant is equal to ``sign * np.exp(logdet)``.
2045 See Also
2046 --------
2047 det
2049 Notes
2050 -----
2052 .. versionadded:: 1.8.0
2054 Broadcasting rules apply, see the `numpy.linalg` documentation for
2055 details.
2057 .. versionadded:: 1.6.0
2059 The determinant is computed via LU factorization using the LAPACK
2060 routine ``z/dgetrf``.
2063 Examples
2064 --------
2065 The determinant of a 2-D array ``[[a, b], [c, d]]`` is ``ad - bc``:
2067 >>> a = np.array([[1, 2], [3, 4]])
2068 >>> (sign, logdet) = np.linalg.slogdet(a)
2069 >>> (sign, logdet)
2070 (-1, 0.69314718055994529) # may vary
2071 >>> sign * np.exp(logdet)
2072 -2.0
2074 Computing log-determinants for a stack of matrices:
2076 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2077 >>> a.shape
2078 (3, 2, 2)
2079 >>> sign, logdet = np.linalg.slogdet(a)
2080 >>> (sign, logdet)
2081 (array([-1., -1., -1.]), array([ 0.69314718, 1.09861229, 2.07944154]))
2082 >>> sign * np.exp(logdet)
2083 array([-2., -3., -8.])
2085 This routine succeeds where ordinary `det` does not:
2087 >>> np.linalg.det(np.eye(500) * 0.1)
2088 0.0
2089 >>> np.linalg.slogdet(np.eye(500) * 0.1)
2090 (1, -1151.2925464970228)
2092 """
2093 a = asarray(a)
2094 _assert_stacked_2d(a)
2095 _assert_stacked_square(a)
2096 t, result_t = _commonType(a)
2097 real_t = _realType(result_t)
2098 signature = 'D->Dd' if isComplexType(t) else 'd->dd'
2099 sign, logdet = _umath_linalg.slogdet(a, signature=signature)
2100 sign = sign.astype(result_t, copy=False)
2101 logdet = logdet.astype(real_t, copy=False)
2102 return sign, logdet
2105@array_function_dispatch(_unary_dispatcher)
2106def det(a):
2107 """
2108 Compute the determinant of an array.
2110 Parameters
2111 ----------
2112 a : (..., M, M) array_like
2113 Input array to compute determinants for.
2115 Returns
2116 -------
2117 det : (...) array_like
2118 Determinant of `a`.
2120 See Also
2121 --------
2122 slogdet : Another way to represent the determinant, more suitable
2123 for large matrices where underflow/overflow may occur.
2124 scipy.linalg.det : Similar function in SciPy.
2126 Notes
2127 -----
2129 .. versionadded:: 1.8.0
2131 Broadcasting rules apply, see the `numpy.linalg` documentation for
2132 details.
2134 The determinant is computed via LU factorization using the LAPACK
2135 routine ``z/dgetrf``.
2137 Examples
2138 --------
2139 The determinant of a 2-D array [[a, b], [c, d]] is ad - bc:
2141 >>> a = np.array([[1, 2], [3, 4]])
2142 >>> np.linalg.det(a)
2143 -2.0 # may vary
2145 Computing determinants for a stack of matrices:
2147 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2148 >>> a.shape
2149 (3, 2, 2)
2150 >>> np.linalg.det(a)
2151 array([-2., -3., -8.])
2153 """
2154 a = asarray(a)
2155 _assert_stacked_2d(a)
2156 _assert_stacked_square(a)
2157 t, result_t = _commonType(a)
2158 signature = 'D->D' if isComplexType(t) else 'd->d'
2159 r = _umath_linalg.det(a, signature=signature)
2160 r = r.astype(result_t, copy=False)
2161 return r
2164# Linear Least Squares
2166def _lstsq_dispatcher(a, b, rcond=None):
2167 return (a, b)
2170@array_function_dispatch(_lstsq_dispatcher)
2171def lstsq(a, b, rcond="warn"):
2172 r"""
2173 Return the least-squares solution to a linear matrix equation.
2175 Computes the vector x that approximatively solves the equation
2176 ``a @ x = b``. The equation may be under-, well-, or over-determined
2177 (i.e., the number of linearly independent rows of `a` can be less than,
2178 equal to, or greater than its number of linearly independent columns).
2179 If `a` is square and of full rank, then `x` (but for round-off error)
2180 is the "exact" solution of the equation. Else, `x` minimizes the
2181 Euclidean 2-norm :math:`|| b - a x ||`.
2183 Parameters
2184 ----------
2185 a : (M, N) array_like
2186 "Coefficient" matrix.
2187 b : {(M,), (M, K)} array_like
2188 Ordinate or "dependent variable" values. If `b` is two-dimensional,
2189 the least-squares solution is calculated for each of the `K` columns
2190 of `b`.
2191 rcond : float, optional
2192 Cut-off ratio for small singular values of `a`.
2193 For the purposes of rank determination, singular values are treated
2194 as zero if they are smaller than `rcond` times the largest singular
2195 value of `a`.
2197 .. versionchanged:: 1.14.0
2198 If not set, a FutureWarning is given. The previous default
2199 of ``-1`` will use the machine precision as `rcond` parameter,
2200 the new default will use the machine precision times `max(M, N)`.
2201 To silence the warning and use the new default, use ``rcond=None``,
2202 to keep using the old behavior, use ``rcond=-1``.
2204 Returns
2205 -------
2206 x : {(N,), (N, K)} ndarray
2207 Least-squares solution. If `b` is two-dimensional,
2208 the solutions are in the `K` columns of `x`.
2209 residuals : {(1,), (K,), (0,)} ndarray
2210 Sums of residuals; squared Euclidean 2-norm for each column in
2211 ``b - a*x``.
2212 If the rank of `a` is < N or M <= N, this is an empty array.
2213 If `b` is 1-dimensional, this is a (1,) shape array.
2214 Otherwise the shape is (K,).
2215 rank : int
2216 Rank of matrix `a`.
2217 s : (min(M, N),) ndarray
2218 Singular values of `a`.
2220 Raises
2221 ------
2222 LinAlgError
2223 If computation does not converge.
2225 See Also
2226 --------
2227 scipy.linalg.lstsq : Similar function in SciPy.
2229 Notes
2230 -----
2231 If `b` is a matrix, then all array results are returned as matrices.
2233 Examples
2234 --------
2235 Fit a line, ``y = mx + c``, through some noisy data-points:
2237 >>> x = np.array([0, 1, 2, 3])
2238 >>> y = np.array([-1, 0.2, 0.9, 2.1])
2240 By examining the coefficients, we see that the line should have a
2241 gradient of roughly 1 and cut the y-axis at, more or less, -1.
2243 We can rewrite the line equation as ``y = Ap``, where ``A = [[x 1]]``
2244 and ``p = [[m], [c]]``. Now use `lstsq` to solve for `p`:
2246 >>> A = np.vstack([x, np.ones(len(x))]).T
2247 >>> A
2248 array([[ 0., 1.],
2249 [ 1., 1.],
2250 [ 2., 1.],
2251 [ 3., 1.]])
2253 >>> m, c = np.linalg.lstsq(A, y, rcond=None)[0]
2254 >>> m, c
2255 (1.0 -0.95) # may vary
2257 Plot the data along with the fitted line:
2259 >>> import matplotlib.pyplot as plt
2260 >>> _ = plt.plot(x, y, 'o', label='Original data', markersize=10)
2261 >>> _ = plt.plot(x, m*x + c, 'r', label='Fitted line')
2262 >>> _ = plt.legend()
2263 >>> plt.show()
2265 """
2266 a, _ = _makearray(a)
2267 b, wrap = _makearray(b)
2268 is_1d = b.ndim == 1
2269 if is_1d:
2270 b = b[:, newaxis]
2271 _assert_2d(a, b)
2272 m, n = a.shape[-2:]
2273 m2, n_rhs = b.shape[-2:]
2274 if m != m2:
2275 raise LinAlgError('Incompatible dimensions')
2277 t, result_t = _commonType(a, b)
2278 # FIXME: real_t is unused
2279 real_t = _linalgRealType(t)
2280 result_real_t = _realType(result_t)
2282 # Determine default rcond value
2283 if rcond == "warn":
2284 # 2017-08-19, 1.14.0
2285 warnings.warn("`rcond` parameter will change to the default of "
2286 "machine precision times ``max(M, N)`` where M and N "
2287 "are the input matrix dimensions.\n"
2288 "To use the future default and silence this warning "
2289 "we advise to pass `rcond=None`, to keep using the old, "
2290 "explicitly pass `rcond=-1`.",
2291 FutureWarning, stacklevel=3)
2292 rcond = -1
2293 if rcond is None:
2294 rcond = finfo(t).eps * max(n, m)
2296 if m <= n:
2297 gufunc = _umath_linalg.lstsq_m
2298 else:
2299 gufunc = _umath_linalg.lstsq_n
2301 signature = 'DDd->Ddid' if isComplexType(t) else 'ddd->ddid'
2302 extobj = get_linalg_error_extobj(_raise_linalgerror_lstsq)
2303 if n_rhs == 0:
2304 # lapack can't handle n_rhs = 0 - so allocate the array one larger in that axis
2305 b = zeros(b.shape[:-2] + (m, n_rhs + 1), dtype=b.dtype)
2306 x, resids, rank, s = gufunc(a, b, rcond, signature=signature, extobj=extobj)
2307 if m == 0:
2308 x[...] = 0
2309 if n_rhs == 0:
2310 # remove the item we added
2311 x = x[..., :n_rhs]
2312 resids = resids[..., :n_rhs]
2314 # remove the axis we added
2315 if is_1d:
2316 x = x.squeeze(axis=-1)
2317 # we probably should squeeze resids too, but we can't
2318 # without breaking compatibility.
2320 # as documented
2321 if rank != n or m <= n:
2322 resids = array([], result_real_t)
2324 # coerce output arrays
2325 s = s.astype(result_real_t, copy=False)
2326 resids = resids.astype(result_real_t, copy=False)
2327 x = x.astype(result_t, copy=True) # Copying lets the memory in r_parts be freed
2328 return wrap(x), wrap(resids), rank, s
2331def _multi_svd_norm(x, row_axis, col_axis, op):
2332 """Compute a function of the singular values of the 2-D matrices in `x`.
2334 This is a private utility function used by `numpy.linalg.norm()`.
2336 Parameters
2337 ----------
2338 x : ndarray
2339 row_axis, col_axis : int
2340 The axes of `x` that hold the 2-D matrices.
2341 op : callable
2342 This should be either numpy.amin or `numpy.amax` or `numpy.sum`.
2344 Returns
2345 -------
2346 result : float or ndarray
2347 If `x` is 2-D, the return values is a float.
2348 Otherwise, it is an array with ``x.ndim - 2`` dimensions.
2349 The return values are either the minimum or maximum or sum of the
2350 singular values of the matrices, depending on whether `op`
2351 is `numpy.amin` or `numpy.amax` or `numpy.sum`.
2353 """
2354 y = moveaxis(x, (row_axis, col_axis), (-2, -1))
2355 result = op(svd(y, compute_uv=False), axis=-1)
2356 return result
2359def _norm_dispatcher(x, ord=None, axis=None, keepdims=None):
2360 return (x,)
2363@array_function_dispatch(_norm_dispatcher)
2364def norm(x, ord=None, axis=None, keepdims=False):
2365 """
2366 Matrix or vector norm.
2368 This function is able to return one of eight different matrix norms,
2369 or one of an infinite number of vector norms (described below), depending
2370 on the value of the ``ord`` parameter.
2372 Parameters
2373 ----------
2374 x : array_like
2375 Input array. If `axis` is None, `x` must be 1-D or 2-D, unless `ord`
2376 is None. If both `axis` and `ord` are None, the 2-norm of
2377 ``x.ravel`` will be returned.
2378 ord : {non-zero int, inf, -inf, 'fro', 'nuc'}, optional
2379 Order of the norm (see table under ``Notes``). inf means numpy's
2380 `inf` object. The default is None.
2381 axis : {None, int, 2-tuple of ints}, optional.
2382 If `axis` is an integer, it specifies the axis of `x` along which to
2383 compute the vector norms. If `axis` is a 2-tuple, it specifies the
2384 axes that hold 2-D matrices, and the matrix norms of these matrices
2385 are computed. If `axis` is None then either a vector norm (when `x`
2386 is 1-D) or a matrix norm (when `x` is 2-D) is returned. The default
2387 is None.
2389 .. versionadded:: 1.8.0
2391 keepdims : bool, optional
2392 If this is set to True, the axes which are normed over are left in the
2393 result as dimensions with size one. With this option the result will
2394 broadcast correctly against the original `x`.
2396 .. versionadded:: 1.10.0
2398 Returns
2399 -------
2400 n : float or ndarray
2401 Norm of the matrix or vector(s).
2403 See Also
2404 --------
2405 scipy.linalg.norm : Similar function in SciPy.
2407 Notes
2408 -----
2409 For values of ``ord < 1``, the result is, strictly speaking, not a
2410 mathematical 'norm', but it may still be useful for various numerical
2411 purposes.
2413 The following norms can be calculated:
2415 ===== ============================ ==========================
2416 ord norm for matrices norm for vectors
2417 ===== ============================ ==========================
2418 None Frobenius norm 2-norm
2419 'fro' Frobenius norm --
2420 'nuc' nuclear norm --
2421 inf max(sum(abs(x), axis=1)) max(abs(x))
2422 -inf min(sum(abs(x), axis=1)) min(abs(x))
2423 0 -- sum(x != 0)
2424 1 max(sum(abs(x), axis=0)) as below
2425 -1 min(sum(abs(x), axis=0)) as below
2426 2 2-norm (largest sing. value) as below
2427 -2 smallest singular value as below
2428 other -- sum(abs(x)**ord)**(1./ord)
2429 ===== ============================ ==========================
2431 The Frobenius norm is given by [1]_:
2433 :math:`||A||_F = [\\sum_{i,j} abs(a_{i,j})^2]^{1/2}`
2435 The nuclear norm is the sum of the singular values.
2437 Both the Frobenius and nuclear norm orders are only defined for
2438 matrices and raise a ValueError when ``x.ndim != 2``.
2440 References
2441 ----------
2442 .. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*,
2443 Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15
2445 Examples
2446 --------
2447 >>> from numpy import linalg as LA
2448 >>> a = np.arange(9) - 4
2449 >>> a
2450 array([-4, -3, -2, ..., 2, 3, 4])
2451 >>> b = a.reshape((3, 3))
2452 >>> b
2453 array([[-4, -3, -2],
2454 [-1, 0, 1],
2455 [ 2, 3, 4]])
2457 >>> LA.norm(a)
2458 7.745966692414834
2459 >>> LA.norm(b)
2460 7.745966692414834
2461 >>> LA.norm(b, 'fro')
2462 7.745966692414834
2463 >>> LA.norm(a, np.inf)
2464 4.0
2465 >>> LA.norm(b, np.inf)
2466 9.0
2467 >>> LA.norm(a, -np.inf)
2468 0.0
2469 >>> LA.norm(b, -np.inf)
2470 2.0
2472 >>> LA.norm(a, 1)
2473 20.0
2474 >>> LA.norm(b, 1)
2475 7.0
2476 >>> LA.norm(a, -1)
2477 -4.6566128774142013e-010
2478 >>> LA.norm(b, -1)
2479 6.0
2480 >>> LA.norm(a, 2)
2481 7.745966692414834
2482 >>> LA.norm(b, 2)
2483 7.3484692283495345
2485 >>> LA.norm(a, -2)
2486 0.0
2487 >>> LA.norm(b, -2)
2488 1.8570331885190563e-016 # may vary
2489 >>> LA.norm(a, 3)
2490 5.8480354764257312 # may vary
2491 >>> LA.norm(a, -3)
2492 0.0
2494 Using the `axis` argument to compute vector norms:
2496 >>> c = np.array([[ 1, 2, 3],
2497 ... [-1, 1, 4]])
2498 >>> LA.norm(c, axis=0)
2499 array([ 1.41421356, 2.23606798, 5. ])
2500 >>> LA.norm(c, axis=1)
2501 array([ 3.74165739, 4.24264069])
2502 >>> LA.norm(c, ord=1, axis=1)
2503 array([ 6., 6.])
2505 Using the `axis` argument to compute matrix norms:
2507 >>> m = np.arange(8).reshape(2,2,2)
2508 >>> LA.norm(m, axis=(1,2))
2509 array([ 3.74165739, 11.22497216])
2510 >>> LA.norm(m[0, :, :]), LA.norm(m[1, :, :])
2511 (3.7416573867739413, 11.224972160321824)
2513 """
2514 x = asarray(x)
2516 if not issubclass(x.dtype.type, (inexact, object_)):
2517 x = x.astype(float)
2519 # Immediately handle some default, simple, fast, and common cases.
2520 if axis is None:
2521 ndim = x.ndim
2522 if ((ord is None) or
2523 (ord in ('f', 'fro') and ndim == 2) or
2524 (ord == 2 and ndim == 1)):
2526 x = x.ravel(order='K')
2527 if isComplexType(x.dtype.type):
2528 sqnorm = dot(x.real, x.real) + dot(x.imag, x.imag)
2529 else:
2530 sqnorm = dot(x, x)
2531 ret = sqrt(sqnorm)
2532 if keepdims:
2533 ret = ret.reshape(ndim*[1])
2534 return ret
2536 # Normalize the `axis` argument to a tuple.
2537 nd = x.ndim
2538 if axis is None:
2539 axis = tuple(range(nd))
2540 elif not isinstance(axis, tuple):
2541 try:
2542 axis = int(axis)
2543 except Exception as e:
2544 raise TypeError("'axis' must be None, an integer or a tuple of integers") from e
2545 axis = (axis,)
2547 if len(axis) == 1:
2548 if ord == Inf:
2549 return abs(x).max(axis=axis, keepdims=keepdims)
2550 elif ord == -Inf:
2551 return abs(x).min(axis=axis, keepdims=keepdims)
2552 elif ord == 0:
2553 # Zero norm
2554 return (x != 0).astype(x.real.dtype).sum(axis=axis, keepdims=keepdims)
2555 elif ord == 1:
2556 # special case for speedup
2557 return add.reduce(abs(x), axis=axis, keepdims=keepdims)
2558 elif ord is None or ord == 2:
2559 # special case for speedup
2560 s = (x.conj() * x).real
2561 return sqrt(add.reduce(s, axis=axis, keepdims=keepdims))
2562 # None of the str-type keywords for ord ('fro', 'nuc')
2563 # are valid for vectors
2564 elif isinstance(ord, str):
2565 raise ValueError(f"Invalid norm order '{ord}' for vectors")
2566 else:
2567 absx = abs(x)
2568 absx **= ord
2569 ret = add.reduce(absx, axis=axis, keepdims=keepdims)
2570 ret **= (1 / ord)
2571 return ret
2572 elif len(axis) == 2:
2573 row_axis, col_axis = axis
2574 row_axis = normalize_axis_index(row_axis, nd)
2575 col_axis = normalize_axis_index(col_axis, nd)
2576 if row_axis == col_axis:
2577 raise ValueError('Duplicate axes given.')
2578 if ord == 2:
2579 ret = _multi_svd_norm(x, row_axis, col_axis, amax)
2580 elif ord == -2:
2581 ret = _multi_svd_norm(x, row_axis, col_axis, amin)
2582 elif ord == 1:
2583 if col_axis > row_axis:
2584 col_axis -= 1
2585 ret = add.reduce(abs(x), axis=row_axis).max(axis=col_axis)
2586 elif ord == Inf:
2587 if row_axis > col_axis:
2588 row_axis -= 1
2589 ret = add.reduce(abs(x), axis=col_axis).max(axis=row_axis)
2590 elif ord == -1:
2591 if col_axis > row_axis:
2592 col_axis -= 1
2593 ret = add.reduce(abs(x), axis=row_axis).min(axis=col_axis)
2594 elif ord == -Inf:
2595 if row_axis > col_axis:
2596 row_axis -= 1
2597 ret = add.reduce(abs(x), axis=col_axis).min(axis=row_axis)
2598 elif ord in [None, 'fro', 'f']:
2599 ret = sqrt(add.reduce((x.conj() * x).real, axis=axis))
2600 elif ord == 'nuc':
2601 ret = _multi_svd_norm(x, row_axis, col_axis, sum)
2602 else:
2603 raise ValueError("Invalid norm order for matrices.")
2604 if keepdims:
2605 ret_shape = list(x.shape)
2606 ret_shape[axis[0]] = 1
2607 ret_shape[axis[1]] = 1
2608 ret = ret.reshape(ret_shape)
2609 return ret
2610 else:
2611 raise ValueError("Improper number of dimensions to norm.")
2614# multi_dot
2616def _multidot_dispatcher(arrays, *, out=None):
2617 yield from arrays
2618 yield out
2621@array_function_dispatch(_multidot_dispatcher)
2622def multi_dot(arrays, *, out=None):
2623 """
2624 Compute the dot product of two or more arrays in a single function call,
2625 while automatically selecting the fastest evaluation order.
2627 `multi_dot` chains `numpy.dot` and uses optimal parenthesization
2628 of the matrices [1]_ [2]_. Depending on the shapes of the matrices,
2629 this can speed up the multiplication a lot.
2631 If the first argument is 1-D it is treated as a row vector.
2632 If the last argument is 1-D it is treated as a column vector.
2633 The other arguments must be 2-D.
2635 Think of `multi_dot` as::
2637 def multi_dot(arrays): return functools.reduce(np.dot, arrays)
2640 Parameters
2641 ----------
2642 arrays : sequence of array_like
2643 If the first argument is 1-D it is treated as row vector.
2644 If the last argument is 1-D it is treated as column vector.
2645 The other arguments must be 2-D.
2646 out : ndarray, optional
2647 Output argument. This must have the exact kind that would be returned
2648 if it was not used. In particular, it must have the right type, must be
2649 C-contiguous, and its dtype must be the dtype that would be returned
2650 for `dot(a, b)`. This is a performance feature. Therefore, if these
2651 conditions are not met, an exception is raised, instead of attempting
2652 to be flexible.
2654 .. versionadded:: 1.19.0
2656 Returns
2657 -------
2658 output : ndarray
2659 Returns the dot product of the supplied arrays.
2661 See Also
2662 --------
2663 dot : dot multiplication with two arguments.
2665 References
2666 ----------
2668 .. [1] Cormen, "Introduction to Algorithms", Chapter 15.2, p. 370-378
2669 .. [2] https://en.wikipedia.org/wiki/Matrix_chain_multiplication
2671 Examples
2672 --------
2673 `multi_dot` allows you to write::
2675 >>> from numpy.linalg import multi_dot
2676 >>> # Prepare some data
2677 >>> A = np.random.random((10000, 100))
2678 >>> B = np.random.random((100, 1000))
2679 >>> C = np.random.random((1000, 5))
2680 >>> D = np.random.random((5, 333))
2681 >>> # the actual dot multiplication
2682 >>> _ = multi_dot([A, B, C, D])
2684 instead of::
2686 >>> _ = np.dot(np.dot(np.dot(A, B), C), D)
2687 >>> # or
2688 >>> _ = A.dot(B).dot(C).dot(D)
2690 Notes
2691 -----
2692 The cost for a matrix multiplication can be calculated with the
2693 following function::
2695 def cost(A, B):
2696 return A.shape[0] * A.shape[1] * B.shape[1]
2698 Assume we have three matrices
2699 :math:`A_{10x100}, B_{100x5}, C_{5x50}`.
2701 The costs for the two different parenthesizations are as follows::
2703 cost((AB)C) = 10*100*5 + 10*5*50 = 5000 + 2500 = 7500
2704 cost(A(BC)) = 10*100*50 + 100*5*50 = 50000 + 25000 = 75000
2706 """
2707 n = len(arrays)
2708 # optimization only makes sense for len(arrays) > 2
2709 if n < 2:
2710 raise ValueError("Expecting at least two arrays.")
2711 elif n == 2:
2712 return dot(arrays[0], arrays[1], out=out)
2714 arrays = [asanyarray(a) for a in arrays]
2716 # save original ndim to reshape the result array into the proper form later
2717 ndim_first, ndim_last = arrays[0].ndim, arrays[-1].ndim
2718 # Explicitly convert vectors to 2D arrays to keep the logic of the internal
2719 # _multi_dot_* functions as simple as possible.
2720 if arrays[0].ndim == 1:
2721 arrays[0] = atleast_2d(arrays[0])
2722 if arrays[-1].ndim == 1:
2723 arrays[-1] = atleast_2d(arrays[-1]).T
2724 _assert_2d(*arrays)
2726 # _multi_dot_three is much faster than _multi_dot_matrix_chain_order
2727 if n == 3:
2728 result = _multi_dot_three(arrays[0], arrays[1], arrays[2], out=out)
2729 else:
2730 order = _multi_dot_matrix_chain_order(arrays)
2731 result = _multi_dot(arrays, order, 0, n - 1, out=out)
2733 # return proper shape
2734 if ndim_first == 1 and ndim_last == 1:
2735 return result[0, 0] # scalar
2736 elif ndim_first == 1 or ndim_last == 1:
2737 return result.ravel() # 1-D
2738 else:
2739 return result
2742def _multi_dot_three(A, B, C, out=None):
2743 """
2744 Find the best order for three arrays and do the multiplication.
2746 For three arguments `_multi_dot_three` is approximately 15 times faster
2747 than `_multi_dot_matrix_chain_order`
2749 """
2750 a0, a1b0 = A.shape
2751 b1c0, c1 = C.shape
2752 # cost1 = cost((AB)C) = a0*a1b0*b1c0 + a0*b1c0*c1
2753 cost1 = a0 * b1c0 * (a1b0 + c1)
2754 # cost2 = cost(A(BC)) = a1b0*b1c0*c1 + a0*a1b0*c1
2755 cost2 = a1b0 * c1 * (a0 + b1c0)
2757 if cost1 < cost2:
2758 return dot(dot(A, B), C, out=out)
2759 else:
2760 return dot(A, dot(B, C), out=out)
2763def _multi_dot_matrix_chain_order(arrays, return_costs=False):
2764 """
2765 Return a np.array that encodes the optimal order of mutiplications.
2767 The optimal order array is then used by `_multi_dot()` to do the
2768 multiplication.
2770 Also return the cost matrix if `return_costs` is `True`
2772 The implementation CLOSELY follows Cormen, "Introduction to Algorithms",
2773 Chapter 15.2, p. 370-378. Note that Cormen uses 1-based indices.
2775 cost[i, j] = min([
2776 cost[prefix] + cost[suffix] + cost_mult(prefix, suffix)
2777 for k in range(i, j)])
2779 """
2780 n = len(arrays)
2781 # p stores the dimensions of the matrices
2782 # Example for p: A_{10x100}, B_{100x5}, C_{5x50} --> p = [10, 100, 5, 50]
2783 p = [a.shape[0] for a in arrays] + [arrays[-1].shape[1]]
2784 # m is a matrix of costs of the subproblems
2785 # m[i,j]: min number of scalar multiplications needed to compute A_{i..j}
2786 m = zeros((n, n), dtype=double)
2787 # s is the actual ordering
2788 # s[i, j] is the value of k at which we split the product A_i..A_j
2789 s = empty((n, n), dtype=intp)
2791 for l in range(1, n):
2792 for i in range(n - l):
2793 j = i + l
2794 m[i, j] = Inf
2795 for k in range(i, j):
2796 q = m[i, k] + m[k+1, j] + p[i]*p[k+1]*p[j+1]
2797 if q < m[i, j]:
2798 m[i, j] = q
2799 s[i, j] = k # Note that Cormen uses 1-based index
2801 return (s, m) if return_costs else s
2804def _multi_dot(arrays, order, i, j, out=None):
2805 """Actually do the multiplication with the given order."""
2806 if i == j:
2807 # the initial call with non-None out should never get here
2808 assert out is None
2810 return arrays[i]
2811 else:
2812 return dot(_multi_dot(arrays, order, i, order[i, j]),
2813 _multi_dot(arrays, order, order[i, j] + 1, j),
2814 out=out)