Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/numpy/linalg/_linalg.py: 20%
752 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-04-09 06:12 +0000
1"""Lite version of scipy.linalg.
3Notes
4-----
5This module is a lite version of the linalg.py module in SciPy which
6contains high-level Python interface to the LAPACK library. The lite
7version only accesses the following LAPACK functions: dgesv, zgesv,
8dgeev, zgeev, dgesdd, zgesdd, dgelsd, zgelsd, dsyevd, zheevd, dgetrf,
9zgetrf, dpotrf, zpotrf, dgeqrf, zgeqrf, zungqr, dorgqr.
10"""
12__all__ = ['matrix_power', 'solve', 'tensorsolve', 'tensorinv', 'inv',
13 'cholesky', 'eigvals', 'eigvalsh', 'pinv', 'slogdet', 'det',
14 'svd', 'svdvals', 'eig', 'eigh', 'lstsq', 'norm', 'qr', 'cond',
15 'matrix_rank', 'LinAlgError', 'multi_dot', 'trace', 'diagonal',
16 'cross', 'outer', 'tensordot', 'matmul', 'matrix_transpose',
17 'matrix_norm', 'vector_norm', 'vecdot']
19import functools
20import operator
21import warnings
22from typing import NamedTuple, Any
24from numpy._utils import set_module
25from numpy._core import (
26 array, asarray, zeros, empty, empty_like, intc, single, double,
27 csingle, cdouble, inexact, complexfloating, newaxis, all, inf, dot,
28 add, multiply, sqrt, sum, isfinite, finfo, errstate, moveaxis, amin,
29 amax, prod, abs, atleast_2d, intp, asanyarray, object_, matmul,
30 swapaxes, divide, count_nonzero, isnan, sign, argsort, sort,
31 reciprocal, overrides, diagonal as _core_diagonal, trace as _core_trace,
32 cross as _core_cross, outer as _core_outer, tensordot as _core_tensordot,
33 matmul as _core_matmul, matrix_transpose as _core_matrix_transpose,
34 transpose as _core_transpose, vecdot as _core_vecdot,
35)
36from numpy._globals import _NoValue
37from numpy.lib._twodim_base_impl import triu, eye
38from numpy.lib.array_utils import normalize_axis_index, normalize_axis_tuple
39from numpy.linalg import _umath_linalg
41from numpy._typing import NDArray
43class EigResult(NamedTuple):
44 eigenvalues: NDArray[Any]
45 eigenvectors: NDArray[Any]
47class EighResult(NamedTuple):
48 eigenvalues: NDArray[Any]
49 eigenvectors: NDArray[Any]
51class QRResult(NamedTuple):
52 Q: NDArray[Any]
53 R: NDArray[Any]
55class SlogdetResult(NamedTuple):
56 sign: NDArray[Any]
57 logabsdet: NDArray[Any]
59class SVDResult(NamedTuple):
60 U: NDArray[Any]
61 S: NDArray[Any]
62 Vh: NDArray[Any]
65array_function_dispatch = functools.partial(
66 overrides.array_function_dispatch, module='numpy.linalg'
67)
70fortran_int = intc
73@set_module('numpy.linalg')
74class LinAlgError(ValueError):
75 """
76 Generic Python-exception-derived object raised by linalg functions.
78 General purpose exception class, derived from Python's ValueError
79 class, programmatically raised in linalg functions when a Linear
80 Algebra-related condition would prevent further correct execution of the
81 function.
83 Parameters
84 ----------
85 None
87 Examples
88 --------
89 >>> from numpy import linalg as LA
90 >>> LA.inv(np.zeros((2,2)))
91 Traceback (most recent call last):
92 File "<stdin>", line 1, in <module>
93 File "...linalg.py", line 350,
94 in inv return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))
95 File "...linalg.py", line 249,
96 in solve
97 raise LinAlgError('Singular matrix')
98 numpy.linalg.LinAlgError: Singular matrix
100 """
103def _raise_linalgerror_singular(err, flag):
104 raise LinAlgError("Singular matrix")
106def _raise_linalgerror_nonposdef(err, flag):
107 raise LinAlgError("Matrix is not positive definite")
109def _raise_linalgerror_eigenvalues_nonconvergence(err, flag):
110 raise LinAlgError("Eigenvalues did not converge")
112def _raise_linalgerror_svd_nonconvergence(err, flag):
113 raise LinAlgError("SVD did not converge")
115def _raise_linalgerror_lstsq(err, flag):
116 raise LinAlgError("SVD did not converge in Linear Least Squares")
118def _raise_linalgerror_qr(err, flag):
119 raise LinAlgError("Incorrect argument found while performing "
120 "QR factorization")
123def _makearray(a):
124 new = asarray(a)
125 wrap = getattr(a, "__array_wrap__", new.__array_wrap__)
126 return new, wrap
128def isComplexType(t):
129 return issubclass(t, complexfloating)
132_real_types_map = {single: single,
133 double: double,
134 csingle: single,
135 cdouble: double}
137_complex_types_map = {single: csingle,
138 double: cdouble,
139 csingle: csingle,
140 cdouble: cdouble}
142def _realType(t, default=double):
143 return _real_types_map.get(t, default)
145def _complexType(t, default=cdouble):
146 return _complex_types_map.get(t, default)
148def _commonType(*arrays):
149 # in lite version, use higher precision (always double or cdouble)
150 result_type = single
151 is_complex = False
152 for a in arrays:
153 type_ = a.dtype.type
154 if issubclass(type_, inexact):
155 if isComplexType(type_):
156 is_complex = True
157 rt = _realType(type_, default=None)
158 if rt is double:
159 result_type = double
160 elif rt is None:
161 # unsupported inexact scalar
162 raise TypeError("array type %s is unsupported in linalg" %
163 (a.dtype.name,))
164 else:
165 result_type = double
166 if is_complex:
167 result_type = _complex_types_map[result_type]
168 return cdouble, result_type
169 else:
170 return double, result_type
173def _to_native_byte_order(*arrays):
174 ret = []
175 for arr in arrays:
176 if arr.dtype.byteorder not in ('=', '|'):
177 ret.append(asarray(arr, dtype=arr.dtype.newbyteorder('=')))
178 else:
179 ret.append(arr)
180 if len(ret) == 1:
181 return ret[0]
182 else:
183 return ret
186def _assert_2d(*arrays):
187 for a in arrays:
188 if a.ndim != 2:
189 raise LinAlgError('%d-dimensional array given. Array must be '
190 'two-dimensional' % a.ndim)
192def _assert_stacked_2d(*arrays):
193 for a in arrays:
194 if a.ndim < 2:
195 raise LinAlgError('%d-dimensional array given. Array must be '
196 'at least two-dimensional' % a.ndim)
198def _assert_stacked_square(*arrays):
199 for a in arrays:
200 m, n = a.shape[-2:]
201 if m != n:
202 raise LinAlgError('Last 2 dimensions of the array must be square')
204def _assert_finite(*arrays):
205 for a in arrays:
206 if not isfinite(a).all():
207 raise LinAlgError("Array must not contain infs or NaNs")
209def _is_empty_2d(arr):
210 # check size first for efficiency
211 return arr.size == 0 and prod(arr.shape[-2:]) == 0
214def transpose(a):
215 """
216 Transpose each matrix in a stack of matrices.
218 Unlike np.transpose, this only swaps the last two axes, rather than all of
219 them
221 Parameters
222 ----------
223 a : (...,M,N) array_like
225 Returns
226 -------
227 aT : (...,N,M) ndarray
228 """
229 return swapaxes(a, -1, -2)
231# Linear equations
233def _tensorsolve_dispatcher(a, b, axes=None):
234 return (a, b)
237@array_function_dispatch(_tensorsolve_dispatcher)
238def tensorsolve(a, b, axes=None):
239 """
240 Solve the tensor equation ``a x = b`` for x.
242 It is assumed that all indices of `x` are summed over in the product,
243 together with the rightmost indices of `a`, as is done in, for example,
244 ``tensordot(a, x, axes=x.ndim)``.
246 Parameters
247 ----------
248 a : array_like
249 Coefficient tensor, of shape ``b.shape + Q``. `Q`, a tuple, equals
250 the shape of that sub-tensor of `a` consisting of the appropriate
251 number of its rightmost indices, and must be such that
252 ``prod(Q) == prod(b.shape)`` (in which sense `a` is said to be
253 'square').
254 b : array_like
255 Right-hand tensor, which can be of any shape.
256 axes : tuple of ints, optional
257 Axes in `a` to reorder to the right, before inversion.
258 If None (default), no reordering is done.
260 Returns
261 -------
262 x : ndarray, shape Q
264 Raises
265 ------
266 LinAlgError
267 If `a` is singular or not 'square' (in the above sense).
269 See Also
270 --------
271 numpy.tensordot, tensorinv, numpy.einsum
273 Examples
274 --------
275 >>> a = np.eye(2*3*4)
276 >>> a.shape = (2*3, 4, 2, 3, 4)
277 >>> b = np.random.randn(2*3, 4)
278 >>> x = np.linalg.tensorsolve(a, b)
279 >>> x.shape
280 (2, 3, 4)
281 >>> np.allclose(np.tensordot(a, x, axes=3), b)
282 True
284 """
285 a, wrap = _makearray(a)
286 b = asarray(b)
287 an = a.ndim
289 if axes is not None:
290 allaxes = list(range(0, an))
291 for k in axes:
292 allaxes.remove(k)
293 allaxes.insert(an, k)
294 a = a.transpose(allaxes)
296 oldshape = a.shape[-(an-b.ndim):]
297 prod = 1
298 for k in oldshape:
299 prod *= k
301 if a.size != prod ** 2:
302 raise LinAlgError(
303 "Input arrays must satisfy the requirement \
304 prod(a.shape[b.ndim:]) == prod(a.shape[:b.ndim])"
305 )
307 a = a.reshape(prod, prod)
308 b = b.ravel()
309 res = wrap(solve(a, b))
310 res.shape = oldshape
311 return res
314def _solve_dispatcher(a, b):
315 return (a, b)
318@array_function_dispatch(_solve_dispatcher)
319def solve(a, b):
320 """
321 Solve a linear matrix equation, or system of linear scalar equations.
323 Computes the "exact" solution, `x`, of the well-determined, i.e., full
324 rank, linear matrix equation `ax = b`.
326 Parameters
327 ----------
328 a : (..., M, M) array_like
329 Coefficient matrix.
330 b : {(M,), (..., M, K)}, array_like
331 Ordinate or "dependent variable" values.
333 Returns
334 -------
335 x : {(..., M,), (..., M, K)} ndarray
336 Solution to the system a x = b. Returned shape is (..., M) if b is
337 shape (M,) and (..., M, K) if b is (..., M, K), where the "..." part is
338 broadcasted between a and 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 .. versionchanged:: 2.0
366 The b array is only treated as a shape (M,) column vector if it is
367 exactly 1-dimensional. In all other instances it is treated as a stack
368 of (M, K) matrices. Previously b would be treated as a stack of (M,)
369 vectors if b.ndim was equal to a.ndim - 1.
371 References
372 ----------
373 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
374 FL, Academic Press, Inc., 1980, pg. 22.
376 Examples
377 --------
378 Solve the system of equations:
379 ``x0 + 2 * x1 = 1`` and
380 ``3 * x0 + 5 * x1 = 2``:
382 >>> a = np.array([[1, 2], [3, 5]])
383 >>> b = np.array([1, 2])
384 >>> x = np.linalg.solve(a, b)
385 >>> x
386 array([-1., 1.])
388 Check that the solution is correct:
390 >>> np.allclose(np.dot(a, x), b)
391 True
393 """
394 a, _ = _makearray(a)
395 _assert_stacked_2d(a)
396 _assert_stacked_square(a)
397 b, wrap = _makearray(b)
398 t, result_t = _commonType(a, b)
400 # We use the b = (..., M,) logic, only if the number of extra dimensions
401 # match exactly
402 if b.ndim == 1:
403 gufunc = _umath_linalg.solve1
404 else:
405 gufunc = _umath_linalg.solve
407 signature = 'DD->D' if isComplexType(t) else 'dd->d'
408 with errstate(call=_raise_linalgerror_singular, invalid='call',
409 over='ignore', divide='ignore', under='ignore'):
410 r = gufunc(a, b, signature=signature)
412 return wrap(r.astype(result_t, copy=False))
415def _tensorinv_dispatcher(a, ind=None):
416 return (a,)
419@array_function_dispatch(_tensorinv_dispatcher)
420def tensorinv(a, ind=2):
421 """
422 Compute the 'inverse' of an N-dimensional array.
424 The result is an inverse for `a` relative to the tensordot operation
425 ``tensordot(a, b, ind)``, i. e., up to floating-point accuracy,
426 ``tensordot(tensorinv(a), a, ind)`` is the "identity" tensor for the
427 tensordot operation.
429 Parameters
430 ----------
431 a : array_like
432 Tensor to 'invert'. Its shape must be 'square', i. e.,
433 ``prod(a.shape[:ind]) == prod(a.shape[ind:])``.
434 ind : int, optional
435 Number of first indices that are involved in the inverse sum.
436 Must be a positive integer, default is 2.
438 Returns
439 -------
440 b : ndarray
441 `a`'s tensordot inverse, shape ``a.shape[ind:] + a.shape[:ind]``.
443 Raises
444 ------
445 LinAlgError
446 If `a` is singular or not 'square' (in the above sense).
448 See Also
449 --------
450 numpy.tensordot, tensorsolve
452 Examples
453 --------
454 >>> a = np.eye(4*6)
455 >>> a.shape = (4, 6, 8, 3)
456 >>> ainv = np.linalg.tensorinv(a, ind=2)
457 >>> ainv.shape
458 (8, 3, 4, 6)
459 >>> b = np.random.randn(4, 6)
460 >>> np.allclose(np.tensordot(ainv, b), np.linalg.tensorsolve(a, b))
461 True
463 >>> a = np.eye(4*6)
464 >>> a.shape = (24, 8, 3)
465 >>> ainv = np.linalg.tensorinv(a, ind=1)
466 >>> ainv.shape
467 (8, 3, 24)
468 >>> b = np.random.randn(24)
469 >>> np.allclose(np.tensordot(ainv, b, 1), np.linalg.tensorsolve(a, b))
470 True
472 """
473 a = asarray(a)
474 oldshape = a.shape
475 prod = 1
476 if ind > 0:
477 invshape = oldshape[ind:] + oldshape[:ind]
478 for k in oldshape[ind:]:
479 prod *= k
480 else:
481 raise ValueError("Invalid ind argument.")
482 a = a.reshape(prod, -1)
483 ia = inv(a)
484 return ia.reshape(*invshape)
487# Matrix inversion
489def _unary_dispatcher(a):
490 return (a,)
493@array_function_dispatch(_unary_dispatcher)
494def inv(a):
495 """
496 Compute the inverse of a matrix.
498 Given a square matrix `a`, return the matrix `ainv` satisfying
499 ``a @ ainv = ainv @ a = eye(a.shape[0])``.
501 Parameters
502 ----------
503 a : (..., M, M) array_like
504 Matrix to be inverted.
506 Returns
507 -------
508 ainv : (..., M, M) ndarray or matrix
509 Inverse of the matrix `a`.
511 Raises
512 ------
513 LinAlgError
514 If `a` is not square or inversion fails.
516 See Also
517 --------
518 scipy.linalg.inv : Similar function in SciPy.
519 numpy.linalg.cond : Compute the condition number of a matrix.
520 numpy.linalg.svd : Compute the singular value decomposition of a matrix.
522 Notes
523 -----
525 .. versionadded:: 1.8.0
527 Broadcasting rules apply, see the `numpy.linalg` documentation for
528 details.
530 If `a` is detected to be singular, a `LinAlgError` is raised. If `a` is
531 ill-conditioned, a `LinAlgError` may or may not be raised, and results may
532 be inaccurate due to floating-point errors.
534 References
535 ----------
536 .. [1] Wikipedia, "Condition number",
537 https://en.wikipedia.org/wiki/Condition_number
539 Examples
540 --------
541 >>> from numpy.linalg import inv
542 >>> a = np.array([[1., 2.], [3., 4.]])
543 >>> ainv = inv(a)
544 >>> np.allclose(a @ ainv, np.eye(2))
545 True
546 >>> np.allclose(ainv @ a, np.eye(2))
547 True
549 If a is a matrix object, then the return value is a matrix as well:
551 >>> ainv = inv(np.matrix(a))
552 >>> ainv
553 matrix([[-2. , 1. ],
554 [ 1.5, -0.5]])
556 Inverses of several matrices can be computed at once:
558 >>> a = np.array([[[1., 2.], [3., 4.]], [[1, 3], [3, 5]]])
559 >>> inv(a)
560 array([[[-2. , 1. ],
561 [ 1.5 , -0.5 ]],
562 [[-1.25, 0.75],
563 [ 0.75, -0.25]]])
565 If a matrix is close to singular, the computed inverse may not satisfy
566 ``a @ ainv = ainv @ a = eye(a.shape[0])`` even if a `LinAlgError`
567 is not raised:
569 >>> a = np.array([[2,4,6],[2,0,2],[6,8,14]])
570 >>> inv(a) # No errors raised
571 array([[-1.12589991e+15, -5.62949953e+14, 5.62949953e+14],
572 [-1.12589991e+15, -5.62949953e+14, 5.62949953e+14],
573 [ 1.12589991e+15, 5.62949953e+14, -5.62949953e+14]])
574 >>> a @ inv(a)
575 array([[ 0. , -0.5 , 0. ], # may vary
576 [-0.5 , 0.625, 0.25 ],
577 [ 0. , 0. , 1. ]])
579 To detect ill-conditioned matrices, you can use `numpy.linalg.cond` to
580 compute its *condition number* [1]_. The larger the condition number, the
581 more ill-conditioned the matrix is. As a rule of thumb, if the condition
582 number ``cond(a) = 10**k``, then you may lose up to ``k`` digits of
583 accuracy on top of what would be lost to the numerical method due to loss
584 of precision from arithmetic methods.
586 >>> from numpy.linalg import cond
587 >>> cond(a)
588 np.float64(8.659885634118668e+17) # may vary
590 It is also possible to detect ill-conditioning by inspecting the matrix's
591 singular values directly. The ratio between the largest and the smallest
592 singular value is the condition number:
594 >>> from numpy.linalg import svd
595 >>> sigma = svd(a, compute_uv=False) # Do not compute singular vectors
596 >>> sigma.max()/sigma.min()
597 8.659885634118668e+17 # may vary
599 """
600 a, wrap = _makearray(a)
601 _assert_stacked_2d(a)
602 _assert_stacked_square(a)
603 t, result_t = _commonType(a)
605 signature = 'D->D' if isComplexType(t) else 'd->d'
606 with errstate(call=_raise_linalgerror_singular, invalid='call',
607 over='ignore', divide='ignore', under='ignore'):
608 ainv = _umath_linalg.inv(a, signature=signature)
609 return wrap(ainv.astype(result_t, copy=False))
612def _matrix_power_dispatcher(a, n):
613 return (a,)
616@array_function_dispatch(_matrix_power_dispatcher)
617def matrix_power(a, n):
618 """
619 Raise a square matrix to the (integer) power `n`.
621 For positive integers `n`, the power is computed by repeated matrix
622 squarings and matrix multiplications. If ``n == 0``, the identity matrix
623 of the same shape as M is returned. If ``n < 0``, the inverse
624 is computed and then raised to the ``abs(n)``.
626 .. note:: Stacks of object matrices are not currently supported.
628 Parameters
629 ----------
630 a : (..., M, M) array_like
631 Matrix to be "powered".
632 n : int
633 The exponent can be any integer or long integer, positive,
634 negative, or zero.
636 Returns
637 -------
638 a**n : (..., M, M) ndarray or matrix object
639 The return value is the same shape and type as `M`;
640 if the exponent is positive or zero then the type of the
641 elements is the same as those of `M`. If the exponent is
642 negative the elements are floating-point.
644 Raises
645 ------
646 LinAlgError
647 For matrices that are not square or that (for negative powers) cannot
648 be inverted numerically.
650 Examples
651 --------
652 >>> from numpy.linalg import matrix_power
653 >>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit
654 >>> matrix_power(i, 3) # should = -i
655 array([[ 0, -1],
656 [ 1, 0]])
657 >>> matrix_power(i, 0)
658 array([[1, 0],
659 [0, 1]])
660 >>> matrix_power(i, -3) # should = 1/(-i) = i, but w/ f.p. elements
661 array([[ 0., 1.],
662 [-1., 0.]])
664 Somewhat more sophisticated example
666 >>> q = np.zeros((4, 4))
667 >>> q[0:2, 0:2] = -i
668 >>> q[2:4, 2:4] = i
669 >>> q # one of the three quaternion units not equal to 1
670 array([[ 0., -1., 0., 0.],
671 [ 1., 0., 0., 0.],
672 [ 0., 0., 0., 1.],
673 [ 0., 0., -1., 0.]])
674 >>> matrix_power(q, 2) # = -np.eye(4)
675 array([[-1., 0., 0., 0.],
676 [ 0., -1., 0., 0.],
677 [ 0., 0., -1., 0.],
678 [ 0., 0., 0., -1.]])
680 """
681 a = asanyarray(a)
682 _assert_stacked_2d(a)
683 _assert_stacked_square(a)
685 try:
686 n = operator.index(n)
687 except TypeError as e:
688 raise TypeError("exponent must be an integer") from e
690 # Fall back on dot for object arrays. Object arrays are not supported by
691 # the current implementation of matmul using einsum
692 if a.dtype != object:
693 fmatmul = matmul
694 elif a.ndim == 2:
695 fmatmul = dot
696 else:
697 raise NotImplementedError(
698 "matrix_power not supported for stacks of object arrays")
700 if n == 0:
701 a = empty_like(a)
702 a[...] = eye(a.shape[-2], dtype=a.dtype)
703 return a
705 elif n < 0:
706 a = inv(a)
707 n = abs(n)
709 # short-cuts.
710 if n == 1:
711 return a
713 elif n == 2:
714 return fmatmul(a, a)
716 elif n == 3:
717 return fmatmul(fmatmul(a, a), a)
719 # Use binary decomposition to reduce the number of matrix multiplications.
720 # Here, we iterate over the bits of n, from LSB to MSB, raise `a` to
721 # increasing powers of 2, and multiply into the result as needed.
722 z = result = None
723 while n > 0:
724 z = a if z is None else fmatmul(z, z)
725 n, bit = divmod(n, 2)
726 if bit:
727 result = z if result is None else fmatmul(result, z)
729 return result
732# Cholesky decomposition
734def _cholesky_dispatcher(a, /, *, upper=None):
735 return (a,)
738@array_function_dispatch(_cholesky_dispatcher)
739def cholesky(a, /, *, upper=False):
740 """
741 Cholesky decomposition.
743 Return the lower or upper Cholesky decomposition, ``L * L.H`` or
744 ``U.H * U``, of the square matrix ``a``, where ``L`` is lower-triangular,
745 ``U`` is upper-triangular, and ``.H`` is the conjugate transpose operator
746 (which is the ordinary transpose if ``a`` is real-valued). ``a`` must be
747 Hermitian (symmetric if real-valued) and positive-definite. No checking is
748 performed to verify whether ``a`` is Hermitian or not. In addition, only
749 the lower or upper-triangular and diagonal elements of ``a`` are used.
750 Only ``L`` or ``U`` is actually returned.
752 Parameters
753 ----------
754 a : (..., M, M) array_like
755 Hermitian (symmetric if all elements are real), positive-definite
756 input matrix.
757 upper : bool
758 If ``True``, the result must be the upper-triangular Cholesky factor.
759 If ``False``, the result must be the lower-triangular Cholesky factor.
760 Default: ``False``.
762 Returns
763 -------
764 L : (..., M, M) array_like
765 Lower or upper-triangular Cholesky factor of `a`. Returns a matrix
766 object if `a` is a matrix object.
768 Raises
769 ------
770 LinAlgError
771 If the decomposition fails, for example, if `a` is not
772 positive-definite.
774 See Also
775 --------
776 scipy.linalg.cholesky : Similar function in SciPy.
777 scipy.linalg.cholesky_banded : Cholesky decompose a banded Hermitian
778 positive-definite matrix.
779 scipy.linalg.cho_factor : Cholesky decomposition of a matrix, to use in
780 `scipy.linalg.cho_solve`.
782 Notes
783 -----
785 .. versionadded:: 1.8.0
787 Broadcasting rules apply, see the `numpy.linalg` documentation for
788 details.
790 The Cholesky decomposition is often used as a fast way of solving
792 .. math:: A \\mathbf{x} = \\mathbf{b}
794 (when `A` is both Hermitian/symmetric and positive-definite).
796 First, we solve for :math:`\\mathbf{y}` in
798 .. math:: L \\mathbf{y} = \\mathbf{b},
800 and then for :math:`\\mathbf{x}` in
802 .. math:: L^{H} \\mathbf{x} = \\mathbf{y}.
804 Examples
805 --------
806 >>> A = np.array([[1,-2j],[2j,5]])
807 >>> A
808 array([[ 1.+0.j, -0.-2.j],
809 [ 0.+2.j, 5.+0.j]])
810 >>> L = np.linalg.cholesky(A)
811 >>> L
812 array([[1.+0.j, 0.+0.j],
813 [0.+2.j, 1.+0.j]])
814 >>> np.dot(L, L.T.conj()) # verify that L * L.H = A
815 array([[1.+0.j, 0.-2.j],
816 [0.+2.j, 5.+0.j]])
817 >>> A = [[1,-2j],[2j,5]] # what happens if A is only array_like?
818 >>> np.linalg.cholesky(A) # an ndarray object is returned
819 array([[1.+0.j, 0.+0.j],
820 [0.+2.j, 1.+0.j]])
821 >>> # But a matrix object is returned if A is a matrix object
822 >>> np.linalg.cholesky(np.matrix(A))
823 matrix([[ 1.+0.j, 0.+0.j],
824 [ 0.+2.j, 1.+0.j]])
825 >>> # The upper-triangular Cholesky factor can also be obtained.
826 >>> np.linalg.cholesky(A, upper=True)
827 array([[1.-0.j, 0.-2.j],
828 [0.-0.j, 1.-0.j]])
830 """
831 gufunc = _umath_linalg.cholesky_up if upper else _umath_linalg.cholesky_lo
832 a, wrap = _makearray(a)
833 _assert_stacked_2d(a)
834 _assert_stacked_square(a)
835 t, result_t = _commonType(a)
836 signature = 'D->D' if isComplexType(t) else 'd->d'
837 with errstate(call=_raise_linalgerror_nonposdef, invalid='call',
838 over='ignore', divide='ignore', under='ignore'):
839 r = gufunc(a, signature=signature)
840 return wrap(r.astype(result_t, copy=False))
843# outer product
846def _outer_dispatcher(x1, x2):
847 return (x1, x2)
850@array_function_dispatch(_outer_dispatcher)
851def outer(x1, x2, /):
852 """
853 Compute the outer product of two vectors.
855 This function is Array API compatible. Compared to ``np.outer``
856 it accepts 1-dimensional inputs only.
858 Parameters
859 ----------
860 x1 : (M,) array_like
861 One-dimensional input array of size ``N``.
862 Must have a numeric data type.
863 x2 : (N,) array_like
864 One-dimensional input array of size ``M``.
865 Must have a numeric data type.
867 Returns
868 -------
869 out : (M, N) ndarray
870 ``out[i, j] = a[i] * b[j]``
872 See also
873 --------
874 outer
876 """
877 x1 = asanyarray(x1)
878 x2 = asanyarray(x2)
879 if x1.ndim != 1 or x2.ndim != 1:
880 raise ValueError(
881 "Input arrays must be one-dimensional, but they are "
882 f"{x1.ndim=} and {x2.ndim=}."
883 )
884 return _core_outer(x1, x2, out=None)
887# QR decomposition
890def _qr_dispatcher(a, mode=None):
891 return (a,)
894@array_function_dispatch(_qr_dispatcher)
895def qr(a, mode='reduced'):
896 """
897 Compute the qr factorization of a matrix.
899 Factor the matrix `a` as *qr*, where `q` is orthonormal and `r` is
900 upper-triangular.
902 Parameters
903 ----------
904 a : array_like, shape (..., M, N)
905 An array-like object with the dimensionality of at least 2.
906 mode : {'reduced', 'complete', 'r', 'raw'}, optional, default: 'reduced'
907 If K = min(M, N), then
909 * 'reduced' : returns Q, R with dimensions (..., M, K), (..., K, N)
910 * 'complete' : returns Q, R with dimensions (..., M, M), (..., M, N)
911 * 'r' : returns R only with dimensions (..., K, N)
912 * 'raw' : returns h, tau with dimensions (..., N, M), (..., K,)
914 The options 'reduced', 'complete, and 'raw' are new in numpy 1.8,
915 see the notes for more information. The default is 'reduced', and to
916 maintain backward compatibility with earlier versions of numpy both
917 it and the old default 'full' can be omitted. Note that array h
918 returned in 'raw' mode is transposed for calling Fortran. The
919 'economic' mode is deprecated. The modes 'full' and 'economic' may
920 be passed using only the first letter for backwards compatibility,
921 but all others must be spelled out. See the Notes for more
922 explanation.
925 Returns
926 -------
927 When mode is 'reduced' or 'complete', the result will be a namedtuple with
928 the attributes `Q` and `R`.
930 Q : ndarray of float or complex, optional
931 A matrix with orthonormal columns. When mode = 'complete' the
932 result is an orthogonal/unitary matrix depending on whether or not
933 a is real/complex. The determinant may be either +/- 1 in that
934 case. In case the number of dimensions in the input array is
935 greater than 2 then a stack of the matrices with above properties
936 is returned.
937 R : ndarray of float or complex, optional
938 The upper-triangular matrix or a stack of upper-triangular
939 matrices if the number of dimensions in the input array is greater
940 than 2.
941 (h, tau) : ndarrays of np.double or np.cdouble, optional
942 The array h contains the Householder reflectors that generate q
943 along with r. The tau array contains scaling factors for the
944 reflectors. In the deprecated 'economic' mode only h is returned.
946 Raises
947 ------
948 LinAlgError
949 If factoring fails.
951 See Also
952 --------
953 scipy.linalg.qr : Similar function in SciPy.
954 scipy.linalg.rq : Compute RQ decomposition of a matrix.
956 Notes
957 -----
958 This is an interface to the LAPACK routines ``dgeqrf``, ``zgeqrf``,
959 ``dorgqr``, and ``zungqr``.
961 For more information on the qr factorization, see for example:
962 https://en.wikipedia.org/wiki/QR_factorization
964 Subclasses of `ndarray` are preserved except for the 'raw' mode. So if
965 `a` is of type `matrix`, all the return values will be matrices too.
967 New 'reduced', 'complete', and 'raw' options for mode were added in
968 NumPy 1.8.0 and the old option 'full' was made an alias of 'reduced'. In
969 addition the options 'full' and 'economic' were deprecated. Because
970 'full' was the previous default and 'reduced' is the new default,
971 backward compatibility can be maintained by letting `mode` default.
972 The 'raw' option was added so that LAPACK routines that can multiply
973 arrays by q using the Householder reflectors can be used. Note that in
974 this case the returned arrays are of type np.double or np.cdouble and
975 the h array is transposed to be FORTRAN compatible. No routines using
976 the 'raw' return are currently exposed by numpy, but some are available
977 in lapack_lite and just await the necessary work.
979 Examples
980 --------
981 >>> a = np.random.randn(9, 6)
982 >>> Q, R = np.linalg.qr(a)
983 >>> np.allclose(a, np.dot(Q, R)) # a does equal QR
984 True
985 >>> R2 = np.linalg.qr(a, mode='r')
986 >>> np.allclose(R, R2) # mode='r' returns the same R as mode='full'
987 True
988 >>> a = np.random.normal(size=(3, 2, 2)) # Stack of 2 x 2 matrices as input
989 >>> Q, R = np.linalg.qr(a)
990 >>> Q.shape
991 (3, 2, 2)
992 >>> R.shape
993 (3, 2, 2)
994 >>> np.allclose(a, np.matmul(Q, R))
995 True
997 Example illustrating a common use of `qr`: solving of least squares
998 problems
1000 What are the least-squares-best `m` and `y0` in ``y = y0 + mx`` for
1001 the following data: {(0,1), (1,0), (1,2), (2,1)}. (Graph the points
1002 and you'll see that it should be y0 = 0, m = 1.) The answer is provided
1003 by solving the over-determined matrix equation ``Ax = b``, where::
1005 A = array([[0, 1], [1, 1], [1, 1], [2, 1]])
1006 x = array([[y0], [m]])
1007 b = array([[1], [0], [2], [1]])
1009 If A = QR such that Q is orthonormal (which is always possible via
1010 Gram-Schmidt), then ``x = inv(R) * (Q.T) * b``. (In numpy practice,
1011 however, we simply use `lstsq`.)
1013 >>> A = np.array([[0, 1], [1, 1], [1, 1], [2, 1]])
1014 >>> A
1015 array([[0, 1],
1016 [1, 1],
1017 [1, 1],
1018 [2, 1]])
1019 >>> b = np.array([1, 2, 2, 3])
1020 >>> Q, R = np.linalg.qr(A)
1021 >>> p = np.dot(Q.T, b)
1022 >>> np.dot(np.linalg.inv(R), p)
1023 array([ 1., 1.])
1025 """
1026 if mode not in ('reduced', 'complete', 'r', 'raw'):
1027 if mode in ('f', 'full'):
1028 # 2013-04-01, 1.8
1029 msg = "".join((
1030 "The 'full' option is deprecated in favor of 'reduced'.\n",
1031 "For backward compatibility let mode default."))
1032 warnings.warn(msg, DeprecationWarning, stacklevel=2)
1033 mode = 'reduced'
1034 elif mode in ('e', 'economic'):
1035 # 2013-04-01, 1.8
1036 msg = "The 'economic' option is deprecated."
1037 warnings.warn(msg, DeprecationWarning, stacklevel=2)
1038 mode = 'economic'
1039 else:
1040 raise ValueError(f"Unrecognized mode '{mode}'")
1042 a, wrap = _makearray(a)
1043 _assert_stacked_2d(a)
1044 m, n = a.shape[-2:]
1045 t, result_t = _commonType(a)
1046 a = a.astype(t, copy=True)
1047 a = _to_native_byte_order(a)
1048 mn = min(m, n)
1050 if m <= n:
1051 gufunc = _umath_linalg.qr_r_raw_m
1052 else:
1053 gufunc = _umath_linalg.qr_r_raw_n
1055 signature = 'D->D' if isComplexType(t) else 'd->d'
1056 with errstate(call=_raise_linalgerror_qr, invalid='call',
1057 over='ignore', divide='ignore', under='ignore'):
1058 tau = gufunc(a, signature=signature)
1060 # handle modes that don't return q
1061 if mode == 'r':
1062 r = triu(a[..., :mn, :])
1063 r = r.astype(result_t, copy=False)
1064 return wrap(r)
1066 if mode == 'raw':
1067 q = transpose(a)
1068 q = q.astype(result_t, copy=False)
1069 tau = tau.astype(result_t, copy=False)
1070 return wrap(q), tau
1072 if mode == 'economic':
1073 a = a.astype(result_t, copy=False)
1074 return wrap(a)
1076 # mc is the number of columns in the resulting q
1077 # matrix. If the mode is complete then it is
1078 # same as number of rows, and if the mode is reduced,
1079 # then it is the minimum of number of rows and columns.
1080 if mode == 'complete' and m > n:
1081 mc = m
1082 gufunc = _umath_linalg.qr_complete
1083 else:
1084 mc = mn
1085 gufunc = _umath_linalg.qr_reduced
1087 signature = 'DD->D' if isComplexType(t) else 'dd->d'
1088 with errstate(call=_raise_linalgerror_qr, invalid='call',
1089 over='ignore', divide='ignore', under='ignore'):
1090 q = gufunc(a, tau, signature=signature)
1091 r = triu(a[..., :mc, :])
1093 q = q.astype(result_t, copy=False)
1094 r = r.astype(result_t, copy=False)
1096 return QRResult(wrap(q), wrap(r))
1098# Eigenvalues
1101@array_function_dispatch(_unary_dispatcher)
1102def eigvals(a):
1103 """
1104 Compute the eigenvalues of a general matrix.
1106 Main difference between `eigvals` and `eig`: the eigenvectors aren't
1107 returned.
1109 Parameters
1110 ----------
1111 a : (..., M, M) array_like
1112 A complex- or real-valued matrix whose eigenvalues will be computed.
1114 Returns
1115 -------
1116 w : (..., M,) ndarray
1117 The eigenvalues, each repeated according to its multiplicity.
1118 They are not necessarily ordered, nor are they necessarily
1119 real for real matrices.
1121 Raises
1122 ------
1123 LinAlgError
1124 If the eigenvalue computation does not converge.
1126 See Also
1127 --------
1128 eig : eigenvalues and right eigenvectors of general arrays
1129 eigvalsh : eigenvalues of real symmetric or complex Hermitian
1130 (conjugate symmetric) arrays.
1131 eigh : eigenvalues and eigenvectors of real symmetric or complex
1132 Hermitian (conjugate symmetric) arrays.
1133 scipy.linalg.eigvals : Similar function in SciPy.
1135 Notes
1136 -----
1138 .. versionadded:: 1.8.0
1140 Broadcasting rules apply, see the `numpy.linalg` documentation for
1141 details.
1143 This is implemented using the ``_geev`` LAPACK routines which compute
1144 the eigenvalues and eigenvectors of general square arrays.
1146 Examples
1147 --------
1148 Illustration, using the fact that the eigenvalues of a diagonal matrix
1149 are its diagonal elements, that multiplying a matrix on the left
1150 by an orthogonal matrix, `Q`, and on the right by `Q.T` (the transpose
1151 of `Q`), preserves the eigenvalues of the "middle" matrix. In other words,
1152 if `Q` is orthogonal, then ``Q * A * Q.T`` has the same eigenvalues as
1153 ``A``:
1155 >>> from numpy import linalg as LA
1156 >>> x = np.random.random()
1157 >>> Q = np.array([[np.cos(x), -np.sin(x)], [np.sin(x), np.cos(x)]])
1158 >>> LA.norm(Q[0, :]), LA.norm(Q[1, :]), np.dot(Q[0, :],Q[1, :])
1159 (1.0, 1.0, 0.0)
1161 Now multiply a diagonal matrix by ``Q`` on one side and
1162 by ``Q.T`` on the other:
1164 >>> D = np.diag((-1,1))
1165 >>> LA.eigvals(D)
1166 array([-1., 1.])
1167 >>> A = np.dot(Q, D)
1168 >>> A = np.dot(A, Q.T)
1169 >>> LA.eigvals(A)
1170 array([ 1., -1.]) # random
1172 """
1173 a, wrap = _makearray(a)
1174 _assert_stacked_2d(a)
1175 _assert_stacked_square(a)
1176 _assert_finite(a)
1177 t, result_t = _commonType(a)
1179 signature = 'D->D' if isComplexType(t) else 'd->D'
1180 with errstate(call=_raise_linalgerror_eigenvalues_nonconvergence,
1181 invalid='call', over='ignore', divide='ignore',
1182 under='ignore'):
1183 w = _umath_linalg.eigvals(a, signature=signature)
1185 if not isComplexType(t):
1186 if all(w.imag == 0):
1187 w = w.real
1188 result_t = _realType(result_t)
1189 else:
1190 result_t = _complexType(result_t)
1192 return w.astype(result_t, copy=False)
1195def _eigvalsh_dispatcher(a, UPLO=None):
1196 return (a,)
1199@array_function_dispatch(_eigvalsh_dispatcher)
1200def eigvalsh(a, UPLO='L'):
1201 """
1202 Compute the eigenvalues of a complex Hermitian or real symmetric matrix.
1204 Main difference from eigh: the eigenvectors are not computed.
1206 Parameters
1207 ----------
1208 a : (..., M, M) array_like
1209 A complex- or real-valued matrix whose eigenvalues are to be
1210 computed.
1211 UPLO : {'L', 'U'}, optional
1212 Specifies whether the calculation is done with the lower triangular
1213 part of `a` ('L', default) or the upper triangular part ('U').
1214 Irrespective of this value only the real parts of the diagonal will
1215 be considered in the computation to preserve the notion of a Hermitian
1216 matrix. It therefore follows that the imaginary part of the diagonal
1217 will always be treated as zero.
1219 Returns
1220 -------
1221 w : (..., M,) ndarray
1222 The eigenvalues in ascending order, each repeated according to
1223 its multiplicity.
1225 Raises
1226 ------
1227 LinAlgError
1228 If the eigenvalue computation does not converge.
1230 See Also
1231 --------
1232 eigh : eigenvalues and eigenvectors of real symmetric or complex Hermitian
1233 (conjugate symmetric) arrays.
1234 eigvals : eigenvalues of general real or complex arrays.
1235 eig : eigenvalues and right eigenvectors of general real or complex
1236 arrays.
1237 scipy.linalg.eigvalsh : Similar function in SciPy.
1239 Notes
1240 -----
1242 .. versionadded:: 1.8.0
1244 Broadcasting rules apply, see the `numpy.linalg` documentation for
1245 details.
1247 The eigenvalues are computed using LAPACK routines ``_syevd``, ``_heevd``.
1249 Examples
1250 --------
1251 >>> from numpy import linalg as LA
1252 >>> a = np.array([[1, -2j], [2j, 5]])
1253 >>> LA.eigvalsh(a)
1254 array([ 0.17157288, 5.82842712]) # may vary
1256 >>> # demonstrate the treatment of the imaginary part of the diagonal
1257 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1258 >>> a
1259 array([[5.+2.j, 9.-2.j],
1260 [0.+2.j, 2.-1.j]])
1261 >>> # with UPLO='L' this is numerically equivalent to using LA.eigvals()
1262 >>> # with:
1263 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1264 >>> b
1265 array([[5.+0.j, 0.-2.j],
1266 [0.+2.j, 2.+0.j]])
1267 >>> wa = LA.eigvalsh(a)
1268 >>> wb = LA.eigvals(b)
1269 >>> wa; wb
1270 array([1., 6.])
1271 array([6.+0.j, 1.+0.j])
1273 """
1274 UPLO = UPLO.upper()
1275 if UPLO not in ('L', 'U'):
1276 raise ValueError("UPLO argument must be 'L' or 'U'")
1278 if UPLO == 'L':
1279 gufunc = _umath_linalg.eigvalsh_lo
1280 else:
1281 gufunc = _umath_linalg.eigvalsh_up
1283 a, wrap = _makearray(a)
1284 _assert_stacked_2d(a)
1285 _assert_stacked_square(a)
1286 t, result_t = _commonType(a)
1287 signature = 'D->d' if isComplexType(t) else 'd->d'
1288 with errstate(call=_raise_linalgerror_eigenvalues_nonconvergence,
1289 invalid='call', over='ignore', divide='ignore',
1290 under='ignore'):
1291 w = gufunc(a, signature=signature)
1292 return w.astype(_realType(result_t), copy=False)
1294def _convertarray(a):
1295 t, result_t = _commonType(a)
1296 a = a.astype(t).T.copy()
1297 return a, t, result_t
1300# Eigenvectors
1303@array_function_dispatch(_unary_dispatcher)
1304def eig(a):
1305 """
1306 Compute the eigenvalues and right eigenvectors of a square array.
1308 Parameters
1309 ----------
1310 a : (..., M, M) array
1311 Matrices for which the eigenvalues and right eigenvectors will
1312 be computed
1314 Returns
1315 -------
1316 A namedtuple with the following attributes:
1318 eigenvalues : (..., M) array
1319 The eigenvalues, each repeated according to its multiplicity.
1320 The eigenvalues are not necessarily ordered. The resulting
1321 array will be of complex type, unless the imaginary part is
1322 zero in which case it will be cast to a real type. When `a`
1323 is real the resulting eigenvalues will be real (0 imaginary
1324 part) or occur in conjugate pairs
1326 eigenvectors : (..., M, M) array
1327 The normalized (unit "length") eigenvectors, such that the
1328 column ``eigenvectors[:,i]`` is the eigenvector corresponding to the
1329 eigenvalue ``eigenvalues[i]``.
1331 Raises
1332 ------
1333 LinAlgError
1334 If the eigenvalue computation does not converge.
1336 See Also
1337 --------
1338 eigvals : eigenvalues of a non-symmetric array.
1339 eigh : eigenvalues and eigenvectors of a real symmetric or complex
1340 Hermitian (conjugate symmetric) array.
1341 eigvalsh : eigenvalues of a real symmetric or complex Hermitian
1342 (conjugate symmetric) array.
1343 scipy.linalg.eig : Similar function in SciPy that also solves the
1344 generalized eigenvalue problem.
1345 scipy.linalg.schur : Best choice for unitary and other non-Hermitian
1346 normal matrices.
1348 Notes
1349 -----
1351 .. versionadded:: 1.8.0
1353 Broadcasting rules apply, see the `numpy.linalg` documentation for
1354 details.
1356 This is implemented using the ``_geev`` LAPACK routines which compute
1357 the eigenvalues and eigenvectors of general square arrays.
1359 The number `w` is an eigenvalue of `a` if there exists a vector `v` such
1360 that ``a @ v = w * v``. Thus, the arrays `a`, `eigenvalues`, and
1361 `eigenvectors` satisfy the equations ``a @ eigenvectors[:,i] =
1362 eigenvalues[i] * eigenvectors[:,i]`` for :math:`i \\in \\{0,...,M-1\\}`.
1364 The array `eigenvectors` may not be of maximum rank, that is, some of the
1365 columns may be linearly dependent, although round-off error may obscure
1366 that fact. If the eigenvalues are all different, then theoretically the
1367 eigenvectors are linearly independent and `a` can be diagonalized by a
1368 similarity transformation using `eigenvectors`, i.e, ``inv(eigenvectors) @
1369 a @ eigenvectors`` is diagonal.
1371 For non-Hermitian normal matrices the SciPy function `scipy.linalg.schur`
1372 is preferred because the matrix `eigenvectors` is guaranteed to be
1373 unitary, which is not the case when using `eig`. The Schur factorization
1374 produces an upper triangular matrix rather than a diagonal matrix, but for
1375 normal matrices only the diagonal of the upper triangular matrix is
1376 needed, the rest is roundoff error.
1378 Finally, it is emphasized that `eigenvectors` consists of the *right* (as
1379 in right-hand side) eigenvectors of `a`. A vector `y` satisfying ``y.T @ a
1380 = z * y.T`` for some number `z` is called a *left* eigenvector of `a`,
1381 and, in general, the left and right eigenvectors of a matrix are not
1382 necessarily the (perhaps conjugate) transposes of each other.
1384 References
1385 ----------
1386 G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, FL,
1387 Academic Press, Inc., 1980, Various pp.
1389 Examples
1390 --------
1391 >>> from numpy import linalg as LA
1393 (Almost) trivial example with real eigenvalues and eigenvectors.
1395 >>> eigenvalues, eigenvectors = LA.eig(np.diag((1, 2, 3)))
1396 >>> eigenvalues
1397 array([1., 2., 3.])
1398 >>> eigenvectors
1399 array([[1., 0., 0.],
1400 [0., 1., 0.],
1401 [0., 0., 1.]])
1403 Real matrix possessing complex eigenvalues and eigenvectors;
1404 note that the eigenvalues are complex conjugates of each other.
1406 >>> eigenvalues, eigenvectors = LA.eig(np.array([[1, -1], [1, 1]]))
1407 >>> eigenvalues
1408 array([1.+1.j, 1.-1.j])
1409 >>> eigenvectors
1410 array([[0.70710678+0.j , 0.70710678-0.j ],
1411 [0. -0.70710678j, 0. +0.70710678j]])
1413 Complex-valued matrix with real eigenvalues (but complex-valued
1414 eigenvectors); note that ``a.conj().T == a``, i.e., `a` is Hermitian.
1416 >>> a = np.array([[1, 1j], [-1j, 1]])
1417 >>> eigenvalues, eigenvectors = LA.eig(a)
1418 >>> eigenvalues
1419 array([2.+0.j, 0.+0.j])
1420 >>> eigenvectors
1421 array([[ 0. +0.70710678j, 0.70710678+0.j ], # may vary
1422 [ 0.70710678+0.j , -0. +0.70710678j]])
1424 Be careful about round-off error!
1426 >>> a = np.array([[1 + 1e-9, 0], [0, 1 - 1e-9]])
1427 >>> # Theor. eigenvalues are 1 +/- 1e-9
1428 >>> eigenvalues, eigenvectors = LA.eig(a)
1429 >>> eigenvalues
1430 array([1., 1.])
1431 >>> eigenvectors
1432 array([[1., 0.],
1433 [0., 1.]])
1435 """
1436 a, wrap = _makearray(a)
1437 _assert_stacked_2d(a)
1438 _assert_stacked_square(a)
1439 _assert_finite(a)
1440 t, result_t = _commonType(a)
1442 signature = 'D->DD' if isComplexType(t) else 'd->DD'
1443 with errstate(call=_raise_linalgerror_eigenvalues_nonconvergence,
1444 invalid='call', over='ignore', divide='ignore',
1445 under='ignore'):
1446 w, vt = _umath_linalg.eig(a, signature=signature)
1448 if not isComplexType(t) and all(w.imag == 0.0):
1449 w = w.real
1450 vt = vt.real
1451 result_t = _realType(result_t)
1452 else:
1453 result_t = _complexType(result_t)
1455 vt = vt.astype(result_t, copy=False)
1456 return EigResult(w.astype(result_t, copy=False), wrap(vt))
1459@array_function_dispatch(_eigvalsh_dispatcher)
1460def eigh(a, UPLO='L'):
1461 """
1462 Return the eigenvalues and eigenvectors of a complex Hermitian
1463 (conjugate symmetric) or a real symmetric matrix.
1465 Returns two objects, a 1-D array containing the eigenvalues of `a`, and
1466 a 2-D square array or matrix (depending on the input type) of the
1467 corresponding eigenvectors (in columns).
1469 Parameters
1470 ----------
1471 a : (..., M, M) array
1472 Hermitian or real symmetric matrices whose eigenvalues and
1473 eigenvectors are to be computed.
1474 UPLO : {'L', 'U'}, optional
1475 Specifies whether the calculation is done with the lower triangular
1476 part of `a` ('L', default) or the upper triangular part ('U').
1477 Irrespective of this value only the real parts of the diagonal will
1478 be considered in the computation to preserve the notion of a Hermitian
1479 matrix. It therefore follows that the imaginary part of the diagonal
1480 will always be treated as zero.
1482 Returns
1483 -------
1484 A namedtuple with the following attributes:
1486 eigenvalues : (..., M) ndarray
1487 The eigenvalues in ascending order, each repeated according to
1488 its multiplicity.
1489 eigenvectors : {(..., M, M) ndarray, (..., M, M) matrix}
1490 The column ``eigenvectors[:, i]`` is the normalized eigenvector
1491 corresponding to the eigenvalue ``eigenvalues[i]``. Will return a
1492 matrix object if `a` is a matrix object.
1494 Raises
1495 ------
1496 LinAlgError
1497 If the eigenvalue computation does not converge.
1499 See Also
1500 --------
1501 eigvalsh : eigenvalues of real symmetric or complex Hermitian
1502 (conjugate symmetric) arrays.
1503 eig : eigenvalues and right eigenvectors for non-symmetric arrays.
1504 eigvals : eigenvalues of non-symmetric arrays.
1505 scipy.linalg.eigh : Similar function in SciPy (but also solves the
1506 generalized eigenvalue problem).
1508 Notes
1509 -----
1511 .. versionadded:: 1.8.0
1513 Broadcasting rules apply, see the `numpy.linalg` documentation for
1514 details.
1516 The eigenvalues/eigenvectors are computed using LAPACK routines ``_syevd``,
1517 ``_heevd``.
1519 The eigenvalues of real symmetric or complex Hermitian matrices are always
1520 real. [1]_ The array `eigenvalues` of (column) eigenvectors is unitary and
1521 `a`, `eigenvalues`, and `eigenvectors` satisfy the equations ``dot(a,
1522 eigenvectors[:, i]) = eigenvalues[i] * eigenvectors[:, i]``.
1524 References
1525 ----------
1526 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
1527 FL, Academic Press, Inc., 1980, pg. 222.
1529 Examples
1530 --------
1531 >>> from numpy import linalg as LA
1532 >>> a = np.array([[1, -2j], [2j, 5]])
1533 >>> a
1534 array([[ 1.+0.j, -0.-2.j],
1535 [ 0.+2.j, 5.+0.j]])
1536 >>> eigenvalues, eigenvectors = LA.eigh(a)
1537 >>> eigenvalues
1538 array([0.17157288, 5.82842712])
1539 >>> eigenvectors
1540 array([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1541 [ 0. +0.38268343j, 0. -0.92387953j]])
1543 >>> (np.dot(a, eigenvectors[:, 0]) -
1544 ... eigenvalues[0] * eigenvectors[:, 0]) # verify 1st eigenval/vec pair
1545 array([5.55111512e-17+0.0000000e+00j, 0.00000000e+00+1.2490009e-16j])
1546 >>> (np.dot(a, eigenvectors[:, 1]) -
1547 ... eigenvalues[1] * eigenvectors[:, 1]) # verify 2nd eigenval/vec pair
1548 array([0.+0.j, 0.+0.j])
1550 >>> A = np.matrix(a) # what happens if input is a matrix object
1551 >>> A
1552 matrix([[ 1.+0.j, -0.-2.j],
1553 [ 0.+2.j, 5.+0.j]])
1554 >>> eigenvalues, eigenvectors = LA.eigh(A)
1555 >>> eigenvalues
1556 array([0.17157288, 5.82842712])
1557 >>> eigenvectors
1558 matrix([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1559 [ 0. +0.38268343j, 0. -0.92387953j]])
1561 >>> # demonstrate the treatment of the imaginary part of the diagonal
1562 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1563 >>> a
1564 array([[5.+2.j, 9.-2.j],
1565 [0.+2.j, 2.-1.j]])
1566 >>> # with UPLO='L' this is numerically equivalent to using LA.eig() with:
1567 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1568 >>> b
1569 array([[5.+0.j, 0.-2.j],
1570 [0.+2.j, 2.+0.j]])
1571 >>> wa, va = LA.eigh(a)
1572 >>> wb, vb = LA.eig(b)
1573 >>> wa; wb
1574 array([1., 6.])
1575 array([6.+0.j, 1.+0.j])
1576 >>> va; vb
1577 array([[-0.4472136 +0.j , -0.89442719+0.j ], # may vary
1578 [ 0. +0.89442719j, 0. -0.4472136j ]])
1579 array([[ 0.89442719+0.j , -0. +0.4472136j],
1580 [-0. +0.4472136j, 0.89442719+0.j ]])
1582 """
1583 UPLO = UPLO.upper()
1584 if UPLO not in ('L', 'U'):
1585 raise ValueError("UPLO argument must be 'L' or 'U'")
1587 a, wrap = _makearray(a)
1588 _assert_stacked_2d(a)
1589 _assert_stacked_square(a)
1590 t, result_t = _commonType(a)
1592 if UPLO == 'L':
1593 gufunc = _umath_linalg.eigh_lo
1594 else:
1595 gufunc = _umath_linalg.eigh_up
1597 signature = 'D->dD' if isComplexType(t) else 'd->dd'
1598 with errstate(call=_raise_linalgerror_eigenvalues_nonconvergence,
1599 invalid='call', over='ignore', divide='ignore',
1600 under='ignore'):
1601 w, vt = gufunc(a, signature=signature)
1602 w = w.astype(_realType(result_t), copy=False)
1603 vt = vt.astype(result_t, copy=False)
1604 return EighResult(w, wrap(vt))
1607# Singular value decomposition
1609def _svd_dispatcher(a, full_matrices=None, compute_uv=None, hermitian=None):
1610 return (a,)
1613@array_function_dispatch(_svd_dispatcher)
1614def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
1615 """
1616 Singular Value Decomposition.
1618 When `a` is a 2D array, and ``full_matrices=False``, then it is
1619 factorized as ``u @ np.diag(s) @ vh = (u * s) @ vh``, where
1620 `u` and the Hermitian transpose of `vh` are 2D arrays with
1621 orthonormal columns and `s` is a 1D array of `a`'s singular
1622 values. When `a` is higher-dimensional, SVD is applied in
1623 stacked mode as explained below.
1625 Parameters
1626 ----------
1627 a : (..., M, N) array_like
1628 A real or complex array with ``a.ndim >= 2``.
1629 full_matrices : bool, optional
1630 If True (default), `u` and `vh` have the shapes ``(..., M, M)`` and
1631 ``(..., N, N)``, respectively. Otherwise, the shapes are
1632 ``(..., M, K)`` and ``(..., K, N)``, respectively, where
1633 ``K = min(M, N)``.
1634 compute_uv : bool, optional
1635 Whether or not to compute `u` and `vh` in addition to `s`. True
1636 by default.
1637 hermitian : bool, optional
1638 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
1639 enabling a more efficient method for finding singular values.
1640 Defaults to False.
1642 .. versionadded:: 1.17.0
1644 Returns
1645 -------
1646 When `compute_uv` is True, the result is a namedtuple with the following
1647 attribute names:
1649 U : { (..., M, M), (..., M, K) } array
1650 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1651 size as those of the input `a`. The size of the last two dimensions
1652 depends on the value of `full_matrices`. Only returned when
1653 `compute_uv` is True.
1654 S : (..., K) array
1655 Vector(s) with the singular values, within each vector sorted in
1656 descending order. The first ``a.ndim - 2`` dimensions have the same
1657 size as those of the input `a`.
1658 Vh : { (..., N, N), (..., K, N) } array
1659 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1660 size as those of the input `a`. The size of the last two dimensions
1661 depends on the value of `full_matrices`. Only returned when
1662 `compute_uv` is True.
1664 Raises
1665 ------
1666 LinAlgError
1667 If SVD computation does not converge.
1669 See Also
1670 --------
1671 scipy.linalg.svd : Similar function in SciPy.
1672 scipy.linalg.svdvals : Compute singular values of a matrix.
1674 Notes
1675 -----
1677 .. versionchanged:: 1.8.0
1678 Broadcasting rules apply, see the `numpy.linalg` documentation for
1679 details.
1681 The decomposition is performed using LAPACK routine ``_gesdd``.
1683 SVD is usually described for the factorization of a 2D matrix :math:`A`.
1684 The higher-dimensional case will be discussed below. In the 2D case, SVD is
1685 written as :math:`A = U S V^H`, where :math:`A = a`, :math:`U= u`,
1686 :math:`S= \\mathtt{np.diag}(s)` and :math:`V^H = vh`. The 1D array `s`
1687 contains the singular values of `a` and `u` and `vh` are unitary. The rows
1688 of `vh` are the eigenvectors of :math:`A^H A` and the columns of `u` are
1689 the eigenvectors of :math:`A A^H`. In both cases the corresponding
1690 (possibly non-zero) eigenvalues are given by ``s**2``.
1692 If `a` has more than two dimensions, then broadcasting rules apply, as
1693 explained in :ref:`routines.linalg-broadcasting`. This means that SVD is
1694 working in "stacked" mode: it iterates over all indices of the first
1695 ``a.ndim - 2`` dimensions and for each combination SVD is applied to the
1696 last two indices. The matrix `a` can be reconstructed from the
1697 decomposition with either ``(u * s[..., None, :]) @ vh`` or
1698 ``u @ (s[..., None] * vh)``. (The ``@`` operator can be replaced by the
1699 function ``np.matmul`` for python versions below 3.5.)
1701 If `a` is a ``matrix`` object (as opposed to an ``ndarray``), then so are
1702 all the return values.
1704 Examples
1705 --------
1706 >>> a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6)
1707 >>> b = np.random.randn(2, 7, 8, 3) + 1j*np.random.randn(2, 7, 8, 3)
1709 Reconstruction based on full SVD, 2D case:
1711 >>> U, S, Vh = np.linalg.svd(a, full_matrices=True)
1712 >>> U.shape, S.shape, Vh.shape
1713 ((9, 9), (6,), (6, 6))
1714 >>> np.allclose(a, np.dot(U[:, :6] * S, Vh))
1715 True
1716 >>> smat = np.zeros((9, 6), dtype=complex)
1717 >>> smat[:6, :6] = np.diag(S)
1718 >>> np.allclose(a, np.dot(U, np.dot(smat, Vh)))
1719 True
1721 Reconstruction based on reduced SVD, 2D case:
1723 >>> U, S, Vh = np.linalg.svd(a, full_matrices=False)
1724 >>> U.shape, S.shape, Vh.shape
1725 ((9, 6), (6,), (6, 6))
1726 >>> np.allclose(a, np.dot(U * S, Vh))
1727 True
1728 >>> smat = np.diag(S)
1729 >>> np.allclose(a, np.dot(U, np.dot(smat, Vh)))
1730 True
1732 Reconstruction based on full SVD, 4D case:
1734 >>> U, S, Vh = np.linalg.svd(b, full_matrices=True)
1735 >>> U.shape, S.shape, Vh.shape
1736 ((2, 7, 8, 8), (2, 7, 3), (2, 7, 3, 3))
1737 >>> np.allclose(b, np.matmul(U[..., :3] * S[..., None, :], Vh))
1738 True
1739 >>> np.allclose(b, np.matmul(U[..., :3], S[..., None] * Vh))
1740 True
1742 Reconstruction based on reduced SVD, 4D case:
1744 >>> U, S, Vh = np.linalg.svd(b, full_matrices=False)
1745 >>> U.shape, S.shape, Vh.shape
1746 ((2, 7, 8, 3), (2, 7, 3), (2, 7, 3, 3))
1747 >>> np.allclose(b, np.matmul(U * S[..., None, :], Vh))
1748 True
1749 >>> np.allclose(b, np.matmul(U, S[..., None] * Vh))
1750 True
1752 """
1753 import numpy as _nx
1754 a, wrap = _makearray(a)
1756 if hermitian:
1757 # note: lapack svd returns eigenvalues with s ** 2 sorted descending,
1758 # but eig returns s sorted ascending, so we re-order the eigenvalues
1759 # and related arrays to have the correct order
1760 if compute_uv:
1761 s, u = eigh(a)
1762 sgn = sign(s)
1763 s = abs(s)
1764 sidx = argsort(s)[..., ::-1]
1765 sgn = _nx.take_along_axis(sgn, sidx, axis=-1)
1766 s = _nx.take_along_axis(s, sidx, axis=-1)
1767 u = _nx.take_along_axis(u, sidx[..., None, :], axis=-1)
1768 # singular values are unsigned, move the sign into v
1769 vt = transpose(u * sgn[..., None, :]).conjugate()
1770 return SVDResult(wrap(u), s, wrap(vt))
1771 else:
1772 s = eigvalsh(a)
1773 s = abs(s)
1774 return sort(s)[..., ::-1]
1776 _assert_stacked_2d(a)
1777 t, result_t = _commonType(a)
1779 m, n = a.shape[-2:]
1780 if compute_uv:
1781 if full_matrices:
1782 if m < n:
1783 gufunc = _umath_linalg.svd_m_f
1784 else:
1785 gufunc = _umath_linalg.svd_n_f
1786 else:
1787 if m < n:
1788 gufunc = _umath_linalg.svd_m_s
1789 else:
1790 gufunc = _umath_linalg.svd_n_s
1792 signature = 'D->DdD' if isComplexType(t) else 'd->ddd'
1793 with errstate(call=_raise_linalgerror_svd_nonconvergence,
1794 invalid='call', over='ignore', divide='ignore',
1795 under='ignore'):
1796 u, s, vh = gufunc(a, signature=signature)
1797 u = u.astype(result_t, copy=False)
1798 s = s.astype(_realType(result_t), copy=False)
1799 vh = vh.astype(result_t, copy=False)
1800 return SVDResult(wrap(u), s, wrap(vh))
1801 else:
1802 if m < n:
1803 gufunc = _umath_linalg.svd_m
1804 else:
1805 gufunc = _umath_linalg.svd_n
1807 signature = 'D->d' if isComplexType(t) else 'd->d'
1808 with errstate(call=_raise_linalgerror_svd_nonconvergence,
1809 invalid='call', over='ignore', divide='ignore',
1810 under='ignore'):
1811 s = gufunc(a, signature=signature)
1812 s = s.astype(_realType(result_t), copy=False)
1813 return s
1816def _svdvals_dispatcher(x):
1817 return (x,)
1820@array_function_dispatch(_svdvals_dispatcher)
1821def svdvals(x, /):
1822 """
1823 Returns the singular values of a matrix (or a stack of matrices) ``x``.
1824 When x is a stack of matrices, the function will compute the singular
1825 values for each matrix in the stack.
1827 This function is Array API compatible.
1829 Calling ``np.svdvals(x)`` to get singular values is the same as
1830 ``np.svd(x, compute_uv=False, hermitian=False)``.
1832 Parameters
1833 ----------
1834 x : (..., M, N) array_like
1835 Input array having shape (..., M, N) and whose last two
1836 dimensions form matrices on which to perform singular value
1837 decomposition. Should have a floating-point data type.
1839 Returns
1840 -------
1841 out : ndarray
1842 An array with shape (..., K) that contains the vector(s)
1843 of singular values of length K, where K = min(M, N).
1845 See Also
1846 --------
1847 scipy.linalg.svdvals : Compute singular values of a matrix.
1849 """
1850 return svd(x, compute_uv=False, hermitian=False)
1853def _cond_dispatcher(x, p=None):
1854 return (x,)
1857@array_function_dispatch(_cond_dispatcher)
1858def cond(x, p=None):
1859 """
1860 Compute the condition number of a matrix.
1862 This function is capable of returning the condition number using
1863 one of seven different norms, depending on the value of `p` (see
1864 Parameters below).
1866 Parameters
1867 ----------
1868 x : (..., M, N) array_like
1869 The matrix whose condition number is sought.
1870 p : {None, 1, -1, 2, -2, inf, -inf, 'fro'}, optional
1871 Order of the norm used in the condition number computation:
1873 ===== ============================
1874 p norm for matrices
1875 ===== ============================
1876 None 2-norm, computed directly using the ``SVD``
1877 'fro' Frobenius norm
1878 inf max(sum(abs(x), axis=1))
1879 -inf min(sum(abs(x), axis=1))
1880 1 max(sum(abs(x), axis=0))
1881 -1 min(sum(abs(x), axis=0))
1882 2 2-norm (largest sing. value)
1883 -2 smallest singular value
1884 ===== ============================
1886 inf means the `numpy.inf` object, and the Frobenius norm is
1887 the root-of-sum-of-squares norm.
1889 Returns
1890 -------
1891 c : {float, inf}
1892 The condition number of the matrix. May be infinite.
1894 See Also
1895 --------
1896 numpy.linalg.norm
1898 Notes
1899 -----
1900 The condition number of `x` is defined as the norm of `x` times the
1901 norm of the inverse of `x` [1]_; the norm can be the usual L2-norm
1902 (root-of-sum-of-squares) or one of a number of other matrix norms.
1904 References
1905 ----------
1906 .. [1] G. Strang, *Linear Algebra and Its Applications*, Orlando, FL,
1907 Academic Press, Inc., 1980, pg. 285.
1909 Examples
1910 --------
1911 >>> from numpy import linalg as LA
1912 >>> a = np.array([[1, 0, -1], [0, 1, 0], [1, 0, 1]])
1913 >>> a
1914 array([[ 1, 0, -1],
1915 [ 0, 1, 0],
1916 [ 1, 0, 1]])
1917 >>> LA.cond(a)
1918 1.4142135623730951
1919 >>> LA.cond(a, 'fro')
1920 3.1622776601683795
1921 >>> LA.cond(a, np.inf)
1922 2.0
1923 >>> LA.cond(a, -np.inf)
1924 1.0
1925 >>> LA.cond(a, 1)
1926 2.0
1927 >>> LA.cond(a, -1)
1928 1.0
1929 >>> LA.cond(a, 2)
1930 1.4142135623730951
1931 >>> LA.cond(a, -2)
1932 0.70710678118654746 # may vary
1933 >>> (min(LA.svd(a, compute_uv=False)) *
1934 ... min(LA.svd(LA.inv(a), compute_uv=False)))
1935 0.70710678118654746 # may vary
1937 """
1938 x = asarray(x) # in case we have a matrix
1939 if _is_empty_2d(x):
1940 raise LinAlgError("cond is not defined on empty arrays")
1941 if p is None or p == 2 or p == -2:
1942 s = svd(x, compute_uv=False)
1943 with errstate(all='ignore'):
1944 if p == -2:
1945 r = s[..., -1] / s[..., 0]
1946 else:
1947 r = s[..., 0] / s[..., -1]
1948 else:
1949 # Call inv(x) ignoring errors. The result array will
1950 # contain nans in the entries where inversion failed.
1951 _assert_stacked_2d(x)
1952 _assert_stacked_square(x)
1953 t, result_t = _commonType(x)
1954 signature = 'D->D' if isComplexType(t) else 'd->d'
1955 with errstate(all='ignore'):
1956 invx = _umath_linalg.inv(x, signature=signature)
1957 r = norm(x, p, axis=(-2, -1)) * norm(invx, p, axis=(-2, -1))
1958 r = r.astype(result_t, copy=False)
1960 # Convert nans to infs unless the original array had nan entries
1961 r = asarray(r)
1962 nan_mask = isnan(r)
1963 if nan_mask.any():
1964 nan_mask &= ~isnan(x).any(axis=(-2, -1))
1965 if r.ndim > 0:
1966 r[nan_mask] = inf
1967 elif nan_mask:
1968 r[()] = inf
1970 # Convention is to return scalars instead of 0d arrays
1971 if r.ndim == 0:
1972 r = r[()]
1974 return r
1977def _matrix_rank_dispatcher(A, tol=None, hermitian=None, *, rtol=None):
1978 return (A,)
1981@array_function_dispatch(_matrix_rank_dispatcher)
1982def matrix_rank(A, tol=None, hermitian=False, *, rtol=None):
1983 """
1984 Return matrix rank of array using SVD method
1986 Rank of the array is the number of singular values of the array that are
1987 greater than `tol`.
1989 .. versionchanged:: 1.14
1990 Can now operate on stacks of matrices
1992 Parameters
1993 ----------
1994 A : {(M,), (..., M, N)} array_like
1995 Input vector or stack of matrices.
1996 tol : (...) array_like, float, optional
1997 Threshold below which SVD values are considered zero. If `tol` is
1998 None, and ``S`` is an array with singular values for `M`, and
1999 ``eps`` is the epsilon value for datatype of ``S``, then `tol` is
2000 set to ``S.max() * max(M, N) * eps``.
2002 .. versionchanged:: 1.14
2003 Broadcasted against the stack of matrices
2004 hermitian : bool, optional
2005 If True, `A` is assumed to be Hermitian (symmetric if real-valued),
2006 enabling a more efficient method for finding singular values.
2007 Defaults to False.
2009 .. versionadded:: 1.14
2010 rtol : (...) array_like, float, optional
2011 Parameter for the relative tolerance component. Only ``tol`` or
2012 ``rtol`` can be set at a time. Defaults to ``max(M, N) * eps``.
2014 .. versionadded:: 2.0.0
2016 Returns
2017 -------
2018 rank : (...) array_like
2019 Rank of A.
2021 Notes
2022 -----
2023 The default threshold to detect rank deficiency is a test on the magnitude
2024 of the singular values of `A`. By default, we identify singular values
2025 less than ``S.max() * max(M, N) * eps`` as indicating rank deficiency
2026 (with the symbols defined above). This is the algorithm MATLAB uses [1].
2027 It also appears in *Numerical recipes* in the discussion of SVD solutions
2028 for linear least squares [2].
2030 This default threshold is designed to detect rank deficiency accounting
2031 for the numerical errors of the SVD computation. Imagine that there
2032 is a column in `A` that is an exact (in floating point) linear combination
2033 of other columns in `A`. Computing the SVD on `A` will not produce
2034 a singular value exactly equal to 0 in general: any difference of
2035 the smallest SVD value from 0 will be caused by numerical imprecision
2036 in the calculation of the SVD. Our threshold for small SVD values takes
2037 this numerical imprecision into account, and the default threshold will
2038 detect such numerical rank deficiency. The threshold may declare a matrix
2039 `A` rank deficient even if the linear combination of some columns of `A`
2040 is not exactly equal to another column of `A` but only numerically very
2041 close to another column of `A`.
2043 We chose our default threshold because it is in wide use. Other thresholds
2044 are possible. For example, elsewhere in the 2007 edition of *Numerical
2045 recipes* there is an alternative threshold of ``S.max() *
2046 np.finfo(A.dtype).eps / 2. * np.sqrt(m + n + 1.)``. The authors describe
2047 this threshold as being based on "expected roundoff error" (p 71).
2049 The thresholds above deal with floating point roundoff error in the
2050 calculation of the SVD. However, you may have more information about
2051 the sources of error in `A` that would make you consider other tolerance
2052 values to detect *effective* rank deficiency. The most useful measure
2053 of the tolerance depends on the operations you intend to use on your
2054 matrix. For example, if your data come from uncertain measurements with
2055 uncertainties greater than floating point epsilon, choosing a tolerance
2056 near that uncertainty may be preferable. The tolerance may be absolute
2057 if the uncertainties are absolute rather than relative.
2059 References
2060 ----------
2061 .. [1] MATLAB reference documentation, "Rank"
2062 https://www.mathworks.com/help/techdoc/ref/rank.html
2063 .. [2] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery,
2064 "Numerical Recipes (3rd edition)", Cambridge University Press, 2007,
2065 page 795.
2067 Examples
2068 --------
2069 >>> from numpy.linalg import matrix_rank
2070 >>> matrix_rank(np.eye(4)) # Full rank matrix
2071 4
2072 >>> I=np.eye(4); I[-1,-1] = 0. # rank deficient matrix
2073 >>> matrix_rank(I)
2074 3
2075 >>> matrix_rank(np.ones((4,))) # 1 dimension - rank 1 unless all 0
2076 1
2077 >>> matrix_rank(np.zeros((4,)))
2078 0
2079 """
2080 if rtol is not None and tol is not None:
2081 raise ValueError("`tol` and `rtol` can't be both set.")
2083 A = asarray(A)
2084 if A.ndim < 2:
2085 return int(not all(A == 0))
2086 S = svd(A, compute_uv=False, hermitian=hermitian)
2088 if tol is None:
2089 if rtol is None:
2090 rtol = max(A.shape[-2:]) * finfo(S.dtype).eps
2091 else:
2092 rtol = asarray(rtol)[..., newaxis]
2093 tol = S.max(axis=-1, keepdims=True) * rtol
2094 else:
2095 tol = asarray(tol)[..., newaxis]
2097 return count_nonzero(S > tol, axis=-1)
2100# Generalized inverse
2102def _pinv_dispatcher(a, rcond=None, hermitian=None, *, rtol=None):
2103 return (a,)
2106@array_function_dispatch(_pinv_dispatcher)
2107def pinv(a, rcond=None, hermitian=False, *, rtol=_NoValue):
2108 """
2109 Compute the (Moore-Penrose) pseudo-inverse of a matrix.
2111 Calculate the generalized inverse of a matrix using its
2112 singular-value decomposition (SVD) and including all
2113 *large* singular values.
2115 .. versionchanged:: 1.14
2116 Can now operate on stacks of matrices
2118 Parameters
2119 ----------
2120 a : (..., M, N) array_like
2121 Matrix or stack of matrices to be pseudo-inverted.
2122 rcond : (...) array_like of float, optional
2123 Cutoff for small singular values.
2124 Singular values less than or equal to
2125 ``rcond * largest_singular_value`` are set to zero.
2126 Broadcasts against the stack of matrices. Default: ``1e-15``.
2127 hermitian : bool, optional
2128 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
2129 enabling a more efficient method for finding singular values.
2130 Defaults to False.
2132 .. versionadded:: 1.17.0
2133 rtol : (...) array_like of float, optional
2134 Same as `rcond`, but it's an Array API compatible parameter name.
2135 Only `rcond` or `rtol` can be set at a time. If none of them are
2136 provided then NumPy's ``1e-15`` default is used. If ``rtol=None``
2137 is passed then the API standard default is used.
2139 .. versionadded:: 2.0.0
2141 Returns
2142 -------
2143 B : (..., N, M) ndarray
2144 The pseudo-inverse of `a`. If `a` is a `matrix` instance, then so
2145 is `B`.
2147 Raises
2148 ------
2149 LinAlgError
2150 If the SVD computation does not converge.
2152 See Also
2153 --------
2154 scipy.linalg.pinv : Similar function in SciPy.
2155 scipy.linalg.pinvh : Compute the (Moore-Penrose) pseudo-inverse of a
2156 Hermitian matrix.
2158 Notes
2159 -----
2160 The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
2161 defined as: "the matrix that 'solves' [the least-squares problem]
2162 :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then
2163 :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`.
2165 It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular
2166 value decomposition of A, then
2167 :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are
2168 orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting
2169 of A's so-called singular values, (followed, typically, by
2170 zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix
2171 consisting of the reciprocals of A's singular values
2172 (again, followed by zeros). [1]_
2174 References
2175 ----------
2176 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
2177 FL, Academic Press, Inc., 1980, pp. 139-142.
2179 Examples
2180 --------
2181 The following example checks that ``a * a+ * a == a`` and
2182 ``a+ * a * a+ == a+``:
2184 >>> a = np.random.randn(9, 6)
2185 >>> B = np.linalg.pinv(a)
2186 >>> np.allclose(a, np.dot(a, np.dot(B, a)))
2187 True
2188 >>> np.allclose(B, np.dot(B, np.dot(a, B)))
2189 True
2191 """
2192 a, wrap = _makearray(a)
2193 if rcond is None:
2194 if rtol is _NoValue:
2195 rcond = 1e-15
2196 elif rtol is None:
2197 rcond = max(a.shape[-2:]) * finfo(a.dtype).eps
2198 else:
2199 rcond = rtol
2200 elif rtol is not _NoValue:
2201 raise ValueError("`rtol` and `rcond` can't be both set.")
2202 else:
2203 # NOTE: Deprecate `rcond` in a few versions.
2204 pass
2206 rcond = asarray(rcond)
2207 if _is_empty_2d(a):
2208 m, n = a.shape[-2:]
2209 res = empty(a.shape[:-2] + (n, m), dtype=a.dtype)
2210 return wrap(res)
2211 a = a.conjugate()
2212 u, s, vt = svd(a, full_matrices=False, hermitian=hermitian)
2214 # discard small singular values
2215 cutoff = rcond[..., newaxis] * amax(s, axis=-1, keepdims=True)
2216 large = s > cutoff
2217 s = divide(1, s, where=large, out=s)
2218 s[~large] = 0
2220 res = matmul(transpose(vt), multiply(s[..., newaxis], transpose(u)))
2221 return wrap(res)
2224# Determinant
2227@array_function_dispatch(_unary_dispatcher)
2228def slogdet(a):
2229 """
2230 Compute the sign and (natural) logarithm of the determinant of an array.
2232 If an array has a very small or very large determinant, then a call to
2233 `det` may overflow or underflow. This routine is more robust against such
2234 issues, because it computes the logarithm of the determinant rather than
2235 the determinant itself.
2237 Parameters
2238 ----------
2239 a : (..., M, M) array_like
2240 Input array, has to be a square 2-D array.
2242 Returns
2243 -------
2244 A namedtuple with the following attributes:
2246 sign : (...) array_like
2247 A number representing the sign of the determinant. For a real matrix,
2248 this is 1, 0, or -1. For a complex matrix, this is a complex number
2249 with absolute value 1 (i.e., it is on the unit circle), or else 0.
2250 logabsdet : (...) array_like
2251 The natural log of the absolute value of the determinant.
2253 If the determinant is zero, then `sign` will be 0 and `logabsdet`
2254 will be -inf. In all cases, the determinant is equal to
2255 ``sign * np.exp(logabsdet)``.
2257 See Also
2258 --------
2259 det
2261 Notes
2262 -----
2264 .. versionadded:: 1.8.0
2266 Broadcasting rules apply, see the `numpy.linalg` documentation for
2267 details.
2269 .. versionadded:: 1.6.0
2271 The determinant is computed via LU factorization using the LAPACK
2272 routine ``z/dgetrf``.
2275 Examples
2276 --------
2277 The determinant of a 2-D array ``[[a, b], [c, d]]`` is ``ad - bc``:
2279 >>> a = np.array([[1, 2], [3, 4]])
2280 >>> (sign, logabsdet) = np.linalg.slogdet(a)
2281 >>> (sign, logabsdet)
2282 (-1, 0.69314718055994529) # may vary
2283 >>> sign * np.exp(logabsdet)
2284 -2.0
2286 Computing log-determinants for a stack of matrices:
2288 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2289 >>> a.shape
2290 (3, 2, 2)
2291 >>> sign, logabsdet = np.linalg.slogdet(a)
2292 >>> (sign, logabsdet)
2293 (array([-1., -1., -1.]), array([ 0.69314718, 1.09861229, 2.07944154]))
2294 >>> sign * np.exp(logabsdet)
2295 array([-2., -3., -8.])
2297 This routine succeeds where ordinary `det` does not:
2299 >>> np.linalg.det(np.eye(500) * 0.1)
2300 0.0
2301 >>> np.linalg.slogdet(np.eye(500) * 0.1)
2302 (1, -1151.2925464970228)
2304 """
2305 a = asarray(a)
2306 _assert_stacked_2d(a)
2307 _assert_stacked_square(a)
2308 t, result_t = _commonType(a)
2309 real_t = _realType(result_t)
2310 signature = 'D->Dd' if isComplexType(t) else 'd->dd'
2311 sign, logdet = _umath_linalg.slogdet(a, signature=signature)
2312 sign = sign.astype(result_t, copy=False)
2313 logdet = logdet.astype(real_t, copy=False)
2314 return SlogdetResult(sign, logdet)
2317@array_function_dispatch(_unary_dispatcher)
2318def det(a):
2319 """
2320 Compute the determinant of an array.
2322 Parameters
2323 ----------
2324 a : (..., M, M) array_like
2325 Input array to compute determinants for.
2327 Returns
2328 -------
2329 det : (...) array_like
2330 Determinant of `a`.
2332 See Also
2333 --------
2334 slogdet : Another way to represent the determinant, more suitable
2335 for large matrices where underflow/overflow may occur.
2336 scipy.linalg.det : Similar function in SciPy.
2338 Notes
2339 -----
2341 .. versionadded:: 1.8.0
2343 Broadcasting rules apply, see the `numpy.linalg` documentation for
2344 details.
2346 The determinant is computed via LU factorization using the LAPACK
2347 routine ``z/dgetrf``.
2349 Examples
2350 --------
2351 The determinant of a 2-D array [[a, b], [c, d]] is ad - bc:
2353 >>> a = np.array([[1, 2], [3, 4]])
2354 >>> np.linalg.det(a)
2355 -2.0 # may vary
2357 Computing determinants for a stack of matrices:
2359 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2360 >>> a.shape
2361 (3, 2, 2)
2362 >>> np.linalg.det(a)
2363 array([-2., -3., -8.])
2365 """
2366 a = asarray(a)
2367 _assert_stacked_2d(a)
2368 _assert_stacked_square(a)
2369 t, result_t = _commonType(a)
2370 signature = 'D->D' if isComplexType(t) else 'd->d'
2371 r = _umath_linalg.det(a, signature=signature)
2372 r = r.astype(result_t, copy=False)
2373 return r
2376# Linear Least Squares
2378def _lstsq_dispatcher(a, b, rcond=None):
2379 return (a, b)
2382@array_function_dispatch(_lstsq_dispatcher)
2383def lstsq(a, b, rcond=None):
2384 r"""
2385 Return the least-squares solution to a linear matrix equation.
2387 Computes the vector `x` that approximately solves the equation
2388 ``a @ x = b``. The equation may be under-, well-, or over-determined
2389 (i.e., the number of linearly independent rows of `a` can be less than,
2390 equal to, or greater than its number of linearly independent columns).
2391 If `a` is square and of full rank, then `x` (but for round-off error)
2392 is the "exact" solution of the equation. Else, `x` minimizes the
2393 Euclidean 2-norm :math:`||b - ax||`. If there are multiple minimizing
2394 solutions, the one with the smallest 2-norm :math:`||x||` is returned.
2396 Parameters
2397 ----------
2398 a : (M, N) array_like
2399 "Coefficient" matrix.
2400 b : {(M,), (M, K)} array_like
2401 Ordinate or "dependent variable" values. If `b` is two-dimensional,
2402 the least-squares solution is calculated for each of the `K` columns
2403 of `b`.
2404 rcond : float, optional
2405 Cut-off ratio for small singular values of `a`.
2406 For the purposes of rank determination, singular values are treated
2407 as zero if they are smaller than `rcond` times the largest singular
2408 value of `a`.
2409 The default uses the machine precision times ``max(M, N)``. Passing
2410 ``-1`` will use machine precision.
2412 .. versionchanged:: 2.0
2413 Previously, the default was ``-1``, but a warning was given that
2414 this would change.
2416 Returns
2417 -------
2418 x : {(N,), (N, K)} ndarray
2419 Least-squares solution. If `b` is two-dimensional,
2420 the solutions are in the `K` columns of `x`.
2421 residuals : {(1,), (K,), (0,)} ndarray
2422 Sums of squared residuals: Squared Euclidean 2-norm for each column in
2423 ``b - a @ x``.
2424 If the rank of `a` is < N or M <= N, this is an empty array.
2425 If `b` is 1-dimensional, this is a (1,) shape array.
2426 Otherwise the shape is (K,).
2427 rank : int
2428 Rank of matrix `a`.
2429 s : (min(M, N),) ndarray
2430 Singular values of `a`.
2432 Raises
2433 ------
2434 LinAlgError
2435 If computation does not converge.
2437 See Also
2438 --------
2439 scipy.linalg.lstsq : Similar function in SciPy.
2441 Notes
2442 -----
2443 If `b` is a matrix, then all array results are returned as matrices.
2445 Examples
2446 --------
2447 Fit a line, ``y = mx + c``, through some noisy data-points:
2449 >>> x = np.array([0, 1, 2, 3])
2450 >>> y = np.array([-1, 0.2, 0.9, 2.1])
2452 By examining the coefficients, we see that the line should have a
2453 gradient of roughly 1 and cut the y-axis at, more or less, -1.
2455 We can rewrite the line equation as ``y = Ap``, where ``A = [[x 1]]``
2456 and ``p = [[m], [c]]``. Now use `lstsq` to solve for `p`:
2458 >>> A = np.vstack([x, np.ones(len(x))]).T
2459 >>> A
2460 array([[ 0., 1.],
2461 [ 1., 1.],
2462 [ 2., 1.],
2463 [ 3., 1.]])
2465 >>> m, c = np.linalg.lstsq(A, y)[0]
2466 >>> m, c
2467 (1.0 -0.95) # may vary
2469 Plot the data along with the fitted line:
2471 >>> import matplotlib.pyplot as plt
2472 >>> _ = plt.plot(x, y, 'o', label='Original data', markersize=10)
2473 >>> _ = plt.plot(x, m*x + c, 'r', label='Fitted line')
2474 >>> _ = plt.legend()
2475 >>> plt.show()
2477 """
2478 a, _ = _makearray(a)
2479 b, wrap = _makearray(b)
2480 is_1d = b.ndim == 1
2481 if is_1d:
2482 b = b[:, newaxis]
2483 _assert_2d(a, b)
2484 m, n = a.shape[-2:]
2485 m2, n_rhs = b.shape[-2:]
2486 if m != m2:
2487 raise LinAlgError('Incompatible dimensions')
2489 t, result_t = _commonType(a, b)
2490 result_real_t = _realType(result_t)
2492 if rcond is None:
2493 rcond = finfo(t).eps * max(n, m)
2495 if m <= n:
2496 gufunc = _umath_linalg.lstsq_m
2497 else:
2498 gufunc = _umath_linalg.lstsq_n
2500 signature = 'DDd->Ddid' if isComplexType(t) else 'ddd->ddid'
2501 if n_rhs == 0:
2502 # lapack can't handle n_rhs = 0 - so allocate
2503 # the array one larger in that axis
2504 b = zeros(b.shape[:-2] + (m, n_rhs + 1), dtype=b.dtype)
2506 with errstate(call=_raise_linalgerror_lstsq, invalid='call',
2507 over='ignore', divide='ignore', under='ignore'):
2508 x, resids, rank, s = gufunc(a, b, rcond, signature=signature)
2509 if m == 0:
2510 x[...] = 0
2511 if n_rhs == 0:
2512 # remove the item we added
2513 x = x[..., :n_rhs]
2514 resids = resids[..., :n_rhs]
2516 # remove the axis we added
2517 if is_1d:
2518 x = x.squeeze(axis=-1)
2519 # we probably should squeeze resids too, but we can't
2520 # without breaking compatibility.
2522 # as documented
2523 if rank != n or m <= n:
2524 resids = array([], result_real_t)
2526 # coerce output arrays
2527 s = s.astype(result_real_t, copy=False)
2528 resids = resids.astype(result_real_t, copy=False)
2529 # Copying lets the memory in r_parts be freed
2530 x = x.astype(result_t, copy=True)
2531 return wrap(x), wrap(resids), rank, s
2534def _multi_svd_norm(x, row_axis, col_axis, op):
2535 """Compute a function of the singular values of the 2-D matrices in `x`.
2537 This is a private utility function used by `numpy.linalg.norm()`.
2539 Parameters
2540 ----------
2541 x : ndarray
2542 row_axis, col_axis : int
2543 The axes of `x` that hold the 2-D matrices.
2544 op : callable
2545 This should be either numpy.amin or `numpy.amax` or `numpy.sum`.
2547 Returns
2548 -------
2549 result : float or ndarray
2550 If `x` is 2-D, the return values is a float.
2551 Otherwise, it is an array with ``x.ndim - 2`` dimensions.
2552 The return values are either the minimum or maximum or sum of the
2553 singular values of the matrices, depending on whether `op`
2554 is `numpy.amin` or `numpy.amax` or `numpy.sum`.
2556 """
2557 y = moveaxis(x, (row_axis, col_axis), (-2, -1))
2558 result = op(svd(y, compute_uv=False), axis=-1)
2559 return result
2562def _norm_dispatcher(x, ord=None, axis=None, keepdims=None):
2563 return (x,)
2566@array_function_dispatch(_norm_dispatcher)
2567def norm(x, ord=None, axis=None, keepdims=False):
2568 """
2569 Matrix or vector norm.
2571 This function is able to return one of eight different matrix norms,
2572 or one of an infinite number of vector norms (described below), depending
2573 on the value of the ``ord`` parameter.
2575 Parameters
2576 ----------
2577 x : array_like
2578 Input array. If `axis` is None, `x` must be 1-D or 2-D, unless `ord`
2579 is None. If both `axis` and `ord` are None, the 2-norm of
2580 ``x.ravel`` will be returned.
2581 ord : {non-zero int, inf, -inf, 'fro', 'nuc'}, optional
2582 Order of the norm (see table under ``Notes``). inf means numpy's
2583 `inf` object. The default is None.
2584 axis : {None, int, 2-tuple of ints}, optional.
2585 If `axis` is an integer, it specifies the axis of `x` along which to
2586 compute the vector norms. If `axis` is a 2-tuple, it specifies the
2587 axes that hold 2-D matrices, and the matrix norms of these matrices
2588 are computed. If `axis` is None then either a vector norm (when `x`
2589 is 1-D) or a matrix norm (when `x` is 2-D) is returned. The default
2590 is None.
2592 .. versionadded:: 1.8.0
2594 keepdims : bool, optional
2595 If this is set to True, the axes which are normed over are left in the
2596 result as dimensions with size one. With this option the result will
2597 broadcast correctly against the original `x`.
2599 .. versionadded:: 1.10.0
2601 Returns
2602 -------
2603 n : float or ndarray
2604 Norm of the matrix or vector(s).
2606 See Also
2607 --------
2608 scipy.linalg.norm : Similar function in SciPy.
2610 Notes
2611 -----
2612 For values of ``ord < 1``, the result is, strictly speaking, not a
2613 mathematical 'norm', but it may still be useful for various numerical
2614 purposes.
2616 The following norms can be calculated:
2618 ===== ============================ ==========================
2619 ord norm for matrices norm for vectors
2620 ===== ============================ ==========================
2621 None Frobenius norm 2-norm
2622 'fro' Frobenius norm --
2623 'nuc' nuclear norm --
2624 inf max(sum(abs(x), axis=1)) max(abs(x))
2625 -inf min(sum(abs(x), axis=1)) min(abs(x))
2626 0 -- sum(x != 0)
2627 1 max(sum(abs(x), axis=0)) as below
2628 -1 min(sum(abs(x), axis=0)) as below
2629 2 2-norm (largest sing. value) as below
2630 -2 smallest singular value as below
2631 other -- sum(abs(x)**ord)**(1./ord)
2632 ===== ============================ ==========================
2634 The Frobenius norm is given by [1]_:
2636 :math:`||A||_F = [\\sum_{i,j} abs(a_{i,j})^2]^{1/2}`
2638 The nuclear norm is the sum of the singular values.
2640 Both the Frobenius and nuclear norm orders are only defined for
2641 matrices and raise a ValueError when ``x.ndim != 2``.
2643 References
2644 ----------
2645 .. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*,
2646 Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15
2648 Examples
2649 --------
2650 >>> from numpy import linalg as LA
2651 >>> a = np.arange(9) - 4
2652 >>> a
2653 array([-4, -3, -2, ..., 2, 3, 4])
2654 >>> b = a.reshape((3, 3))
2655 >>> b
2656 array([[-4, -3, -2],
2657 [-1, 0, 1],
2658 [ 2, 3, 4]])
2660 >>> LA.norm(a)
2661 7.745966692414834
2662 >>> LA.norm(b)
2663 7.745966692414834
2664 >>> LA.norm(b, 'fro')
2665 7.745966692414834
2666 >>> LA.norm(a, np.inf)
2667 4.0
2668 >>> LA.norm(b, np.inf)
2669 9.0
2670 >>> LA.norm(a, -np.inf)
2671 0.0
2672 >>> LA.norm(b, -np.inf)
2673 2.0
2675 >>> LA.norm(a, 1)
2676 20.0
2677 >>> LA.norm(b, 1)
2678 7.0
2679 >>> LA.norm(a, -1)
2680 -4.6566128774142013e-010
2681 >>> LA.norm(b, -1)
2682 6.0
2683 >>> LA.norm(a, 2)
2684 7.745966692414834
2685 >>> LA.norm(b, 2)
2686 7.3484692283495345
2688 >>> LA.norm(a, -2)
2689 0.0
2690 >>> LA.norm(b, -2)
2691 1.8570331885190563e-016 # may vary
2692 >>> LA.norm(a, 3)
2693 5.8480354764257312 # may vary
2694 >>> LA.norm(a, -3)
2695 0.0
2697 Using the `axis` argument to compute vector norms:
2699 >>> c = np.array([[ 1, 2, 3],
2700 ... [-1, 1, 4]])
2701 >>> LA.norm(c, axis=0)
2702 array([ 1.41421356, 2.23606798, 5. ])
2703 >>> LA.norm(c, axis=1)
2704 array([ 3.74165739, 4.24264069])
2705 >>> LA.norm(c, ord=1, axis=1)
2706 array([ 6., 6.])
2708 Using the `axis` argument to compute matrix norms:
2710 >>> m = np.arange(8).reshape(2,2,2)
2711 >>> LA.norm(m, axis=(1,2))
2712 array([ 3.74165739, 11.22497216])
2713 >>> LA.norm(m[0, :, :]), LA.norm(m[1, :, :])
2714 (3.7416573867739413, 11.224972160321824)
2716 """
2717 x = asarray(x)
2719 if not issubclass(x.dtype.type, (inexact, object_)):
2720 x = x.astype(float)
2722 # Immediately handle some default, simple, fast, and common cases.
2723 if axis is None:
2724 ndim = x.ndim
2725 if (
2726 (ord is None) or
2727 (ord in ('f', 'fro') and ndim == 2) or
2728 (ord == 2 and ndim == 1)
2729 ):
2730 x = x.ravel(order='K')
2731 if isComplexType(x.dtype.type):
2732 x_real = x.real
2733 x_imag = x.imag
2734 sqnorm = x_real.dot(x_real) + x_imag.dot(x_imag)
2735 else:
2736 sqnorm = x.dot(x)
2737 ret = sqrt(sqnorm)
2738 if keepdims:
2739 ret = ret.reshape(ndim*[1])
2740 return ret
2742 # Normalize the `axis` argument to a tuple.
2743 nd = x.ndim
2744 if axis is None:
2745 axis = tuple(range(nd))
2746 elif not isinstance(axis, tuple):
2747 try:
2748 axis = int(axis)
2749 except Exception as e:
2750 raise TypeError(
2751 "'axis' must be None, an integer or a tuple of integers"
2752 ) from e
2753 axis = (axis,)
2755 if len(axis) == 1:
2756 if ord == inf:
2757 return abs(x).max(axis=axis, keepdims=keepdims)
2758 elif ord == -inf:
2759 return abs(x).min(axis=axis, keepdims=keepdims)
2760 elif ord == 0:
2761 # Zero norm
2762 return (
2763 (x != 0)
2764 .astype(x.real.dtype)
2765 .sum(axis=axis, keepdims=keepdims)
2766 )
2767 elif ord == 1:
2768 # special case for speedup
2769 return add.reduce(abs(x), axis=axis, keepdims=keepdims)
2770 elif ord is None or ord == 2:
2771 # special case for speedup
2772 s = (x.conj() * x).real
2773 return sqrt(add.reduce(s, axis=axis, keepdims=keepdims))
2774 # None of the str-type keywords for ord ('fro', 'nuc')
2775 # are valid for vectors
2776 elif isinstance(ord, str):
2777 raise ValueError(f"Invalid norm order '{ord}' for vectors")
2778 else:
2779 absx = abs(x)
2780 absx **= ord
2781 ret = add.reduce(absx, axis=axis, keepdims=keepdims)
2782 ret **= reciprocal(ord, dtype=ret.dtype)
2783 return ret
2784 elif len(axis) == 2:
2785 row_axis, col_axis = axis
2786 row_axis = normalize_axis_index(row_axis, nd)
2787 col_axis = normalize_axis_index(col_axis, nd)
2788 if row_axis == col_axis:
2789 raise ValueError('Duplicate axes given.')
2790 if ord == 2:
2791 ret = _multi_svd_norm(x, row_axis, col_axis, amax)
2792 elif ord == -2:
2793 ret = _multi_svd_norm(x, row_axis, col_axis, amin)
2794 elif ord == 1:
2795 if col_axis > row_axis:
2796 col_axis -= 1
2797 ret = add.reduce(abs(x), axis=row_axis).max(axis=col_axis)
2798 elif ord == inf:
2799 if row_axis > col_axis:
2800 row_axis -= 1
2801 ret = add.reduce(abs(x), axis=col_axis).max(axis=row_axis)
2802 elif ord == -1:
2803 if col_axis > row_axis:
2804 col_axis -= 1
2805 ret = add.reduce(abs(x), axis=row_axis).min(axis=col_axis)
2806 elif ord == -inf:
2807 if row_axis > col_axis:
2808 row_axis -= 1
2809 ret = add.reduce(abs(x), axis=col_axis).min(axis=row_axis)
2810 elif ord in [None, 'fro', 'f']:
2811 ret = sqrt(add.reduce((x.conj() * x).real, axis=axis))
2812 elif ord == 'nuc':
2813 ret = _multi_svd_norm(x, row_axis, col_axis, sum)
2814 else:
2815 raise ValueError("Invalid norm order for matrices.")
2816 if keepdims:
2817 ret_shape = list(x.shape)
2818 ret_shape[axis[0]] = 1
2819 ret_shape[axis[1]] = 1
2820 ret = ret.reshape(ret_shape)
2821 return ret
2822 else:
2823 raise ValueError("Improper number of dimensions to norm.")
2826# multi_dot
2828def _multidot_dispatcher(arrays, *, out=None):
2829 yield from arrays
2830 yield out
2833@array_function_dispatch(_multidot_dispatcher)
2834def multi_dot(arrays, *, out=None):
2835 """
2836 Compute the dot product of two or more arrays in a single function call,
2837 while automatically selecting the fastest evaluation order.
2839 `multi_dot` chains `numpy.dot` and uses optimal parenthesization
2840 of the matrices [1]_ [2]_. Depending on the shapes of the matrices,
2841 this can speed up the multiplication a lot.
2843 If the first argument is 1-D it is treated as a row vector.
2844 If the last argument is 1-D it is treated as a column vector.
2845 The other arguments must be 2-D.
2847 Think of `multi_dot` as::
2849 def multi_dot(arrays): return functools.reduce(np.dot, arrays)
2852 Parameters
2853 ----------
2854 arrays : sequence of array_like
2855 If the first argument is 1-D it is treated as row vector.
2856 If the last argument is 1-D it is treated as column vector.
2857 The other arguments must be 2-D.
2858 out : ndarray, optional
2859 Output argument. This must have the exact kind that would be returned
2860 if it was not used. In particular, it must have the right type, must be
2861 C-contiguous, and its dtype must be the dtype that would be returned
2862 for `dot(a, b)`. This is a performance feature. Therefore, if these
2863 conditions are not met, an exception is raised, instead of attempting
2864 to be flexible.
2866 .. versionadded:: 1.19.0
2868 Returns
2869 -------
2870 output : ndarray
2871 Returns the dot product of the supplied arrays.
2873 See Also
2874 --------
2875 numpy.dot : dot multiplication with two arguments.
2877 References
2878 ----------
2880 .. [1] Cormen, "Introduction to Algorithms", Chapter 15.2, p. 370-378
2881 .. [2] https://en.wikipedia.org/wiki/Matrix_chain_multiplication
2883 Examples
2884 --------
2885 `multi_dot` allows you to write::
2887 >>> from numpy.linalg import multi_dot
2888 >>> # Prepare some data
2889 >>> A = np.random.random((10000, 100))
2890 >>> B = np.random.random((100, 1000))
2891 >>> C = np.random.random((1000, 5))
2892 >>> D = np.random.random((5, 333))
2893 >>> # the actual dot multiplication
2894 >>> _ = multi_dot([A, B, C, D])
2896 instead of::
2898 >>> _ = np.dot(np.dot(np.dot(A, B), C), D)
2899 >>> # or
2900 >>> _ = A.dot(B).dot(C).dot(D)
2902 Notes
2903 -----
2904 The cost for a matrix multiplication can be calculated with the
2905 following function::
2907 def cost(A, B):
2908 return A.shape[0] * A.shape[1] * B.shape[1]
2910 Assume we have three matrices
2911 :math:`A_{10x100}, B_{100x5}, C_{5x50}`.
2913 The costs for the two different parenthesizations are as follows::
2915 cost((AB)C) = 10*100*5 + 10*5*50 = 5000 + 2500 = 7500
2916 cost(A(BC)) = 10*100*50 + 100*5*50 = 50000 + 25000 = 75000
2918 """
2919 n = len(arrays)
2920 # optimization only makes sense for len(arrays) > 2
2921 if n < 2:
2922 raise ValueError("Expecting at least two arrays.")
2923 elif n == 2:
2924 return dot(arrays[0], arrays[1], out=out)
2926 arrays = [asanyarray(a) for a in arrays]
2928 # save original ndim to reshape the result array into the proper form later
2929 ndim_first, ndim_last = arrays[0].ndim, arrays[-1].ndim
2930 # Explicitly convert vectors to 2D arrays to keep the logic of the internal
2931 # _multi_dot_* functions as simple as possible.
2932 if arrays[0].ndim == 1:
2933 arrays[0] = atleast_2d(arrays[0])
2934 if arrays[-1].ndim == 1:
2935 arrays[-1] = atleast_2d(arrays[-1]).T
2936 _assert_2d(*arrays)
2938 # _multi_dot_three is much faster than _multi_dot_matrix_chain_order
2939 if n == 3:
2940 result = _multi_dot_three(arrays[0], arrays[1], arrays[2], out=out)
2941 else:
2942 order = _multi_dot_matrix_chain_order(arrays)
2943 result = _multi_dot(arrays, order, 0, n - 1, out=out)
2945 # return proper shape
2946 if ndim_first == 1 and ndim_last == 1:
2947 return result[0, 0] # scalar
2948 elif ndim_first == 1 or ndim_last == 1:
2949 return result.ravel() # 1-D
2950 else:
2951 return result
2954def _multi_dot_three(A, B, C, out=None):
2955 """
2956 Find the best order for three arrays and do the multiplication.
2958 For three arguments `_multi_dot_three` is approximately 15 times faster
2959 than `_multi_dot_matrix_chain_order`
2961 """
2962 a0, a1b0 = A.shape
2963 b1c0, c1 = C.shape
2964 # cost1 = cost((AB)C) = a0*a1b0*b1c0 + a0*b1c0*c1
2965 cost1 = a0 * b1c0 * (a1b0 + c1)
2966 # cost2 = cost(A(BC)) = a1b0*b1c0*c1 + a0*a1b0*c1
2967 cost2 = a1b0 * c1 * (a0 + b1c0)
2969 if cost1 < cost2:
2970 return dot(dot(A, B), C, out=out)
2971 else:
2972 return dot(A, dot(B, C), out=out)
2975def _multi_dot_matrix_chain_order(arrays, return_costs=False):
2976 """
2977 Return a np.array that encodes the optimal order of mutiplications.
2979 The optimal order array is then used by `_multi_dot()` to do the
2980 multiplication.
2982 Also return the cost matrix if `return_costs` is `True`
2984 The implementation CLOSELY follows Cormen, "Introduction to Algorithms",
2985 Chapter 15.2, p. 370-378. Note that Cormen uses 1-based indices.
2987 cost[i, j] = min([
2988 cost[prefix] + cost[suffix] + cost_mult(prefix, suffix)
2989 for k in range(i, j)])
2991 """
2992 n = len(arrays)
2993 # p stores the dimensions of the matrices
2994 # Example for p: A_{10x100}, B_{100x5}, C_{5x50} --> p = [10, 100, 5, 50]
2995 p = [a.shape[0] for a in arrays] + [arrays[-1].shape[1]]
2996 # m is a matrix of costs of the subproblems
2997 # m[i,j]: min number of scalar multiplications needed to compute A_{i..j}
2998 m = zeros((n, n), dtype=double)
2999 # s is the actual ordering
3000 # s[i, j] is the value of k at which we split the product A_i..A_j
3001 s = empty((n, n), dtype=intp)
3003 for l in range(1, n):
3004 for i in range(n - l):
3005 j = i + l
3006 m[i, j] = inf
3007 for k in range(i, j):
3008 q = m[i, k] + m[k+1, j] + p[i]*p[k+1]*p[j+1]
3009 if q < m[i, j]:
3010 m[i, j] = q
3011 s[i, j] = k # Note that Cormen uses 1-based index
3013 return (s, m) if return_costs else s
3016def _multi_dot(arrays, order, i, j, out=None):
3017 """Actually do the multiplication with the given order."""
3018 if i == j:
3019 # the initial call with non-None out should never get here
3020 assert out is None
3022 return arrays[i]
3023 else:
3024 return dot(_multi_dot(arrays, order, i, order[i, j]),
3025 _multi_dot(arrays, order, order[i, j] + 1, j),
3026 out=out)
3029# diagonal
3031def _diagonal_dispatcher(x, /, *, offset=None):
3032 return (x,)
3035@array_function_dispatch(_diagonal_dispatcher)
3036def diagonal(x, /, *, offset=0):
3037 """
3038 Returns specified diagonals of a matrix (or a stack of matrices) ``x``.
3040 This function is Array API compatible, contrary to
3041 :py:func:`numpy.diagonal`, the matrix is assumed
3042 to be defined by the last two dimensions.
3044 Parameters
3045 ----------
3046 x : (...,M,N) array_like
3047 Input array having shape (..., M, N) and whose innermost two
3048 dimensions form MxN matrices.
3049 offset : int, optional
3050 Offset specifying the off-diagonal relative to the main diagonal,
3051 where::
3053 * offset = 0: the main diagonal.
3054 * offset > 0: off-diagonal above the main diagonal.
3055 * offset < 0: off-diagonal below the main diagonal.
3057 Returns
3058 -------
3059 out : (...,min(N,M)) ndarray
3060 An array containing the diagonals and whose shape is determined by
3061 removing the last two dimensions and appending a dimension equal to
3062 the size of the resulting diagonals. The returned array must have
3063 the same data type as ``x``.
3065 See Also
3066 --------
3067 numpy.diagonal
3069 """
3070 return _core_diagonal(x, offset, axis1=-2, axis2=-1)
3073# trace
3075def _trace_dispatcher(x, /, *, offset=None, dtype=None):
3076 return (x,)
3079@array_function_dispatch(_trace_dispatcher)
3080def trace(x, /, *, offset=0, dtype=None):
3081 """
3082 Returns the sum along the specified diagonals of a matrix
3083 (or a stack of matrices) ``x``.
3085 This function is Array API compatible, contrary to
3086 :py:func:`numpy.trace`.
3088 Parameters
3089 ----------
3090 x : (...,M,N) array_like
3091 Input array having shape (..., M, N) and whose innermost two
3092 dimensions form MxN matrices.
3093 offset : int, optional
3094 Offset specifying the off-diagonal relative to the main diagonal,
3095 where::
3097 * offset = 0: the main diagonal.
3098 * offset > 0: off-diagonal above the main diagonal.
3099 * offset < 0: off-diagonal below the main diagonal.
3101 dtype : dtype, optional
3102 Data type of the returned array.
3104 Returns
3105 -------
3106 out : ndarray
3107 An array containing the traces and whose shape is determined by
3108 removing the last two dimensions and storing the traces in the last
3109 array dimension. For example, if x has rank k and shape:
3110 (I, J, K, ..., L, M, N), then an output array has rank k-2 and shape:
3111 (I, J, K, ..., L) where::
3113 out[i, j, k, ..., l] = trace(a[i, j, k, ..., l, :, :])
3115 The returned array must have a data type as described by the dtype
3116 parameter above.
3118 See Also
3119 --------
3120 numpy.trace
3122 """
3123 return _core_trace(x, offset, axis1=-2, axis2=-1, dtype=dtype)
3126# cross
3128def _cross_dispatcher(x1, x2, /, *, axis=None):
3129 return (x1, x2,)
3132@array_function_dispatch(_cross_dispatcher)
3133def cross(x1, x2, /, *, axis=-1):
3134 """
3135 Returns the cross product of 3-element vectors.
3137 If ``x1`` and/or ``x2`` are multi-dimensional arrays, then
3138 the cross-product of each pair of corresponding 3-element vectors
3139 is independently computed.
3141 This function is Array API compatible, contrary to
3142 :func:`numpy.cross`.
3144 Parameters
3145 ----------
3146 x1 : array_like
3147 The first input array.
3148 x2 : array_like
3149 The second input array. Must be compatible with ``x1`` for all
3150 non-compute axes. The size of the axis over which to compute
3151 the cross-product must be the same size as the respective axis
3152 in ``x1``.
3153 axis : int, optional
3154 The axis (dimension) of ``x1`` and ``x2`` containing the vectors for
3155 which to compute the cross-product. Default: ``-1``.
3157 Returns
3158 -------
3159 out : ndarray
3160 An array containing the cross products.
3162 See Also
3163 --------
3164 numpy.cross
3166 """
3167 if x1.shape[axis] != 3 or x2.shape[axis] != 3:
3168 raise ValueError(
3169 "Both input arrays must be (arrays of) 3-dimensional vectors, "
3170 f"but they are {x1.shape[axis]} and {x2.shape[axis]} "
3171 "dimensional instead."
3172 )
3174 return _core_cross(x1, x2, axis=axis)
3177# matmul
3179def _matmul_dispatcher(x1, x2, /):
3180 return (x1, x2)
3183@array_function_dispatch(_matmul_dispatcher)
3184def matmul(x1, x2, /):
3185 """
3186 Computes the matrix product.
3188 This function is Array API compatible, contrary to
3189 :func:`numpy.matmul`.
3191 Parameters
3192 ----------
3193 x1 : array_like
3194 The first input array.
3195 x2 : array_like
3196 The second input array.
3198 Returns
3199 -------
3200 out : ndarray
3201 The matrix product of the inputs.
3202 This is a scalar only when both ``x1``, ``x2`` are 1-d vectors.
3204 Raises
3205 ------
3206 ValueError
3207 If the last dimension of ``x1`` is not the same size as
3208 the second-to-last dimension of ``x2``.
3210 If a scalar value is passed in.
3212 See Also
3213 --------
3214 numpy.matmul
3216 """
3217 return _core_matmul(x1, x2)
3220# tensordot
3222def _tensordot_dispatcher(x1, x2, /, *, axes=None):
3223 return (x1, x2)
3226@array_function_dispatch(_tensordot_dispatcher)
3227def tensordot(x1, x2, /, *, axes=2):
3228 return _core_tensordot(x1, x2, axes=axes)
3231tensordot.__doc__ = _core_tensordot.__doc__
3234# matrix_transpose
3236def _matrix_transpose_dispatcher(x):
3237 return (x,)
3239@array_function_dispatch(_matrix_transpose_dispatcher)
3240def matrix_transpose(x, /):
3241 return _core_matrix_transpose(x)
3244matrix_transpose.__doc__ = _core_matrix_transpose.__doc__
3247# matrix_norm
3249def _matrix_norm_dispatcher(x, /, *, keepdims=None, ord=None):
3250 return (x,)
3252@array_function_dispatch(_matrix_norm_dispatcher)
3253def matrix_norm(x, /, *, keepdims=False, ord="fro"):
3254 """
3255 Computes the matrix norm of a matrix (or a stack of matrices) ``x``.
3257 This function is Array API compatible.
3259 Parameters
3260 ----------
3261 x : array_like
3262 Input array having shape (..., M, N) and whose two innermost
3263 dimensions form ``MxN`` matrices.
3264 keepdims : bool, optional
3265 If this is set to True, the axes which are normed over are left in
3266 the result as dimensions with size one. Default: False.
3267 ord : {1, -1, 2, -2, inf, -inf, 'fro', 'nuc'}, optional
3268 The order of the norm. For details see the table under ``Notes``
3269 in `numpy.linalg.norm`.
3271 See Also
3272 --------
3273 numpy.linalg.norm : Generic norm function
3275 """
3276 x = asanyarray(x)
3277 return norm(x, axis=(-2, -1), keepdims=keepdims, ord=ord)
3280# vector_norm
3282def _vector_norm_dispatcher(x, /, *, axis=None, keepdims=None, ord=None):
3283 return (x,)
3285@array_function_dispatch(_vector_norm_dispatcher)
3286def vector_norm(x, /, *, axis=None, keepdims=False, ord=2):
3287 """
3288 Computes the vector norm of a vector (or batch of vectors) ``x``.
3290 This function is Array API compatible.
3292 Parameters
3293 ----------
3294 x : array_like
3295 Input array.
3296 axis : {None, int, 2-tuple of ints}, optional
3297 If an integer, ``axis`` specifies the axis (dimension) along which
3298 to compute vector norms. If an n-tuple, ``axis`` specifies the axes
3299 (dimensions) along which to compute batched vector norms. If ``None``,
3300 the vector norm must be computed over all array values (i.e.,
3301 equivalent to computing the vector norm of a flattened array).
3302 Default: ``None``.
3303 keepdims : bool, optional
3304 If this is set to True, the axes which are normed over are left in
3305 the result as dimensions with size one. Default: False.
3306 ord : {1, -1, 2, -2, inf, -inf, 'fro', 'nuc'}, optional
3307 The order of the norm. For details see the table under ``Notes``
3308 in `numpy.linalg.norm`.
3310 See Also
3311 --------
3312 numpy.linalg.norm : Generic norm function
3314 """
3315 x = asanyarray(x)
3316 shape = list(x.shape)
3317 if axis is None:
3318 # Note: np.linalg.norm() doesn't handle 0-D arrays
3319 x = x.ravel()
3320 _axis = 0
3321 elif isinstance(axis, tuple):
3322 # Note: The axis argument supports any number of axes, whereas
3323 # np.linalg.norm() only supports a single axis for vector norm.
3324 normalized_axis = normalize_axis_tuple(axis, x.ndim)
3325 rest = tuple(i for i in range(x.ndim) if i not in normalized_axis)
3326 newshape = axis + rest
3327 x = _core_transpose(x, newshape).reshape(
3328 (
3329 prod([x.shape[i] for i in axis], dtype=int),
3330 *[x.shape[i] for i in rest]
3331 )
3332 )
3333 _axis = 0
3334 else:
3335 _axis = axis
3337 res = norm(x, axis=_axis, ord=ord)
3339 if keepdims:
3340 # We can't reuse np.linalg.norm(keepdims) because of the reshape hacks
3341 # above to avoid matrix norm logic.
3342 _axis = normalize_axis_tuple(
3343 range(len(shape)) if axis is None else axis, len(shape)
3344 )
3345 for i in _axis:
3346 shape[i] = 1
3347 res = res.reshape(tuple(shape))
3349 return res
3352# vecdot
3354def _vecdot_dispatcher(x1, x2, /, *, axis=None):
3355 return (x1, x2)
3357@array_function_dispatch(_vecdot_dispatcher)
3358def vecdot(x1, x2, /, *, axis=-1):
3359 """
3360 Computes the vector dot product.
3362 This function is restricted to arguments compatible with the Array API,
3363 contrary to :func:`numpy.vecdot`.
3365 Let :math:`\\mathbf{a}` be a vector in ``x1`` and :math:`\\mathbf{b}` be
3366 a corresponding vector in ``x2``. The dot product is defined as:
3368 .. math::
3369 \\mathbf{a} \\cdot \\mathbf{b} = \\sum_{i=0}^{n-1} \\overline{a_i}b_i
3371 over the dimension specified by ``axis`` and where :math:`\\overline{a_i}`
3372 denotes the complex conjugate if :math:`a_i` is complex and the identity
3373 otherwise.
3375 Parameters
3376 ----------
3377 x1 : array_like
3378 First input array.
3379 x2 : array_like
3380 Second input array.
3381 axis : int, optional
3382 Axis over which to compute the dot product. Default: ``-1``.
3384 Returns
3385 -------
3386 output : ndarray
3387 The vector dot product of the input.
3389 See Also
3390 --------
3391 numpy.vecdot
3393 """
3394 return _core_vecdot(x1, x2, axis=axis)