Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.10/site-packages/numpy/linalg/_linalg.py: 20%
Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""Lite version of scipy.linalg.
3Notes
4-----
5This module is a lite version of the linalg.py module in SciPy which
6contains high-level Python interface to the LAPACK library. The lite
7version only accesses the following LAPACK functions: dgesv, zgesv,
8dgeev, zgeev, dgesdd, zgesdd, dgelsd, zgelsd, dsyevd, zheevd, dgetrf,
9zgetrf, dpotrf, zpotrf, dgeqrf, zgeqrf, zungqr, dorgqr.
10"""
12__all__ = ['matrix_power', 'solve', 'tensorsolve', 'tensorinv', 'inv',
13 'cholesky', 'eigvals', 'eigvalsh', 'pinv', 'slogdet', 'det',
14 'svd', '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_,
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 >>> import numpy as np
276 >>> a = np.eye(2*3*4)
277 >>> a.shape = (2*3, 4, 2, 3, 4)
278 >>> rng = np.random.default_rng()
279 >>> b = rng.normal(size=(2*3, 4))
280 >>> x = np.linalg.tensorsolve(a, b)
281 >>> x.shape
282 (2, 3, 4)
283 >>> np.allclose(np.tensordot(a, x, axes=3), b)
284 True
286 """
287 a, wrap = _makearray(a)
288 b = asarray(b)
289 an = a.ndim
291 if axes is not None:
292 allaxes = list(range(0, an))
293 for k in axes:
294 allaxes.remove(k)
295 allaxes.insert(an, k)
296 a = a.transpose(allaxes)
298 oldshape = a.shape[-(an-b.ndim):]
299 prod = 1
300 for k in oldshape:
301 prod *= k
303 if a.size != prod ** 2:
304 raise LinAlgError(
305 "Input arrays must satisfy the requirement \
306 prod(a.shape[b.ndim:]) == prod(a.shape[:b.ndim])"
307 )
309 a = a.reshape(prod, prod)
310 b = b.ravel()
311 res = wrap(solve(a, b))
312 res.shape = oldshape
313 return res
316def _solve_dispatcher(a, b):
317 return (a, b)
320@array_function_dispatch(_solve_dispatcher)
321def solve(a, b):
322 """
323 Solve a linear matrix equation, or system of linear scalar equations.
325 Computes the "exact" solution, `x`, of the well-determined, i.e., full
326 rank, linear matrix equation `ax = b`.
328 Parameters
329 ----------
330 a : (..., M, M) array_like
331 Coefficient matrix.
332 b : {(M,), (..., M, K)}, array_like
333 Ordinate or "dependent variable" values.
335 Returns
336 -------
337 x : {(..., M,), (..., M, K)} ndarray
338 Solution to the system a x = b. Returned shape is (..., M) if b is
339 shape (M,) and (..., M, K) if b is (..., M, K), where the "..." part is
340 broadcasted between a and b.
342 Raises
343 ------
344 LinAlgError
345 If `a` is singular or not square.
347 See Also
348 --------
349 scipy.linalg.solve : Similar function in SciPy.
351 Notes
352 -----
353 Broadcasting rules apply, see the `numpy.linalg` documentation for
354 details.
356 The solutions are computed using LAPACK routine ``_gesv``.
358 `a` must be square and of full-rank, i.e., all rows (or, equivalently,
359 columns) must be linearly independent; if either is not true, use
360 `lstsq` for the least-squares best "solution" of the
361 system/equation.
363 .. versionchanged:: 2.0
365 The b array is only treated as a shape (M,) column vector if it is
366 exactly 1-dimensional. In all other instances it is treated as a stack
367 of (M, K) matrices. Previously b would be treated as a stack of (M,)
368 vectors if b.ndim was equal to a.ndim - 1.
370 References
371 ----------
372 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
373 FL, Academic Press, Inc., 1980, pg. 22.
375 Examples
376 --------
377 Solve the system of equations:
378 ``x0 + 2 * x1 = 1`` and
379 ``3 * x0 + 5 * x1 = 2``:
381 >>> import numpy as np
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 >>> import numpy as np
455 >>> a = np.eye(4*6)
456 >>> a.shape = (4, 6, 8, 3)
457 >>> ainv = np.linalg.tensorinv(a, ind=2)
458 >>> ainv.shape
459 (8, 3, 4, 6)
460 >>> rng = np.random.default_rng()
461 >>> b = rng.normal(size=(4, 6))
462 >>> np.allclose(np.tensordot(ainv, b), np.linalg.tensorsolve(a, b))
463 True
465 >>> a = np.eye(4*6)
466 >>> a.shape = (24, 8, 3)
467 >>> ainv = np.linalg.tensorinv(a, ind=1)
468 >>> ainv.shape
469 (8, 3, 24)
470 >>> rng = np.random.default_rng()
471 >>> b = rng.normal(size=24)
472 >>> np.allclose(np.tensordot(ainv, b, 1), np.linalg.tensorsolve(a, b))
473 True
475 """
476 a = asarray(a)
477 oldshape = a.shape
478 prod = 1
479 if ind > 0:
480 invshape = oldshape[ind:] + oldshape[:ind]
481 for k in oldshape[ind:]:
482 prod *= k
483 else:
484 raise ValueError("Invalid ind argument.")
485 a = a.reshape(prod, -1)
486 ia = inv(a)
487 return ia.reshape(*invshape)
490# Matrix inversion
492def _unary_dispatcher(a):
493 return (a,)
496@array_function_dispatch(_unary_dispatcher)
497def inv(a):
498 """
499 Compute the inverse of a matrix.
501 Given a square matrix `a`, return the matrix `ainv` satisfying
502 ``a @ ainv = ainv @ a = eye(a.shape[0])``.
504 Parameters
505 ----------
506 a : (..., M, M) array_like
507 Matrix to be inverted.
509 Returns
510 -------
511 ainv : (..., M, M) ndarray or matrix
512 Inverse of the matrix `a`.
514 Raises
515 ------
516 LinAlgError
517 If `a` is not square or inversion fails.
519 See Also
520 --------
521 scipy.linalg.inv : Similar function in SciPy.
522 numpy.linalg.cond : Compute the condition number of a matrix.
523 numpy.linalg.svd : Compute the singular value decomposition of a matrix.
525 Notes
526 -----
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 >>> import numpy as np
542 >>> from numpy.linalg import inv
543 >>> a = np.array([[1., 2.], [3., 4.]])
544 >>> ainv = inv(a)
545 >>> np.allclose(a @ ainv, np.eye(2))
546 True
547 >>> np.allclose(ainv @ a, np.eye(2))
548 True
550 If a is a matrix object, then the return value is a matrix as well:
552 >>> ainv = inv(np.matrix(a))
553 >>> ainv
554 matrix([[-2. , 1. ],
555 [ 1.5, -0.5]])
557 Inverses of several matrices can be computed at once:
559 >>> a = np.array([[[1., 2.], [3., 4.]], [[1, 3], [3, 5]]])
560 >>> inv(a)
561 array([[[-2. , 1. ],
562 [ 1.5 , -0.5 ]],
563 [[-1.25, 0.75],
564 [ 0.75, -0.25]]])
566 If a matrix is close to singular, the computed inverse may not satisfy
567 ``a @ ainv = ainv @ a = eye(a.shape[0])`` even if a `LinAlgError`
568 is not raised:
570 >>> a = np.array([[2,4,6],[2,0,2],[6,8,14]])
571 >>> inv(a) # No errors raised
572 array([[-1.12589991e+15, -5.62949953e+14, 5.62949953e+14],
573 [-1.12589991e+15, -5.62949953e+14, 5.62949953e+14],
574 [ 1.12589991e+15, 5.62949953e+14, -5.62949953e+14]])
575 >>> a @ inv(a)
576 array([[ 0. , -0.5 , 0. ], # may vary
577 [-0.5 , 0.625, 0.25 ],
578 [ 0. , 0. , 1. ]])
580 To detect ill-conditioned matrices, you can use `numpy.linalg.cond` to
581 compute its *condition number* [1]_. The larger the condition number, the
582 more ill-conditioned the matrix is. As a rule of thumb, if the condition
583 number ``cond(a) = 10**k``, then you may lose up to ``k`` digits of
584 accuracy on top of what would be lost to the numerical method due to loss
585 of precision from arithmetic methods.
587 >>> from numpy.linalg import cond
588 >>> cond(a)
589 np.float64(8.659885634118668e+17) # may vary
591 It is also possible to detect ill-conditioning by inspecting the matrix's
592 singular values directly. The ratio between the largest and the smallest
593 singular value is the condition number:
595 >>> from numpy.linalg import svd
596 >>> sigma = svd(a, compute_uv=False) # Do not compute singular vectors
597 >>> sigma.max()/sigma.min()
598 8.659885634118668e+17 # may vary
600 """
601 a, wrap = _makearray(a)
602 _assert_stacked_2d(a)
603 _assert_stacked_square(a)
604 t, result_t = _commonType(a)
606 signature = 'D->D' if isComplexType(t) else 'd->d'
607 with errstate(call=_raise_linalgerror_singular, invalid='call',
608 over='ignore', divide='ignore', under='ignore'):
609 ainv = _umath_linalg.inv(a, signature=signature)
610 return wrap(ainv.astype(result_t, copy=False))
613def _matrix_power_dispatcher(a, n):
614 return (a,)
617@array_function_dispatch(_matrix_power_dispatcher)
618def matrix_power(a, n):
619 """
620 Raise a square matrix to the (integer) power `n`.
622 For positive integers `n`, the power is computed by repeated matrix
623 squarings and matrix multiplications. If ``n == 0``, the identity matrix
624 of the same shape as M is returned. If ``n < 0``, the inverse
625 is computed and then raised to the ``abs(n)``.
627 .. note:: Stacks of object matrices are not currently supported.
629 Parameters
630 ----------
631 a : (..., M, M) array_like
632 Matrix to be "powered".
633 n : int
634 The exponent can be any integer or long integer, positive,
635 negative, or zero.
637 Returns
638 -------
639 a**n : (..., M, M) ndarray or matrix object
640 The return value is the same shape and type as `M`;
641 if the exponent is positive or zero then the type of the
642 elements is the same as those of `M`. If the exponent is
643 negative the elements are floating-point.
645 Raises
646 ------
647 LinAlgError
648 For matrices that are not square or that (for negative powers) cannot
649 be inverted numerically.
651 Examples
652 --------
653 >>> import numpy as np
654 >>> from numpy.linalg import matrix_power
655 >>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit
656 >>> matrix_power(i, 3) # should = -i
657 array([[ 0, -1],
658 [ 1, 0]])
659 >>> matrix_power(i, 0)
660 array([[1, 0],
661 [0, 1]])
662 >>> matrix_power(i, -3) # should = 1/(-i) = i, but w/ f.p. elements
663 array([[ 0., 1.],
664 [-1., 0.]])
666 Somewhat more sophisticated example
668 >>> q = np.zeros((4, 4))
669 >>> q[0:2, 0:2] = -i
670 >>> q[2:4, 2:4] = i
671 >>> q # one of the three quaternion units not equal to 1
672 array([[ 0., -1., 0., 0.],
673 [ 1., 0., 0., 0.],
674 [ 0., 0., 0., 1.],
675 [ 0., 0., -1., 0.]])
676 >>> matrix_power(q, 2) # = -np.eye(4)
677 array([[-1., 0., 0., 0.],
678 [ 0., -1., 0., 0.],
679 [ 0., 0., -1., 0.],
680 [ 0., 0., 0., -1.]])
682 """
683 a = asanyarray(a)
684 _assert_stacked_2d(a)
685 _assert_stacked_square(a)
687 try:
688 n = operator.index(n)
689 except TypeError as e:
690 raise TypeError("exponent must be an integer") from e
692 # Fall back on dot for object arrays. Object arrays are not supported by
693 # the current implementation of matmul using einsum
694 if a.dtype != object:
695 fmatmul = matmul
696 elif a.ndim == 2:
697 fmatmul = dot
698 else:
699 raise NotImplementedError(
700 "matrix_power not supported for stacks of object arrays")
702 if n == 0:
703 a = empty_like(a)
704 a[...] = eye(a.shape[-2], dtype=a.dtype)
705 return a
707 elif n < 0:
708 a = inv(a)
709 n = abs(n)
711 # short-cuts.
712 if n == 1:
713 return a
715 elif n == 2:
716 return fmatmul(a, a)
718 elif n == 3:
719 return fmatmul(fmatmul(a, a), a)
721 # Use binary decomposition to reduce the number of matrix multiplications.
722 # Here, we iterate over the bits of n, from LSB to MSB, raise `a` to
723 # increasing powers of 2, and multiply into the result as needed.
724 z = result = None
725 while n > 0:
726 z = a if z is None else fmatmul(z, z)
727 n, bit = divmod(n, 2)
728 if bit:
729 result = z if result is None else fmatmul(result, z)
731 return result
734# Cholesky decomposition
736def _cholesky_dispatcher(a, /, *, upper=None):
737 return (a,)
740@array_function_dispatch(_cholesky_dispatcher)
741def cholesky(a, /, *, upper=False):
742 """
743 Cholesky decomposition.
745 Return the lower or upper Cholesky decomposition, ``L * L.H`` or
746 ``U.H * U``, of the square matrix ``a``, where ``L`` is lower-triangular,
747 ``U`` is upper-triangular, and ``.H`` is the conjugate transpose operator
748 (which is the ordinary transpose if ``a`` is real-valued). ``a`` must be
749 Hermitian (symmetric if real-valued) and positive-definite. No checking is
750 performed to verify whether ``a`` is Hermitian or not. In addition, only
751 the lower or upper-triangular and diagonal elements of ``a`` are used.
752 Only ``L`` or ``U`` is actually returned.
754 Parameters
755 ----------
756 a : (..., M, M) array_like
757 Hermitian (symmetric if all elements are real), positive-definite
758 input matrix.
759 upper : bool
760 If ``True``, the result must be the upper-triangular Cholesky factor.
761 If ``False``, the result must be the lower-triangular Cholesky factor.
762 Default: ``False``.
764 Returns
765 -------
766 L : (..., M, M) array_like
767 Lower or upper-triangular Cholesky factor of `a`. Returns a matrix
768 object if `a` is a matrix object.
770 Raises
771 ------
772 LinAlgError
773 If the decomposition fails, for example, if `a` is not
774 positive-definite.
776 See Also
777 --------
778 scipy.linalg.cholesky : Similar function in SciPy.
779 scipy.linalg.cholesky_banded : Cholesky decompose a banded Hermitian
780 positive-definite matrix.
781 scipy.linalg.cho_factor : Cholesky decomposition of a matrix, to use in
782 `scipy.linalg.cho_solve`.
784 Notes
785 -----
786 Broadcasting rules apply, see the `numpy.linalg` documentation for
787 details.
789 The Cholesky decomposition is often used as a fast way of solving
791 .. math:: A \\mathbf{x} = \\mathbf{b}
793 (when `A` is both Hermitian/symmetric and positive-definite).
795 First, we solve for :math:`\\mathbf{y}` in
797 .. math:: L \\mathbf{y} = \\mathbf{b},
799 and then for :math:`\\mathbf{x}` in
801 .. math:: L^{H} \\mathbf{x} = \\mathbf{y}.
803 Examples
804 --------
805 >>> import numpy as np
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 Examples
877 --------
878 Make a (*very* coarse) grid for computing a Mandelbrot set:
880 >>> rl = np.linalg.outer(np.ones((5,)), np.linspace(-2, 2, 5))
881 >>> rl
882 array([[-2., -1., 0., 1., 2.],
883 [-2., -1., 0., 1., 2.],
884 [-2., -1., 0., 1., 2.],
885 [-2., -1., 0., 1., 2.],
886 [-2., -1., 0., 1., 2.]])
887 >>> im = np.linalg.outer(1j*np.linspace(2, -2, 5), np.ones((5,)))
888 >>> im
889 array([[0.+2.j, 0.+2.j, 0.+2.j, 0.+2.j, 0.+2.j],
890 [0.+1.j, 0.+1.j, 0.+1.j, 0.+1.j, 0.+1.j],
891 [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
892 [0.-1.j, 0.-1.j, 0.-1.j, 0.-1.j, 0.-1.j],
893 [0.-2.j, 0.-2.j, 0.-2.j, 0.-2.j, 0.-2.j]])
894 >>> grid = rl + im
895 >>> grid
896 array([[-2.+2.j, -1.+2.j, 0.+2.j, 1.+2.j, 2.+2.j],
897 [-2.+1.j, -1.+1.j, 0.+1.j, 1.+1.j, 2.+1.j],
898 [-2.+0.j, -1.+0.j, 0.+0.j, 1.+0.j, 2.+0.j],
899 [-2.-1.j, -1.-1.j, 0.-1.j, 1.-1.j, 2.-1.j],
900 [-2.-2.j, -1.-2.j, 0.-2.j, 1.-2.j, 2.-2.j]])
902 An example using a "vector" of letters:
904 >>> x = np.array(['a', 'b', 'c'], dtype=object)
905 >>> np.linalg.outer(x, [1, 2, 3])
906 array([['a', 'aa', 'aaa'],
907 ['b', 'bb', 'bbb'],
908 ['c', 'cc', 'ccc']], dtype=object)
910 """
911 x1 = asanyarray(x1)
912 x2 = asanyarray(x2)
913 if x1.ndim != 1 or x2.ndim != 1:
914 raise ValueError(
915 "Input arrays must be one-dimensional, but they are "
916 f"{x1.ndim=} and {x2.ndim=}."
917 )
918 return _core_outer(x1, x2, out=None)
921# QR decomposition
924def _qr_dispatcher(a, mode=None):
925 return (a,)
928@array_function_dispatch(_qr_dispatcher)
929def qr(a, mode='reduced'):
930 """
931 Compute the qr factorization of a matrix.
933 Factor the matrix `a` as *qr*, where `q` is orthonormal and `r` is
934 upper-triangular.
936 Parameters
937 ----------
938 a : array_like, shape (..., M, N)
939 An array-like object with the dimensionality of at least 2.
940 mode : {'reduced', 'complete', 'r', 'raw'}, optional, default: 'reduced'
941 If K = min(M, N), then
943 * 'reduced' : returns Q, R with dimensions (..., M, K), (..., K, N)
944 * 'complete' : returns Q, R with dimensions (..., M, M), (..., M, N)
945 * 'r' : returns R only with dimensions (..., K, N)
946 * 'raw' : returns h, tau with dimensions (..., N, M), (..., K,)
948 The options 'reduced', 'complete, and 'raw' are new in numpy 1.8,
949 see the notes for more information. The default is 'reduced', and to
950 maintain backward compatibility with earlier versions of numpy both
951 it and the old default 'full' can be omitted. Note that array h
952 returned in 'raw' mode is transposed for calling Fortran. The
953 'economic' mode is deprecated. The modes 'full' and 'economic' may
954 be passed using only the first letter for backwards compatibility,
955 but all others must be spelled out. See the Notes for more
956 explanation.
959 Returns
960 -------
961 When mode is 'reduced' or 'complete', the result will be a namedtuple with
962 the attributes `Q` and `R`.
964 Q : ndarray of float or complex, optional
965 A matrix with orthonormal columns. When mode = 'complete' the
966 result is an orthogonal/unitary matrix depending on whether or not
967 a is real/complex. The determinant may be either +/- 1 in that
968 case. In case the number of dimensions in the input array is
969 greater than 2 then a stack of the matrices with above properties
970 is returned.
971 R : ndarray of float or complex, optional
972 The upper-triangular matrix or a stack of upper-triangular
973 matrices if the number of dimensions in the input array is greater
974 than 2.
975 (h, tau) : ndarrays of np.double or np.cdouble, optional
976 The array h contains the Householder reflectors that generate q
977 along with r. The tau array contains scaling factors for the
978 reflectors. In the deprecated 'economic' mode only h is returned.
980 Raises
981 ------
982 LinAlgError
983 If factoring fails.
985 See Also
986 --------
987 scipy.linalg.qr : Similar function in SciPy.
988 scipy.linalg.rq : Compute RQ decomposition of a matrix.
990 Notes
991 -----
992 This is an interface to the LAPACK routines ``dgeqrf``, ``zgeqrf``,
993 ``dorgqr``, and ``zungqr``.
995 For more information on the qr factorization, see for example:
996 https://en.wikipedia.org/wiki/QR_factorization
998 Subclasses of `ndarray` are preserved except for the 'raw' mode. So if
999 `a` is of type `matrix`, all the return values will be matrices too.
1001 New 'reduced', 'complete', and 'raw' options for mode were added in
1002 NumPy 1.8.0 and the old option 'full' was made an alias of 'reduced'. In
1003 addition the options 'full' and 'economic' were deprecated. Because
1004 'full' was the previous default and 'reduced' is the new default,
1005 backward compatibility can be maintained by letting `mode` default.
1006 The 'raw' option was added so that LAPACK routines that can multiply
1007 arrays by q using the Householder reflectors can be used. Note that in
1008 this case the returned arrays are of type np.double or np.cdouble and
1009 the h array is transposed to be FORTRAN compatible. No routines using
1010 the 'raw' return are currently exposed by numpy, but some are available
1011 in lapack_lite and just await the necessary work.
1013 Examples
1014 --------
1015 >>> import numpy as np
1016 >>> rng = np.random.default_rng()
1017 >>> a = rng.normal(size=(9, 6))
1018 >>> Q, R = np.linalg.qr(a)
1019 >>> np.allclose(a, np.dot(Q, R)) # a does equal QR
1020 True
1021 >>> R2 = np.linalg.qr(a, mode='r')
1022 >>> np.allclose(R, R2) # mode='r' returns the same R as mode='full'
1023 True
1024 >>> a = np.random.normal(size=(3, 2, 2)) # Stack of 2 x 2 matrices as input
1025 >>> Q, R = np.linalg.qr(a)
1026 >>> Q.shape
1027 (3, 2, 2)
1028 >>> R.shape
1029 (3, 2, 2)
1030 >>> np.allclose(a, np.matmul(Q, R))
1031 True
1033 Example illustrating a common use of `qr`: solving of least squares
1034 problems
1036 What are the least-squares-best `m` and `y0` in ``y = y0 + mx`` for
1037 the following data: {(0,1), (1,0), (1,2), (2,1)}. (Graph the points
1038 and you'll see that it should be y0 = 0, m = 1.) The answer is provided
1039 by solving the over-determined matrix equation ``Ax = b``, where::
1041 A = array([[0, 1], [1, 1], [1, 1], [2, 1]])
1042 x = array([[y0], [m]])
1043 b = array([[1], [0], [2], [1]])
1045 If A = QR such that Q is orthonormal (which is always possible via
1046 Gram-Schmidt), then ``x = inv(R) * (Q.T) * b``. (In numpy practice,
1047 however, we simply use `lstsq`.)
1049 >>> A = np.array([[0, 1], [1, 1], [1, 1], [2, 1]])
1050 >>> A
1051 array([[0, 1],
1052 [1, 1],
1053 [1, 1],
1054 [2, 1]])
1055 >>> b = np.array([1, 2, 2, 3])
1056 >>> Q, R = np.linalg.qr(A)
1057 >>> p = np.dot(Q.T, b)
1058 >>> np.dot(np.linalg.inv(R), p)
1059 array([ 1., 1.])
1061 """
1062 if mode not in ('reduced', 'complete', 'r', 'raw'):
1063 if mode in ('f', 'full'):
1064 # 2013-04-01, 1.8
1065 msg = (
1066 "The 'full' option is deprecated in favor of 'reduced'.\n"
1067 "For backward compatibility let mode default."
1068 )
1069 warnings.warn(msg, DeprecationWarning, stacklevel=2)
1070 mode = 'reduced'
1071 elif mode in ('e', 'economic'):
1072 # 2013-04-01, 1.8
1073 msg = "The 'economic' option is deprecated."
1074 warnings.warn(msg, DeprecationWarning, stacklevel=2)
1075 mode = 'economic'
1076 else:
1077 raise ValueError(f"Unrecognized mode '{mode}'")
1079 a, wrap = _makearray(a)
1080 _assert_stacked_2d(a)
1081 m, n = a.shape[-2:]
1082 t, result_t = _commonType(a)
1083 a = a.astype(t, copy=True)
1084 a = _to_native_byte_order(a)
1085 mn = min(m, n)
1087 signature = 'D->D' if isComplexType(t) else 'd->d'
1088 with errstate(call=_raise_linalgerror_qr, invalid='call',
1089 over='ignore', divide='ignore', under='ignore'):
1090 tau = _umath_linalg.qr_r_raw(a, signature=signature)
1092 # handle modes that don't return q
1093 if mode == 'r':
1094 r = triu(a[..., :mn, :])
1095 r = r.astype(result_t, copy=False)
1096 return wrap(r)
1098 if mode == 'raw':
1099 q = transpose(a)
1100 q = q.astype(result_t, copy=False)
1101 tau = tau.astype(result_t, copy=False)
1102 return wrap(q), tau
1104 if mode == 'economic':
1105 a = a.astype(result_t, copy=False)
1106 return wrap(a)
1108 # mc is the number of columns in the resulting q
1109 # matrix. If the mode is complete then it is
1110 # same as number of rows, and if the mode is reduced,
1111 # then it is the minimum of number of rows and columns.
1112 if mode == 'complete' and m > n:
1113 mc = m
1114 gufunc = _umath_linalg.qr_complete
1115 else:
1116 mc = mn
1117 gufunc = _umath_linalg.qr_reduced
1119 signature = 'DD->D' if isComplexType(t) else 'dd->d'
1120 with errstate(call=_raise_linalgerror_qr, invalid='call',
1121 over='ignore', divide='ignore', under='ignore'):
1122 q = gufunc(a, tau, signature=signature)
1123 r = triu(a[..., :mc, :])
1125 q = q.astype(result_t, copy=False)
1126 r = r.astype(result_t, copy=False)
1128 return QRResult(wrap(q), wrap(r))
1130# Eigenvalues
1133@array_function_dispatch(_unary_dispatcher)
1134def eigvals(a):
1135 """
1136 Compute the eigenvalues of a general matrix.
1138 Main difference between `eigvals` and `eig`: the eigenvectors aren't
1139 returned.
1141 Parameters
1142 ----------
1143 a : (..., M, M) array_like
1144 A complex- or real-valued matrix whose eigenvalues will be computed.
1146 Returns
1147 -------
1148 w : (..., M,) ndarray
1149 The eigenvalues, each repeated according to its multiplicity.
1150 They are not necessarily ordered, nor are they necessarily
1151 real for real matrices.
1153 Raises
1154 ------
1155 LinAlgError
1156 If the eigenvalue computation does not converge.
1158 See Also
1159 --------
1160 eig : eigenvalues and right eigenvectors of general arrays
1161 eigvalsh : eigenvalues of real symmetric or complex Hermitian
1162 (conjugate symmetric) arrays.
1163 eigh : eigenvalues and eigenvectors of real symmetric or complex
1164 Hermitian (conjugate symmetric) arrays.
1165 scipy.linalg.eigvals : Similar function in SciPy.
1167 Notes
1168 -----
1169 Broadcasting rules apply, see the `numpy.linalg` documentation for
1170 details.
1172 This is implemented using the ``_geev`` LAPACK routines which compute
1173 the eigenvalues and eigenvectors of general square arrays.
1175 Examples
1176 --------
1177 Illustration, using the fact that the eigenvalues of a diagonal matrix
1178 are its diagonal elements, that multiplying a matrix on the left
1179 by an orthogonal matrix, `Q`, and on the right by `Q.T` (the transpose
1180 of `Q`), preserves the eigenvalues of the "middle" matrix. In other words,
1181 if `Q` is orthogonal, then ``Q * A * Q.T`` has the same eigenvalues as
1182 ``A``:
1184 >>> import numpy as np
1185 >>> from numpy import linalg as LA
1186 >>> x = np.random.random()
1187 >>> Q = np.array([[np.cos(x), -np.sin(x)], [np.sin(x), np.cos(x)]])
1188 >>> LA.norm(Q[0, :]), LA.norm(Q[1, :]), np.dot(Q[0, :],Q[1, :])
1189 (1.0, 1.0, 0.0)
1191 Now multiply a diagonal matrix by ``Q`` on one side and
1192 by ``Q.T`` on the other:
1194 >>> D = np.diag((-1,1))
1195 >>> LA.eigvals(D)
1196 array([-1., 1.])
1197 >>> A = np.dot(Q, D)
1198 >>> A = np.dot(A, Q.T)
1199 >>> LA.eigvals(A)
1200 array([ 1., -1.]) # random
1202 """
1203 a, wrap = _makearray(a)
1204 _assert_stacked_2d(a)
1205 _assert_stacked_square(a)
1206 _assert_finite(a)
1207 t, result_t = _commonType(a)
1209 signature = 'D->D' if isComplexType(t) else 'd->D'
1210 with errstate(call=_raise_linalgerror_eigenvalues_nonconvergence,
1211 invalid='call', over='ignore', divide='ignore',
1212 under='ignore'):
1213 w = _umath_linalg.eigvals(a, signature=signature)
1215 if not isComplexType(t):
1216 if all(w.imag == 0):
1217 w = w.real
1218 result_t = _realType(result_t)
1219 else:
1220 result_t = _complexType(result_t)
1222 return w.astype(result_t, copy=False)
1225def _eigvalsh_dispatcher(a, UPLO=None):
1226 return (a,)
1229@array_function_dispatch(_eigvalsh_dispatcher)
1230def eigvalsh(a, UPLO='L'):
1231 """
1232 Compute the eigenvalues of a complex Hermitian or real symmetric matrix.
1234 Main difference from eigh: the eigenvectors are not computed.
1236 Parameters
1237 ----------
1238 a : (..., M, M) array_like
1239 A complex- or real-valued matrix whose eigenvalues are to be
1240 computed.
1241 UPLO : {'L', 'U'}, optional
1242 Specifies whether the calculation is done with the lower triangular
1243 part of `a` ('L', default) or the upper triangular part ('U').
1244 Irrespective of this value only the real parts of the diagonal will
1245 be considered in the computation to preserve the notion of a Hermitian
1246 matrix. It therefore follows that the imaginary part of the diagonal
1247 will always be treated as zero.
1249 Returns
1250 -------
1251 w : (..., M,) ndarray
1252 The eigenvalues in ascending order, each repeated according to
1253 its multiplicity.
1255 Raises
1256 ------
1257 LinAlgError
1258 If the eigenvalue computation does not converge.
1260 See Also
1261 --------
1262 eigh : eigenvalues and eigenvectors of real symmetric or complex Hermitian
1263 (conjugate symmetric) arrays.
1264 eigvals : eigenvalues of general real or complex arrays.
1265 eig : eigenvalues and right eigenvectors of general real or complex
1266 arrays.
1267 scipy.linalg.eigvalsh : Similar function in SciPy.
1269 Notes
1270 -----
1271 Broadcasting rules apply, see the `numpy.linalg` documentation for
1272 details.
1274 The eigenvalues are computed using LAPACK routines ``_syevd``, ``_heevd``.
1276 Examples
1277 --------
1278 >>> import numpy as np
1279 >>> from numpy import linalg as LA
1280 >>> a = np.array([[1, -2j], [2j, 5]])
1281 >>> LA.eigvalsh(a)
1282 array([ 0.17157288, 5.82842712]) # may vary
1284 >>> # demonstrate the treatment of the imaginary part of the diagonal
1285 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1286 >>> a
1287 array([[5.+2.j, 9.-2.j],
1288 [0.+2.j, 2.-1.j]])
1289 >>> # with UPLO='L' this is numerically equivalent to using LA.eigvals()
1290 >>> # with:
1291 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1292 >>> b
1293 array([[5.+0.j, 0.-2.j],
1294 [0.+2.j, 2.+0.j]])
1295 >>> wa = LA.eigvalsh(a)
1296 >>> wb = LA.eigvals(b)
1297 >>> wa; wb
1298 array([1., 6.])
1299 array([6.+0.j, 1.+0.j])
1301 """
1302 UPLO = UPLO.upper()
1303 if UPLO not in ('L', 'U'):
1304 raise ValueError("UPLO argument must be 'L' or 'U'")
1306 if UPLO == 'L':
1307 gufunc = _umath_linalg.eigvalsh_lo
1308 else:
1309 gufunc = _umath_linalg.eigvalsh_up
1311 a, wrap = _makearray(a)
1312 _assert_stacked_2d(a)
1313 _assert_stacked_square(a)
1314 t, result_t = _commonType(a)
1315 signature = 'D->d' if isComplexType(t) else 'd->d'
1316 with errstate(call=_raise_linalgerror_eigenvalues_nonconvergence,
1317 invalid='call', over='ignore', divide='ignore',
1318 under='ignore'):
1319 w = gufunc(a, signature=signature)
1320 return w.astype(_realType(result_t), copy=False)
1322def _convertarray(a):
1323 t, result_t = _commonType(a)
1324 a = a.astype(t).T.copy()
1325 return a, t, result_t
1328# Eigenvectors
1331@array_function_dispatch(_unary_dispatcher)
1332def eig(a):
1333 """
1334 Compute the eigenvalues and right eigenvectors of a square array.
1336 Parameters
1337 ----------
1338 a : (..., M, M) array
1339 Matrices for which the eigenvalues and right eigenvectors will
1340 be computed
1342 Returns
1343 -------
1344 A namedtuple with the following attributes:
1346 eigenvalues : (..., M) array
1347 The eigenvalues, each repeated according to its multiplicity.
1348 The eigenvalues are not necessarily ordered. The resulting
1349 array will be of complex type, unless the imaginary part is
1350 zero in which case it will be cast to a real type. When `a`
1351 is real the resulting eigenvalues will be real (0 imaginary
1352 part) or occur in conjugate pairs
1354 eigenvectors : (..., M, M) array
1355 The normalized (unit "length") eigenvectors, such that the
1356 column ``eigenvectors[:,i]`` is the eigenvector corresponding to the
1357 eigenvalue ``eigenvalues[i]``.
1359 Raises
1360 ------
1361 LinAlgError
1362 If the eigenvalue computation does not converge.
1364 See Also
1365 --------
1366 eigvals : eigenvalues of a non-symmetric array.
1367 eigh : eigenvalues and eigenvectors of a real symmetric or complex
1368 Hermitian (conjugate symmetric) array.
1369 eigvalsh : eigenvalues of a real symmetric or complex Hermitian
1370 (conjugate symmetric) array.
1371 scipy.linalg.eig : Similar function in SciPy that also solves the
1372 generalized eigenvalue problem.
1373 scipy.linalg.schur : Best choice for unitary and other non-Hermitian
1374 normal matrices.
1376 Notes
1377 -----
1378 Broadcasting rules apply, see the `numpy.linalg` documentation for
1379 details.
1381 This is implemented using the ``_geev`` LAPACK routines which compute
1382 the eigenvalues and eigenvectors of general square arrays.
1384 The number `w` is an eigenvalue of `a` if there exists a vector `v` such
1385 that ``a @ v = w * v``. Thus, the arrays `a`, `eigenvalues`, and
1386 `eigenvectors` satisfy the equations ``a @ eigenvectors[:,i] =
1387 eigenvalues[i] * eigenvectors[:,i]`` for :math:`i \\in \\{0,...,M-1\\}`.
1389 The array `eigenvectors` may not be of maximum rank, that is, some of the
1390 columns may be linearly dependent, although round-off error may obscure
1391 that fact. If the eigenvalues are all different, then theoretically the
1392 eigenvectors are linearly independent and `a` can be diagonalized by a
1393 similarity transformation using `eigenvectors`, i.e, ``inv(eigenvectors) @
1394 a @ eigenvectors`` is diagonal.
1396 For non-Hermitian normal matrices the SciPy function `scipy.linalg.schur`
1397 is preferred because the matrix `eigenvectors` is guaranteed to be
1398 unitary, which is not the case when using `eig`. The Schur factorization
1399 produces an upper triangular matrix rather than a diagonal matrix, but for
1400 normal matrices only the diagonal of the upper triangular matrix is
1401 needed, the rest is roundoff error.
1403 Finally, it is emphasized that `eigenvectors` consists of the *right* (as
1404 in right-hand side) eigenvectors of `a`. A vector `y` satisfying ``y.T @ a
1405 = z * y.T`` for some number `z` is called a *left* eigenvector of `a`,
1406 and, in general, the left and right eigenvectors of a matrix are not
1407 necessarily the (perhaps conjugate) transposes of each other.
1409 References
1410 ----------
1411 G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, FL,
1412 Academic Press, Inc., 1980, Various pp.
1414 Examples
1415 --------
1416 >>> import numpy as np
1417 >>> from numpy import linalg as LA
1419 (Almost) trivial example with real eigenvalues and eigenvectors.
1421 >>> eigenvalues, eigenvectors = LA.eig(np.diag((1, 2, 3)))
1422 >>> eigenvalues
1423 array([1., 2., 3.])
1424 >>> eigenvectors
1425 array([[1., 0., 0.],
1426 [0., 1., 0.],
1427 [0., 0., 1.]])
1429 Real matrix possessing complex eigenvalues and eigenvectors;
1430 note that the eigenvalues are complex conjugates of each other.
1432 >>> eigenvalues, eigenvectors = LA.eig(np.array([[1, -1], [1, 1]]))
1433 >>> eigenvalues
1434 array([1.+1.j, 1.-1.j])
1435 >>> eigenvectors
1436 array([[0.70710678+0.j , 0.70710678-0.j ],
1437 [0. -0.70710678j, 0. +0.70710678j]])
1439 Complex-valued matrix with real eigenvalues (but complex-valued
1440 eigenvectors); note that ``a.conj().T == a``, i.e., `a` is Hermitian.
1442 >>> a = np.array([[1, 1j], [-1j, 1]])
1443 >>> eigenvalues, eigenvectors = LA.eig(a)
1444 >>> eigenvalues
1445 array([2.+0.j, 0.+0.j])
1446 >>> eigenvectors
1447 array([[ 0. +0.70710678j, 0.70710678+0.j ], # may vary
1448 [ 0.70710678+0.j , -0. +0.70710678j]])
1450 Be careful about round-off error!
1452 >>> a = np.array([[1 + 1e-9, 0], [0, 1 - 1e-9]])
1453 >>> # Theor. eigenvalues are 1 +/- 1e-9
1454 >>> eigenvalues, eigenvectors = LA.eig(a)
1455 >>> eigenvalues
1456 array([1., 1.])
1457 >>> eigenvectors
1458 array([[1., 0.],
1459 [0., 1.]])
1461 """
1462 a, wrap = _makearray(a)
1463 _assert_stacked_2d(a)
1464 _assert_stacked_square(a)
1465 _assert_finite(a)
1466 t, result_t = _commonType(a)
1468 signature = 'D->DD' if isComplexType(t) else 'd->DD'
1469 with errstate(call=_raise_linalgerror_eigenvalues_nonconvergence,
1470 invalid='call', over='ignore', divide='ignore',
1471 under='ignore'):
1472 w, vt = _umath_linalg.eig(a, signature=signature)
1474 if not isComplexType(t) and all(w.imag == 0.0):
1475 w = w.real
1476 vt = vt.real
1477 result_t = _realType(result_t)
1478 else:
1479 result_t = _complexType(result_t)
1481 vt = vt.astype(result_t, copy=False)
1482 return EigResult(w.astype(result_t, copy=False), wrap(vt))
1485@array_function_dispatch(_eigvalsh_dispatcher)
1486def eigh(a, UPLO='L'):
1487 """
1488 Return the eigenvalues and eigenvectors of a complex Hermitian
1489 (conjugate symmetric) or a real symmetric matrix.
1491 Returns two objects, a 1-D array containing the eigenvalues of `a`, and
1492 a 2-D square array or matrix (depending on the input type) of the
1493 corresponding eigenvectors (in columns).
1495 Parameters
1496 ----------
1497 a : (..., M, M) array
1498 Hermitian or real symmetric matrices whose eigenvalues and
1499 eigenvectors are to be computed.
1500 UPLO : {'L', 'U'}, optional
1501 Specifies whether the calculation is done with the lower triangular
1502 part of `a` ('L', default) or the upper triangular part ('U').
1503 Irrespective of this value only the real parts of the diagonal will
1504 be considered in the computation to preserve the notion of a Hermitian
1505 matrix. It therefore follows that the imaginary part of the diagonal
1506 will always be treated as zero.
1508 Returns
1509 -------
1510 A namedtuple with the following attributes:
1512 eigenvalues : (..., M) ndarray
1513 The eigenvalues in ascending order, each repeated according to
1514 its multiplicity.
1515 eigenvectors : {(..., M, M) ndarray, (..., M, M) matrix}
1516 The column ``eigenvectors[:, i]`` is the normalized eigenvector
1517 corresponding to the eigenvalue ``eigenvalues[i]``. Will return a
1518 matrix object if `a` is a matrix object.
1520 Raises
1521 ------
1522 LinAlgError
1523 If the eigenvalue computation does not converge.
1525 See Also
1526 --------
1527 eigvalsh : eigenvalues of real symmetric or complex Hermitian
1528 (conjugate symmetric) arrays.
1529 eig : eigenvalues and right eigenvectors for non-symmetric arrays.
1530 eigvals : eigenvalues of non-symmetric arrays.
1531 scipy.linalg.eigh : Similar function in SciPy (but also solves the
1532 generalized eigenvalue problem).
1534 Notes
1535 -----
1536 Broadcasting rules apply, see the `numpy.linalg` documentation for
1537 details.
1539 The eigenvalues/eigenvectors are computed using LAPACK routines ``_syevd``,
1540 ``_heevd``.
1542 The eigenvalues of real symmetric or complex Hermitian matrices are always
1543 real. [1]_ The array `eigenvalues` of (column) eigenvectors is unitary and
1544 `a`, `eigenvalues`, and `eigenvectors` satisfy the equations ``dot(a,
1545 eigenvectors[:, i]) = eigenvalues[i] * eigenvectors[:, i]``.
1547 References
1548 ----------
1549 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
1550 FL, Academic Press, Inc., 1980, pg. 222.
1552 Examples
1553 --------
1554 >>> import numpy as np
1555 >>> from numpy import linalg as LA
1556 >>> a = np.array([[1, -2j], [2j, 5]])
1557 >>> a
1558 array([[ 1.+0.j, -0.-2.j],
1559 [ 0.+2.j, 5.+0.j]])
1560 >>> eigenvalues, eigenvectors = LA.eigh(a)
1561 >>> eigenvalues
1562 array([0.17157288, 5.82842712])
1563 >>> eigenvectors
1564 array([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1565 [ 0. +0.38268343j, 0. -0.92387953j]])
1567 >>> (np.dot(a, eigenvectors[:, 0]) -
1568 ... eigenvalues[0] * eigenvectors[:, 0]) # verify 1st eigenval/vec pair
1569 array([5.55111512e-17+0.0000000e+00j, 0.00000000e+00+1.2490009e-16j])
1570 >>> (np.dot(a, eigenvectors[:, 1]) -
1571 ... eigenvalues[1] * eigenvectors[:, 1]) # verify 2nd eigenval/vec pair
1572 array([0.+0.j, 0.+0.j])
1574 >>> A = np.matrix(a) # what happens if input is a matrix object
1575 >>> A
1576 matrix([[ 1.+0.j, -0.-2.j],
1577 [ 0.+2.j, 5.+0.j]])
1578 >>> eigenvalues, eigenvectors = LA.eigh(A)
1579 >>> eigenvalues
1580 array([0.17157288, 5.82842712])
1581 >>> eigenvectors
1582 matrix([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1583 [ 0. +0.38268343j, 0. -0.92387953j]])
1585 >>> # demonstrate the treatment of the imaginary part of the diagonal
1586 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1587 >>> a
1588 array([[5.+2.j, 9.-2.j],
1589 [0.+2.j, 2.-1.j]])
1590 >>> # with UPLO='L' this is numerically equivalent to using LA.eig() with:
1591 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1592 >>> b
1593 array([[5.+0.j, 0.-2.j],
1594 [0.+2.j, 2.+0.j]])
1595 >>> wa, va = LA.eigh(a)
1596 >>> wb, vb = LA.eig(b)
1597 >>> wa
1598 array([1., 6.])
1599 >>> wb
1600 array([6.+0.j, 1.+0.j])
1601 >>> va
1602 array([[-0.4472136 +0.j , -0.89442719+0.j ], # may vary
1603 [ 0. +0.89442719j, 0. -0.4472136j ]])
1604 >>> vb
1605 array([[ 0.89442719+0.j , -0. +0.4472136j],
1606 [-0. +0.4472136j, 0.89442719+0.j ]])
1608 """
1609 UPLO = UPLO.upper()
1610 if UPLO not in ('L', 'U'):
1611 raise ValueError("UPLO argument must be 'L' or 'U'")
1613 a, wrap = _makearray(a)
1614 _assert_stacked_2d(a)
1615 _assert_stacked_square(a)
1616 t, result_t = _commonType(a)
1618 if UPLO == 'L':
1619 gufunc = _umath_linalg.eigh_lo
1620 else:
1621 gufunc = _umath_linalg.eigh_up
1623 signature = 'D->dD' if isComplexType(t) else 'd->dd'
1624 with errstate(call=_raise_linalgerror_eigenvalues_nonconvergence,
1625 invalid='call', over='ignore', divide='ignore',
1626 under='ignore'):
1627 w, vt = gufunc(a, signature=signature)
1628 w = w.astype(_realType(result_t), copy=False)
1629 vt = vt.astype(result_t, copy=False)
1630 return EighResult(w, wrap(vt))
1633# Singular value decomposition
1635def _svd_dispatcher(a, full_matrices=None, compute_uv=None, hermitian=None):
1636 return (a,)
1639@array_function_dispatch(_svd_dispatcher)
1640def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
1641 """
1642 Singular Value Decomposition.
1644 When `a` is a 2D array, and ``full_matrices=False``, then it is
1645 factorized as ``u @ np.diag(s) @ vh = (u * s) @ vh``, where
1646 `u` and the Hermitian transpose of `vh` are 2D arrays with
1647 orthonormal columns and `s` is a 1D array of `a`'s singular
1648 values. When `a` is higher-dimensional, SVD is applied in
1649 stacked mode as explained below.
1651 Parameters
1652 ----------
1653 a : (..., M, N) array_like
1654 A real or complex array with ``a.ndim >= 2``.
1655 full_matrices : bool, optional
1656 If True (default), `u` and `vh` have the shapes ``(..., M, M)`` and
1657 ``(..., N, N)``, respectively. Otherwise, the shapes are
1658 ``(..., M, K)`` and ``(..., K, N)``, respectively, where
1659 ``K = min(M, N)``.
1660 compute_uv : bool, optional
1661 Whether or not to compute `u` and `vh` in addition to `s`. True
1662 by default.
1663 hermitian : bool, optional
1664 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
1665 enabling a more efficient method for finding singular values.
1666 Defaults to False.
1668 Returns
1669 -------
1670 When `compute_uv` is True, the result is a namedtuple with the following
1671 attribute names:
1673 U : { (..., M, M), (..., M, K) } array
1674 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1675 size as those of the input `a`. The size of the last two dimensions
1676 depends on the value of `full_matrices`. Only returned when
1677 `compute_uv` is True.
1678 S : (..., K) array
1679 Vector(s) with the singular values, within each vector sorted in
1680 descending order. The first ``a.ndim - 2`` dimensions have the same
1681 size as those of the input `a`.
1682 Vh : { (..., N, N), (..., K, N) } array
1683 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1684 size as those of the input `a`. The size of the last two dimensions
1685 depends on the value of `full_matrices`. Only returned when
1686 `compute_uv` is True.
1688 Raises
1689 ------
1690 LinAlgError
1691 If SVD computation does not converge.
1693 See Also
1694 --------
1695 scipy.linalg.svd : Similar function in SciPy.
1696 scipy.linalg.svdvals : Compute singular values of a matrix.
1698 Notes
1699 -----
1700 The decomposition is performed using LAPACK routine ``_gesdd``.
1702 SVD is usually described for the factorization of a 2D matrix :math:`A`.
1703 The higher-dimensional case will be discussed below. In the 2D case, SVD is
1704 written as :math:`A = U S V^H`, where :math:`A = a`, :math:`U= u`,
1705 :math:`S= \\mathtt{np.diag}(s)` and :math:`V^H = vh`. The 1D array `s`
1706 contains the singular values of `a` and `u` and `vh` are unitary. The rows
1707 of `vh` are the eigenvectors of :math:`A^H A` and the columns of `u` are
1708 the eigenvectors of :math:`A A^H`. In both cases the corresponding
1709 (possibly non-zero) eigenvalues are given by ``s**2``.
1711 If `a` has more than two dimensions, then broadcasting rules apply, as
1712 explained in :ref:`routines.linalg-broadcasting`. This means that SVD is
1713 working in "stacked" mode: it iterates over all indices of the first
1714 ``a.ndim - 2`` dimensions and for each combination SVD is applied to the
1715 last two indices. The matrix `a` can be reconstructed from the
1716 decomposition with either ``(u * s[..., None, :]) @ vh`` or
1717 ``u @ (s[..., None] * vh)``. (The ``@`` operator can be replaced by the
1718 function ``np.matmul`` for python versions below 3.5.)
1720 If `a` is a ``matrix`` object (as opposed to an ``ndarray``), then so are
1721 all the return values.
1723 Examples
1724 --------
1725 >>> import numpy as np
1726 >>> rng = np.random.default_rng()
1727 >>> a = rng.normal(size=(9, 6)) + 1j*rng.normal(size=(9, 6))
1728 >>> b = rng.normal(size=(2, 7, 8, 3)) + 1j*rng.normal(size=(2, 7, 8, 3))
1731 Reconstruction based on full SVD, 2D case:
1733 >>> U, S, Vh = np.linalg.svd(a, full_matrices=True)
1734 >>> U.shape, S.shape, Vh.shape
1735 ((9, 9), (6,), (6, 6))
1736 >>> np.allclose(a, np.dot(U[:, :6] * S, Vh))
1737 True
1738 >>> smat = np.zeros((9, 6), dtype=complex)
1739 >>> smat[:6, :6] = np.diag(S)
1740 >>> np.allclose(a, np.dot(U, np.dot(smat, Vh)))
1741 True
1743 Reconstruction based on reduced SVD, 2D case:
1745 >>> U, S, Vh = np.linalg.svd(a, full_matrices=False)
1746 >>> U.shape, S.shape, Vh.shape
1747 ((9, 6), (6,), (6, 6))
1748 >>> np.allclose(a, np.dot(U * S, Vh))
1749 True
1750 >>> smat = np.diag(S)
1751 >>> np.allclose(a, np.dot(U, np.dot(smat, Vh)))
1752 True
1754 Reconstruction based on full SVD, 4D case:
1756 >>> U, S, Vh = np.linalg.svd(b, full_matrices=True)
1757 >>> U.shape, S.shape, Vh.shape
1758 ((2, 7, 8, 8), (2, 7, 3), (2, 7, 3, 3))
1759 >>> np.allclose(b, np.matmul(U[..., :3] * S[..., None, :], Vh))
1760 True
1761 >>> np.allclose(b, np.matmul(U[..., :3], S[..., None] * Vh))
1762 True
1764 Reconstruction based on reduced SVD, 4D case:
1766 >>> U, S, Vh = np.linalg.svd(b, full_matrices=False)
1767 >>> U.shape, S.shape, Vh.shape
1768 ((2, 7, 8, 3), (2, 7, 3), (2, 7, 3, 3))
1769 >>> np.allclose(b, np.matmul(U * S[..., None, :], Vh))
1770 True
1771 >>> np.allclose(b, np.matmul(U, S[..., None] * Vh))
1772 True
1774 """
1775 import numpy as _nx
1776 a, wrap = _makearray(a)
1778 if hermitian:
1779 # note: lapack svd returns eigenvalues with s ** 2 sorted descending,
1780 # but eig returns s sorted ascending, so we re-order the eigenvalues
1781 # and related arrays to have the correct order
1782 if compute_uv:
1783 s, u = eigh(a)
1784 sgn = sign(s)
1785 s = abs(s)
1786 sidx = argsort(s)[..., ::-1]
1787 sgn = _nx.take_along_axis(sgn, sidx, axis=-1)
1788 s = _nx.take_along_axis(s, sidx, axis=-1)
1789 u = _nx.take_along_axis(u, sidx[..., None, :], axis=-1)
1790 # singular values are unsigned, move the sign into v
1791 vt = transpose(u * sgn[..., None, :]).conjugate()
1792 return SVDResult(wrap(u), s, wrap(vt))
1793 else:
1794 s = eigvalsh(a)
1795 s = abs(s)
1796 return sort(s)[..., ::-1]
1798 _assert_stacked_2d(a)
1799 t, result_t = _commonType(a)
1801 m, n = a.shape[-2:]
1802 if compute_uv:
1803 if full_matrices:
1804 gufunc = _umath_linalg.svd_f
1805 else:
1806 gufunc = _umath_linalg.svd_s
1808 signature = 'D->DdD' if isComplexType(t) else 'd->ddd'
1809 with errstate(call=_raise_linalgerror_svd_nonconvergence,
1810 invalid='call', over='ignore', divide='ignore',
1811 under='ignore'):
1812 u, s, vh = gufunc(a, signature=signature)
1813 u = u.astype(result_t, copy=False)
1814 s = s.astype(_realType(result_t), copy=False)
1815 vh = vh.astype(result_t, copy=False)
1816 return SVDResult(wrap(u), s, wrap(vh))
1817 else:
1818 signature = 'D->d' if isComplexType(t) else 'd->d'
1819 with errstate(call=_raise_linalgerror_svd_nonconvergence,
1820 invalid='call', over='ignore', divide='ignore',
1821 under='ignore'):
1822 s = _umath_linalg.svd(a, signature=signature)
1823 s = s.astype(_realType(result_t), copy=False)
1824 return s
1827def _svdvals_dispatcher(x):
1828 return (x,)
1831@array_function_dispatch(_svdvals_dispatcher)
1832def svdvals(x, /):
1833 """
1834 Returns the singular values of a matrix (or a stack of matrices) ``x``.
1835 When x is a stack of matrices, the function will compute the singular
1836 values for each matrix in the stack.
1838 This function is Array API compatible.
1840 Calling ``np.svdvals(x)`` to get singular values is the same as
1841 ``np.svd(x, compute_uv=False, hermitian=False)``.
1843 Parameters
1844 ----------
1845 x : (..., M, N) array_like
1846 Input array having shape (..., M, N) and whose last two
1847 dimensions form matrices on which to perform singular value
1848 decomposition. Should have a floating-point data type.
1850 Returns
1851 -------
1852 out : ndarray
1853 An array with shape (..., K) that contains the vector(s)
1854 of singular values of length K, where K = min(M, N).
1856 See Also
1857 --------
1858 scipy.linalg.svdvals : Compute singular values of a matrix.
1860 Examples
1861 --------
1863 >>> np.linalg.svdvals([[1, 2, 3, 4, 5],
1864 ... [1, 4, 9, 16, 25],
1865 ... [1, 8, 27, 64, 125]])
1866 array([146.68862757, 5.57510612, 0.60393245])
1868 Determine the rank of a matrix using singular values:
1870 >>> s = np.linalg.svdvals([[1, 2, 3],
1871 ... [2, 4, 6],
1872 ... [-1, 1, -1]]); s
1873 array([8.38434191e+00, 1.64402274e+00, 2.31534378e-16])
1874 >>> np.count_nonzero(s > 1e-10) # Matrix of rank 2
1875 2
1877 """
1878 return svd(x, compute_uv=False, hermitian=False)
1881def _cond_dispatcher(x, p=None):
1882 return (x,)
1885@array_function_dispatch(_cond_dispatcher)
1886def cond(x, p=None):
1887 """
1888 Compute the condition number of a matrix.
1890 This function is capable of returning the condition number using
1891 one of seven different norms, depending on the value of `p` (see
1892 Parameters below).
1894 Parameters
1895 ----------
1896 x : (..., M, N) array_like
1897 The matrix whose condition number is sought.
1898 p : {None, 1, -1, 2, -2, inf, -inf, 'fro'}, optional
1899 Order of the norm used in the condition number computation:
1901 ===== ============================
1902 p norm for matrices
1903 ===== ============================
1904 None 2-norm, computed directly using the ``SVD``
1905 'fro' Frobenius norm
1906 inf max(sum(abs(x), axis=1))
1907 -inf min(sum(abs(x), axis=1))
1908 1 max(sum(abs(x), axis=0))
1909 -1 min(sum(abs(x), axis=0))
1910 2 2-norm (largest sing. value)
1911 -2 smallest singular value
1912 ===== ============================
1914 inf means the `numpy.inf` object, and the Frobenius norm is
1915 the root-of-sum-of-squares norm.
1917 Returns
1918 -------
1919 c : {float, inf}
1920 The condition number of the matrix. May be infinite.
1922 See Also
1923 --------
1924 numpy.linalg.norm
1926 Notes
1927 -----
1928 The condition number of `x` is defined as the norm of `x` times the
1929 norm of the inverse of `x` [1]_; the norm can be the usual L2-norm
1930 (root-of-sum-of-squares) or one of a number of other matrix norms.
1932 References
1933 ----------
1934 .. [1] G. Strang, *Linear Algebra and Its Applications*, Orlando, FL,
1935 Academic Press, Inc., 1980, pg. 285.
1937 Examples
1938 --------
1939 >>> import numpy as np
1940 >>> from numpy import linalg as LA
1941 >>> a = np.array([[1, 0, -1], [0, 1, 0], [1, 0, 1]])
1942 >>> a
1943 array([[ 1, 0, -1],
1944 [ 0, 1, 0],
1945 [ 1, 0, 1]])
1946 >>> LA.cond(a)
1947 1.4142135623730951
1948 >>> LA.cond(a, 'fro')
1949 3.1622776601683795
1950 >>> LA.cond(a, np.inf)
1951 2.0
1952 >>> LA.cond(a, -np.inf)
1953 1.0
1954 >>> LA.cond(a, 1)
1955 2.0
1956 >>> LA.cond(a, -1)
1957 1.0
1958 >>> LA.cond(a, 2)
1959 1.4142135623730951
1960 >>> LA.cond(a, -2)
1961 0.70710678118654746 # may vary
1962 >>> (min(LA.svd(a, compute_uv=False)) *
1963 ... min(LA.svd(LA.inv(a), compute_uv=False)))
1964 0.70710678118654746 # may vary
1966 """
1967 x = asarray(x) # in case we have a matrix
1968 if _is_empty_2d(x):
1969 raise LinAlgError("cond is not defined on empty arrays")
1970 if p is None or p == 2 or p == -2:
1971 s = svd(x, compute_uv=False)
1972 with errstate(all='ignore'):
1973 if p == -2:
1974 r = s[..., -1] / s[..., 0]
1975 else:
1976 r = s[..., 0] / s[..., -1]
1977 else:
1978 # Call inv(x) ignoring errors. The result array will
1979 # contain nans in the entries where inversion failed.
1980 _assert_stacked_2d(x)
1981 _assert_stacked_square(x)
1982 t, result_t = _commonType(x)
1983 signature = 'D->D' if isComplexType(t) else 'd->d'
1984 with errstate(all='ignore'):
1985 invx = _umath_linalg.inv(x, signature=signature)
1986 r = norm(x, p, axis=(-2, -1)) * norm(invx, p, axis=(-2, -1))
1987 r = r.astype(result_t, copy=False)
1989 # Convert nans to infs unless the original array had nan entries
1990 r = asarray(r)
1991 nan_mask = isnan(r)
1992 if nan_mask.any():
1993 nan_mask &= ~isnan(x).any(axis=(-2, -1))
1994 if r.ndim > 0:
1995 r[nan_mask] = inf
1996 elif nan_mask:
1997 r[()] = inf
1999 # Convention is to return scalars instead of 0d arrays
2000 if r.ndim == 0:
2001 r = r[()]
2003 return r
2006def _matrix_rank_dispatcher(A, tol=None, hermitian=None, *, rtol=None):
2007 return (A,)
2010@array_function_dispatch(_matrix_rank_dispatcher)
2011def matrix_rank(A, tol=None, hermitian=False, *, rtol=None):
2012 """
2013 Return matrix rank of array using SVD method
2015 Rank of the array is the number of singular values of the array that are
2016 greater than `tol`.
2018 Parameters
2019 ----------
2020 A : {(M,), (..., M, N)} array_like
2021 Input vector or stack of matrices.
2022 tol : (...) array_like, float, optional
2023 Threshold below which SVD values are considered zero. If `tol` is
2024 None, and ``S`` is an array with singular values for `M`, and
2025 ``eps`` is the epsilon value for datatype of ``S``, then `tol` is
2026 set to ``S.max() * max(M, N) * eps``.
2027 hermitian : bool, optional
2028 If True, `A` is assumed to be Hermitian (symmetric if real-valued),
2029 enabling a more efficient method for finding singular values.
2030 Defaults to False.
2031 rtol : (...) array_like, float, optional
2032 Parameter for the relative tolerance component. Only ``tol`` or
2033 ``rtol`` can be set at a time. Defaults to ``max(M, N) * eps``.
2035 .. versionadded:: 2.0.0
2037 Returns
2038 -------
2039 rank : (...) array_like
2040 Rank of A.
2042 Notes
2043 -----
2044 The default threshold to detect rank deficiency is a test on the magnitude
2045 of the singular values of `A`. By default, we identify singular values
2046 less than ``S.max() * max(M, N) * eps`` as indicating rank deficiency
2047 (with the symbols defined above). This is the algorithm MATLAB uses [1].
2048 It also appears in *Numerical recipes* in the discussion of SVD solutions
2049 for linear least squares [2].
2051 This default threshold is designed to detect rank deficiency accounting
2052 for the numerical errors of the SVD computation. Imagine that there
2053 is a column in `A` that is an exact (in floating point) linear combination
2054 of other columns in `A`. Computing the SVD on `A` will not produce
2055 a singular value exactly equal to 0 in general: any difference of
2056 the smallest SVD value from 0 will be caused by numerical imprecision
2057 in the calculation of the SVD. Our threshold for small SVD values takes
2058 this numerical imprecision into account, and the default threshold will
2059 detect such numerical rank deficiency. The threshold may declare a matrix
2060 `A` rank deficient even if the linear combination of some columns of `A`
2061 is not exactly equal to another column of `A` but only numerically very
2062 close to another column of `A`.
2064 We chose our default threshold because it is in wide use. Other thresholds
2065 are possible. For example, elsewhere in the 2007 edition of *Numerical
2066 recipes* there is an alternative threshold of ``S.max() *
2067 np.finfo(A.dtype).eps / 2. * np.sqrt(m + n + 1.)``. The authors describe
2068 this threshold as being based on "expected roundoff error" (p 71).
2070 The thresholds above deal with floating point roundoff error in the
2071 calculation of the SVD. However, you may have more information about
2072 the sources of error in `A` that would make you consider other tolerance
2073 values to detect *effective* rank deficiency. The most useful measure
2074 of the tolerance depends on the operations you intend to use on your
2075 matrix. For example, if your data come from uncertain measurements with
2076 uncertainties greater than floating point epsilon, choosing a tolerance
2077 near that uncertainty may be preferable. The tolerance may be absolute
2078 if the uncertainties are absolute rather than relative.
2080 References
2081 ----------
2082 .. [1] MATLAB reference documentation, "Rank"
2083 https://www.mathworks.com/help/techdoc/ref/rank.html
2084 .. [2] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery,
2085 "Numerical Recipes (3rd edition)", Cambridge University Press, 2007,
2086 page 795.
2088 Examples
2089 --------
2090 >>> import numpy as np
2091 >>> from numpy.linalg import matrix_rank
2092 >>> matrix_rank(np.eye(4)) # Full rank matrix
2093 4
2094 >>> I=np.eye(4); I[-1,-1] = 0. # rank deficient matrix
2095 >>> matrix_rank(I)
2096 3
2097 >>> matrix_rank(np.ones((4,))) # 1 dimension - rank 1 unless all 0
2098 1
2099 >>> matrix_rank(np.zeros((4,)))
2100 0
2101 """
2102 if rtol is not None and tol is not None:
2103 raise ValueError("`tol` and `rtol` can't be both set.")
2105 A = asarray(A)
2106 if A.ndim < 2:
2107 return int(not all(A == 0))
2108 S = svd(A, compute_uv=False, hermitian=hermitian)
2110 if tol is None:
2111 if rtol is None:
2112 rtol = max(A.shape[-2:]) * finfo(S.dtype).eps
2113 else:
2114 rtol = asarray(rtol)[..., newaxis]
2115 tol = S.max(axis=-1, keepdims=True) * rtol
2116 else:
2117 tol = asarray(tol)[..., newaxis]
2119 return count_nonzero(S > tol, axis=-1)
2122# Generalized inverse
2124def _pinv_dispatcher(a, rcond=None, hermitian=None, *, rtol=None):
2125 return (a,)
2128@array_function_dispatch(_pinv_dispatcher)
2129def pinv(a, rcond=None, hermitian=False, *, rtol=_NoValue):
2130 """
2131 Compute the (Moore-Penrose) pseudo-inverse of a matrix.
2133 Calculate the generalized inverse of a matrix using its
2134 singular-value decomposition (SVD) and including all
2135 *large* singular values.
2137 Parameters
2138 ----------
2139 a : (..., M, N) array_like
2140 Matrix or stack of matrices to be pseudo-inverted.
2141 rcond : (...) array_like of float, optional
2142 Cutoff for small singular values.
2143 Singular values less than or equal to
2144 ``rcond * largest_singular_value`` are set to zero.
2145 Broadcasts against the stack of matrices. Default: ``1e-15``.
2146 hermitian : bool, optional
2147 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
2148 enabling a more efficient method for finding singular values.
2149 Defaults to False.
2150 rtol : (...) array_like of float, optional
2151 Same as `rcond`, but it's an Array API compatible parameter name.
2152 Only `rcond` or `rtol` can be set at a time. If none of them are
2153 provided then NumPy's ``1e-15`` default is used. If ``rtol=None``
2154 is passed then the API standard default is used.
2156 .. versionadded:: 2.0.0
2158 Returns
2159 -------
2160 B : (..., N, M) ndarray
2161 The pseudo-inverse of `a`. If `a` is a `matrix` instance, then so
2162 is `B`.
2164 Raises
2165 ------
2166 LinAlgError
2167 If the SVD computation does not converge.
2169 See Also
2170 --------
2171 scipy.linalg.pinv : Similar function in SciPy.
2172 scipy.linalg.pinvh : Compute the (Moore-Penrose) pseudo-inverse of a
2173 Hermitian matrix.
2175 Notes
2176 -----
2177 The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
2178 defined as: "the matrix that 'solves' [the least-squares problem]
2179 :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then
2180 :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`.
2182 It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular
2183 value decomposition of A, then
2184 :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are
2185 orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting
2186 of A's so-called singular values, (followed, typically, by
2187 zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix
2188 consisting of the reciprocals of A's singular values
2189 (again, followed by zeros). [1]_
2191 References
2192 ----------
2193 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
2194 FL, Academic Press, Inc., 1980, pp. 139-142.
2196 Examples
2197 --------
2198 The following example checks that ``a * a+ * a == a`` and
2199 ``a+ * a * a+ == a+``:
2201 >>> import numpy as np
2202 >>> rng = np.random.default_rng()
2203 >>> a = rng.normal(size=(9, 6))
2204 >>> B = np.linalg.pinv(a)
2205 >>> np.allclose(a, np.dot(a, np.dot(B, a)))
2206 True
2207 >>> np.allclose(B, np.dot(B, np.dot(a, B)))
2208 True
2210 """
2211 a, wrap = _makearray(a)
2212 if rcond is None:
2213 if rtol is _NoValue:
2214 rcond = 1e-15
2215 elif rtol is None:
2216 rcond = max(a.shape[-2:]) * finfo(a.dtype).eps
2217 else:
2218 rcond = rtol
2219 elif rtol is not _NoValue:
2220 raise ValueError("`rtol` and `rcond` can't be both set.")
2221 else:
2222 # NOTE: Deprecate `rcond` in a few versions.
2223 pass
2225 rcond = asarray(rcond)
2226 if _is_empty_2d(a):
2227 m, n = a.shape[-2:]
2228 res = empty(a.shape[:-2] + (n, m), dtype=a.dtype)
2229 return wrap(res)
2230 a = a.conjugate()
2231 u, s, vt = svd(a, full_matrices=False, hermitian=hermitian)
2233 # discard small singular values
2234 cutoff = rcond[..., newaxis] * amax(s, axis=-1, keepdims=True)
2235 large = s > cutoff
2236 s = divide(1, s, where=large, out=s)
2237 s[~large] = 0
2239 res = matmul(transpose(vt), multiply(s[..., newaxis], transpose(u)))
2240 return wrap(res)
2243# Determinant
2246@array_function_dispatch(_unary_dispatcher)
2247def slogdet(a):
2248 """
2249 Compute the sign and (natural) logarithm of the determinant of an array.
2251 If an array has a very small or very large determinant, then a call to
2252 `det` may overflow or underflow. This routine is more robust against such
2253 issues, because it computes the logarithm of the determinant rather than
2254 the determinant itself.
2256 Parameters
2257 ----------
2258 a : (..., M, M) array_like
2259 Input array, has to be a square 2-D array.
2261 Returns
2262 -------
2263 A namedtuple with the following attributes:
2265 sign : (...) array_like
2266 A number representing the sign of the determinant. For a real matrix,
2267 this is 1, 0, or -1. For a complex matrix, this is a complex number
2268 with absolute value 1 (i.e., it is on the unit circle), or else 0.
2269 logabsdet : (...) array_like
2270 The natural log of the absolute value of the determinant.
2272 If the determinant is zero, then `sign` will be 0 and `logabsdet`
2273 will be -inf. In all cases, the determinant is equal to
2274 ``sign * np.exp(logabsdet)``.
2276 See Also
2277 --------
2278 det
2280 Notes
2281 -----
2282 Broadcasting rules apply, see the `numpy.linalg` documentation for
2283 details.
2285 The determinant is computed via LU factorization using the LAPACK
2286 routine ``z/dgetrf``.
2288 Examples
2289 --------
2290 The determinant of a 2-D array ``[[a, b], [c, d]]`` is ``ad - bc``:
2292 >>> import numpy as np
2293 >>> a = np.array([[1, 2], [3, 4]])
2294 >>> (sign, logabsdet) = np.linalg.slogdet(a)
2295 >>> (sign, logabsdet)
2296 (-1, 0.69314718055994529) # may vary
2297 >>> sign * np.exp(logabsdet)
2298 -2.0
2300 Computing log-determinants for a stack of matrices:
2302 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2303 >>> a.shape
2304 (3, 2, 2)
2305 >>> sign, logabsdet = np.linalg.slogdet(a)
2306 >>> (sign, logabsdet)
2307 (array([-1., -1., -1.]), array([ 0.69314718, 1.09861229, 2.07944154]))
2308 >>> sign * np.exp(logabsdet)
2309 array([-2., -3., -8.])
2311 This routine succeeds where ordinary `det` does not:
2313 >>> np.linalg.det(np.eye(500) * 0.1)
2314 0.0
2315 >>> np.linalg.slogdet(np.eye(500) * 0.1)
2316 (1, -1151.2925464970228)
2318 """
2319 a = asarray(a)
2320 _assert_stacked_2d(a)
2321 _assert_stacked_square(a)
2322 t, result_t = _commonType(a)
2323 real_t = _realType(result_t)
2324 signature = 'D->Dd' if isComplexType(t) else 'd->dd'
2325 sign, logdet = _umath_linalg.slogdet(a, signature=signature)
2326 sign = sign.astype(result_t, copy=False)
2327 logdet = logdet.astype(real_t, copy=False)
2328 return SlogdetResult(sign, logdet)
2331@array_function_dispatch(_unary_dispatcher)
2332def det(a):
2333 """
2334 Compute the determinant of an array.
2336 Parameters
2337 ----------
2338 a : (..., M, M) array_like
2339 Input array to compute determinants for.
2341 Returns
2342 -------
2343 det : (...) array_like
2344 Determinant of `a`.
2346 See Also
2347 --------
2348 slogdet : Another way to represent the determinant, more suitable
2349 for large matrices where underflow/overflow may occur.
2350 scipy.linalg.det : Similar function in SciPy.
2352 Notes
2353 -----
2354 Broadcasting rules apply, see the `numpy.linalg` documentation for
2355 details.
2357 The determinant is computed via LU factorization using the LAPACK
2358 routine ``z/dgetrf``.
2360 Examples
2361 --------
2362 The determinant of a 2-D array [[a, b], [c, d]] is ad - bc:
2364 >>> import numpy as np
2365 >>> a = np.array([[1, 2], [3, 4]])
2366 >>> np.linalg.det(a)
2367 -2.0 # may vary
2369 Computing determinants for a stack of matrices:
2371 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2372 >>> a.shape
2373 (3, 2, 2)
2374 >>> np.linalg.det(a)
2375 array([-2., -3., -8.])
2377 """
2378 a = asarray(a)
2379 _assert_stacked_2d(a)
2380 _assert_stacked_square(a)
2381 t, result_t = _commonType(a)
2382 signature = 'D->D' if isComplexType(t) else 'd->d'
2383 r = _umath_linalg.det(a, signature=signature)
2384 r = r.astype(result_t, copy=False)
2385 return r
2388# Linear Least Squares
2390def _lstsq_dispatcher(a, b, rcond=None):
2391 return (a, b)
2394@array_function_dispatch(_lstsq_dispatcher)
2395def lstsq(a, b, rcond=None):
2396 r"""
2397 Return the least-squares solution to a linear matrix equation.
2399 Computes the vector `x` that approximately solves the equation
2400 ``a @ x = b``. The equation may be under-, well-, or over-determined
2401 (i.e., the number of linearly independent rows of `a` can be less than,
2402 equal to, or greater than its number of linearly independent columns).
2403 If `a` is square and of full rank, then `x` (but for round-off error)
2404 is the "exact" solution of the equation. Else, `x` minimizes the
2405 Euclidean 2-norm :math:`||b - ax||`. If there are multiple minimizing
2406 solutions, the one with the smallest 2-norm :math:`||x||` is returned.
2408 Parameters
2409 ----------
2410 a : (M, N) array_like
2411 "Coefficient" matrix.
2412 b : {(M,), (M, K)} array_like
2413 Ordinate or "dependent variable" values. If `b` is two-dimensional,
2414 the least-squares solution is calculated for each of the `K` columns
2415 of `b`.
2416 rcond : float, optional
2417 Cut-off ratio for small singular values of `a`.
2418 For the purposes of rank determination, singular values are treated
2419 as zero if they are smaller than `rcond` times the largest singular
2420 value of `a`.
2421 The default uses the machine precision times ``max(M, N)``. Passing
2422 ``-1`` will use machine precision.
2424 .. versionchanged:: 2.0
2425 Previously, the default was ``-1``, but a warning was given that
2426 this would change.
2428 Returns
2429 -------
2430 x : {(N,), (N, K)} ndarray
2431 Least-squares solution. If `b` is two-dimensional,
2432 the solutions are in the `K` columns of `x`.
2433 residuals : {(1,), (K,), (0,)} ndarray
2434 Sums of squared residuals: Squared Euclidean 2-norm for each column in
2435 ``b - a @ x``.
2436 If the rank of `a` is < N or M <= N, this is an empty array.
2437 If `b` is 1-dimensional, this is a (1,) shape array.
2438 Otherwise the shape is (K,).
2439 rank : int
2440 Rank of matrix `a`.
2441 s : (min(M, N),) ndarray
2442 Singular values of `a`.
2444 Raises
2445 ------
2446 LinAlgError
2447 If computation does not converge.
2449 See Also
2450 --------
2451 scipy.linalg.lstsq : Similar function in SciPy.
2453 Notes
2454 -----
2455 If `b` is a matrix, then all array results are returned as matrices.
2457 Examples
2458 --------
2459 Fit a line, ``y = mx + c``, through some noisy data-points:
2461 >>> import numpy as np
2462 >>> x = np.array([0, 1, 2, 3])
2463 >>> y = np.array([-1, 0.2, 0.9, 2.1])
2465 By examining the coefficients, we see that the line should have a
2466 gradient of roughly 1 and cut the y-axis at, more or less, -1.
2468 We can rewrite the line equation as ``y = Ap``, where ``A = [[x 1]]``
2469 and ``p = [[m], [c]]``. Now use `lstsq` to solve for `p`:
2471 >>> A = np.vstack([x, np.ones(len(x))]).T
2472 >>> A
2473 array([[ 0., 1.],
2474 [ 1., 1.],
2475 [ 2., 1.],
2476 [ 3., 1.]])
2478 >>> m, c = np.linalg.lstsq(A, y)[0]
2479 >>> m, c
2480 (1.0 -0.95) # may vary
2482 Plot the data along with the fitted line:
2484 >>> import matplotlib.pyplot as plt
2485 >>> _ = plt.plot(x, y, 'o', label='Original data', markersize=10)
2486 >>> _ = plt.plot(x, m*x + c, 'r', label='Fitted line')
2487 >>> _ = plt.legend()
2488 >>> plt.show()
2490 """
2491 a, _ = _makearray(a)
2492 b, wrap = _makearray(b)
2493 is_1d = b.ndim == 1
2494 if is_1d:
2495 b = b[:, newaxis]
2496 _assert_2d(a, b)
2497 m, n = a.shape[-2:]
2498 m2, n_rhs = b.shape[-2:]
2499 if m != m2:
2500 raise LinAlgError('Incompatible dimensions')
2502 t, result_t = _commonType(a, b)
2503 result_real_t = _realType(result_t)
2505 if rcond is None:
2506 rcond = finfo(t).eps * max(n, m)
2508 signature = 'DDd->Ddid' if isComplexType(t) else 'ddd->ddid'
2509 if n_rhs == 0:
2510 # lapack can't handle n_rhs = 0 - so allocate
2511 # the array one larger in that axis
2512 b = zeros(b.shape[:-2] + (m, n_rhs + 1), dtype=b.dtype)
2514 with errstate(call=_raise_linalgerror_lstsq, invalid='call',
2515 over='ignore', divide='ignore', under='ignore'):
2516 x, resids, rank, s = _umath_linalg.lstsq(a, b, rcond,
2517 signature=signature)
2518 if m == 0:
2519 x[...] = 0
2520 if n_rhs == 0:
2521 # remove the item we added
2522 x = x[..., :n_rhs]
2523 resids = resids[..., :n_rhs]
2525 # remove the axis we added
2526 if is_1d:
2527 x = x.squeeze(axis=-1)
2528 # we probably should squeeze resids too, but we can't
2529 # without breaking compatibility.
2531 # as documented
2532 if rank != n or m <= n:
2533 resids = array([], result_real_t)
2535 # coerce output arrays
2536 s = s.astype(result_real_t, copy=False)
2537 resids = resids.astype(result_real_t, copy=False)
2538 # Copying lets the memory in r_parts be freed
2539 x = x.astype(result_t, copy=True)
2540 return wrap(x), wrap(resids), rank, s
2543def _multi_svd_norm(x, row_axis, col_axis, op):
2544 """Compute a function of the singular values of the 2-D matrices in `x`.
2546 This is a private utility function used by `numpy.linalg.norm()`.
2548 Parameters
2549 ----------
2550 x : ndarray
2551 row_axis, col_axis : int
2552 The axes of `x` that hold the 2-D matrices.
2553 op : callable
2554 This should be either numpy.amin or `numpy.amax` or `numpy.sum`.
2556 Returns
2557 -------
2558 result : float or ndarray
2559 If `x` is 2-D, the return values is a float.
2560 Otherwise, it is an array with ``x.ndim - 2`` dimensions.
2561 The return values are either the minimum or maximum or sum of the
2562 singular values of the matrices, depending on whether `op`
2563 is `numpy.amin` or `numpy.amax` or `numpy.sum`.
2565 """
2566 y = moveaxis(x, (row_axis, col_axis), (-2, -1))
2567 result = op(svd(y, compute_uv=False), axis=-1)
2568 return result
2571def _norm_dispatcher(x, ord=None, axis=None, keepdims=None):
2572 return (x,)
2575@array_function_dispatch(_norm_dispatcher)
2576def norm(x, ord=None, axis=None, keepdims=False):
2577 """
2578 Matrix or vector norm.
2580 This function is able to return one of eight different matrix norms,
2581 or one of an infinite number of vector norms (described below), depending
2582 on the value of the ``ord`` parameter.
2584 Parameters
2585 ----------
2586 x : array_like
2587 Input array. If `axis` is None, `x` must be 1-D or 2-D, unless `ord`
2588 is None. If both `axis` and `ord` are None, the 2-norm of
2589 ``x.ravel`` will be returned.
2590 ord : {int, float, inf, -inf, 'fro', 'nuc'}, optional
2591 Order of the norm (see table under ``Notes`` for what values are
2592 supported for matrices and vectors respectively). inf means numpy's
2593 `inf` object. The default is None.
2594 axis : {None, int, 2-tuple of ints}, optional.
2595 If `axis` is an integer, it specifies the axis of `x` along which to
2596 compute the vector norms. If `axis` is a 2-tuple, it specifies the
2597 axes that hold 2-D matrices, and the matrix norms of these matrices
2598 are computed. If `axis` is None then either a vector norm (when `x`
2599 is 1-D) or a matrix norm (when `x` is 2-D) is returned. The default
2600 is None.
2602 keepdims : bool, optional
2603 If this is set to True, the axes which are normed over are left in the
2604 result as dimensions with size one. With this option the result will
2605 broadcast correctly against the original `x`.
2607 Returns
2608 -------
2609 n : float or ndarray
2610 Norm of the matrix or vector(s).
2612 See Also
2613 --------
2614 scipy.linalg.norm : Similar function in SciPy.
2616 Notes
2617 -----
2618 For values of ``ord < 1``, the result is, strictly speaking, not a
2619 mathematical 'norm', but it may still be useful for various numerical
2620 purposes.
2622 The following norms can be calculated:
2624 ===== ============================ ==========================
2625 ord norm for matrices norm for vectors
2626 ===== ============================ ==========================
2627 None Frobenius norm 2-norm
2628 'fro' Frobenius norm --
2629 'nuc' nuclear norm --
2630 inf max(sum(abs(x), axis=1)) max(abs(x))
2631 -inf min(sum(abs(x), axis=1)) min(abs(x))
2632 0 -- sum(x != 0)
2633 1 max(sum(abs(x), axis=0)) as below
2634 -1 min(sum(abs(x), axis=0)) as below
2635 2 2-norm (largest sing. value) as below
2636 -2 smallest singular value as below
2637 other -- sum(abs(x)**ord)**(1./ord)
2638 ===== ============================ ==========================
2640 The Frobenius norm is given by [1]_:
2642 :math:`||A||_F = [\\sum_{i,j} abs(a_{i,j})^2]^{1/2}`
2644 The nuclear norm is the sum of the singular values.
2646 Both the Frobenius and nuclear norm orders are only defined for
2647 matrices and raise a ValueError when ``x.ndim != 2``.
2649 References
2650 ----------
2651 .. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*,
2652 Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15
2654 Examples
2655 --------
2657 >>> import numpy as np
2658 >>> from numpy import linalg as LA
2659 >>> a = np.arange(9) - 4
2660 >>> a
2661 array([-4, -3, -2, ..., 2, 3, 4])
2662 >>> b = a.reshape((3, 3))
2663 >>> b
2664 array([[-4, -3, -2],
2665 [-1, 0, 1],
2666 [ 2, 3, 4]])
2668 >>> LA.norm(a)
2669 7.745966692414834
2670 >>> LA.norm(b)
2671 7.745966692414834
2672 >>> LA.norm(b, 'fro')
2673 7.745966692414834
2674 >>> LA.norm(a, np.inf)
2675 4.0
2676 >>> LA.norm(b, np.inf)
2677 9.0
2678 >>> LA.norm(a, -np.inf)
2679 0.0
2680 >>> LA.norm(b, -np.inf)
2681 2.0
2683 >>> LA.norm(a, 1)
2684 20.0
2685 >>> LA.norm(b, 1)
2686 7.0
2687 >>> LA.norm(a, -1)
2688 -4.6566128774142013e-010
2689 >>> LA.norm(b, -1)
2690 6.0
2691 >>> LA.norm(a, 2)
2692 7.745966692414834
2693 >>> LA.norm(b, 2)
2694 7.3484692283495345
2696 >>> LA.norm(a, -2)
2697 0.0
2698 >>> LA.norm(b, -2)
2699 1.8570331885190563e-016 # may vary
2700 >>> LA.norm(a, 3)
2701 5.8480354764257312 # may vary
2702 >>> LA.norm(a, -3)
2703 0.0
2705 Using the `axis` argument to compute vector norms:
2707 >>> c = np.array([[ 1, 2, 3],
2708 ... [-1, 1, 4]])
2709 >>> LA.norm(c, axis=0)
2710 array([ 1.41421356, 2.23606798, 5. ])
2711 >>> LA.norm(c, axis=1)
2712 array([ 3.74165739, 4.24264069])
2713 >>> LA.norm(c, ord=1, axis=1)
2714 array([ 6., 6.])
2716 Using the `axis` argument to compute matrix norms:
2718 >>> m = np.arange(8).reshape(2,2,2)
2719 >>> LA.norm(m, axis=(1,2))
2720 array([ 3.74165739, 11.22497216])
2721 >>> LA.norm(m[0, :, :]), LA.norm(m[1, :, :])
2722 (3.7416573867739413, 11.224972160321824)
2724 """
2725 x = asarray(x)
2727 if not issubclass(x.dtype.type, (inexact, object_)):
2728 x = x.astype(float)
2730 # Immediately handle some default, simple, fast, and common cases.
2731 if axis is None:
2732 ndim = x.ndim
2733 if (
2734 (ord is None) or
2735 (ord in ('f', 'fro') and ndim == 2) or
2736 (ord == 2 and ndim == 1)
2737 ):
2738 x = x.ravel(order='K')
2739 if isComplexType(x.dtype.type):
2740 x_real = x.real
2741 x_imag = x.imag
2742 sqnorm = x_real.dot(x_real) + x_imag.dot(x_imag)
2743 else:
2744 sqnorm = x.dot(x)
2745 ret = sqrt(sqnorm)
2746 if keepdims:
2747 ret = ret.reshape(ndim*[1])
2748 return ret
2750 # Normalize the `axis` argument to a tuple.
2751 nd = x.ndim
2752 if axis is None:
2753 axis = tuple(range(nd))
2754 elif not isinstance(axis, tuple):
2755 try:
2756 axis = int(axis)
2757 except Exception as e:
2758 raise TypeError(
2759 "'axis' must be None, an integer or a tuple of integers"
2760 ) from e
2761 axis = (axis,)
2763 if len(axis) == 1:
2764 if ord == inf:
2765 return abs(x).max(axis=axis, keepdims=keepdims)
2766 elif ord == -inf:
2767 return abs(x).min(axis=axis, keepdims=keepdims)
2768 elif ord == 0:
2769 # Zero norm
2770 return (
2771 (x != 0)
2772 .astype(x.real.dtype)
2773 .sum(axis=axis, keepdims=keepdims)
2774 )
2775 elif ord == 1:
2776 # special case for speedup
2777 return add.reduce(abs(x), axis=axis, keepdims=keepdims)
2778 elif ord is None or ord == 2:
2779 # special case for speedup
2780 s = (x.conj() * x).real
2781 return sqrt(add.reduce(s, axis=axis, keepdims=keepdims))
2782 # None of the str-type keywords for ord ('fro', 'nuc')
2783 # are valid for vectors
2784 elif isinstance(ord, str):
2785 raise ValueError(f"Invalid norm order '{ord}' for vectors")
2786 else:
2787 absx = abs(x)
2788 absx **= ord
2789 ret = add.reduce(absx, axis=axis, keepdims=keepdims)
2790 ret **= reciprocal(ord, dtype=ret.dtype)
2791 return ret
2792 elif len(axis) == 2:
2793 row_axis, col_axis = axis
2794 row_axis = normalize_axis_index(row_axis, nd)
2795 col_axis = normalize_axis_index(col_axis, nd)
2796 if row_axis == col_axis:
2797 raise ValueError('Duplicate axes given.')
2798 if ord == 2:
2799 ret = _multi_svd_norm(x, row_axis, col_axis, amax)
2800 elif ord == -2:
2801 ret = _multi_svd_norm(x, row_axis, col_axis, amin)
2802 elif ord == 1:
2803 if col_axis > row_axis:
2804 col_axis -= 1
2805 ret = add.reduce(abs(x), axis=row_axis).max(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).max(axis=row_axis)
2810 elif ord == -1:
2811 if col_axis > row_axis:
2812 col_axis -= 1
2813 ret = add.reduce(abs(x), axis=row_axis).min(axis=col_axis)
2814 elif ord == -inf:
2815 if row_axis > col_axis:
2816 row_axis -= 1
2817 ret = add.reduce(abs(x), axis=col_axis).min(axis=row_axis)
2818 elif ord in [None, 'fro', 'f']:
2819 ret = sqrt(add.reduce((x.conj() * x).real, axis=axis))
2820 elif ord == 'nuc':
2821 ret = _multi_svd_norm(x, row_axis, col_axis, sum)
2822 else:
2823 raise ValueError("Invalid norm order for matrices.")
2824 if keepdims:
2825 ret_shape = list(x.shape)
2826 ret_shape[axis[0]] = 1
2827 ret_shape[axis[1]] = 1
2828 ret = ret.reshape(ret_shape)
2829 return ret
2830 else:
2831 raise ValueError("Improper number of dimensions to norm.")
2834# multi_dot
2836def _multidot_dispatcher(arrays, *, out=None):
2837 yield from arrays
2838 yield out
2841@array_function_dispatch(_multidot_dispatcher)
2842def multi_dot(arrays, *, out=None):
2843 """
2844 Compute the dot product of two or more arrays in a single function call,
2845 while automatically selecting the fastest evaluation order.
2847 `multi_dot` chains `numpy.dot` and uses optimal parenthesization
2848 of the matrices [1]_ [2]_. Depending on the shapes of the matrices,
2849 this can speed up the multiplication a lot.
2851 If the first argument is 1-D it is treated as a row vector.
2852 If the last argument is 1-D it is treated as a column vector.
2853 The other arguments must be 2-D.
2855 Think of `multi_dot` as::
2857 def multi_dot(arrays): return functools.reduce(np.dot, arrays)
2860 Parameters
2861 ----------
2862 arrays : sequence of array_like
2863 If the first argument is 1-D it is treated as row vector.
2864 If the last argument is 1-D it is treated as column vector.
2865 The other arguments must be 2-D.
2866 out : ndarray, optional
2867 Output argument. This must have the exact kind that would be returned
2868 if it was not used. In particular, it must have the right type, must be
2869 C-contiguous, and its dtype must be the dtype that would be returned
2870 for `dot(a, b)`. This is a performance feature. Therefore, if these
2871 conditions are not met, an exception is raised, instead of attempting
2872 to be flexible.
2874 Returns
2875 -------
2876 output : ndarray
2877 Returns the dot product of the supplied arrays.
2879 See Also
2880 --------
2881 numpy.dot : dot multiplication with two arguments.
2883 References
2884 ----------
2886 .. [1] Cormen, "Introduction to Algorithms", Chapter 15.2, p. 370-378
2887 .. [2] https://en.wikipedia.org/wiki/Matrix_chain_multiplication
2889 Examples
2890 --------
2891 `multi_dot` allows you to write::
2893 >>> import numpy as np
2894 >>> from numpy.linalg import multi_dot
2895 >>> # Prepare some data
2896 >>> A = np.random.random((10000, 100))
2897 >>> B = np.random.random((100, 1000))
2898 >>> C = np.random.random((1000, 5))
2899 >>> D = np.random.random((5, 333))
2900 >>> # the actual dot multiplication
2901 >>> _ = multi_dot([A, B, C, D])
2903 instead of::
2905 >>> _ = np.dot(np.dot(np.dot(A, B), C), D)
2906 >>> # or
2907 >>> _ = A.dot(B).dot(C).dot(D)
2909 Notes
2910 -----
2911 The cost for a matrix multiplication can be calculated with the
2912 following function::
2914 def cost(A, B):
2915 return A.shape[0] * A.shape[1] * B.shape[1]
2917 Assume we have three matrices
2918 :math:`A_{10x100}, B_{100x5}, C_{5x50}`.
2920 The costs for the two different parenthesizations are as follows::
2922 cost((AB)C) = 10*100*5 + 10*5*50 = 5000 + 2500 = 7500
2923 cost(A(BC)) = 10*100*50 + 100*5*50 = 50000 + 25000 = 75000
2925 """
2926 n = len(arrays)
2927 # optimization only makes sense for len(arrays) > 2
2928 if n < 2:
2929 raise ValueError("Expecting at least two arrays.")
2930 elif n == 2:
2931 return dot(arrays[0], arrays[1], out=out)
2933 arrays = [asanyarray(a) for a in arrays]
2935 # save original ndim to reshape the result array into the proper form later
2936 ndim_first, ndim_last = arrays[0].ndim, arrays[-1].ndim
2937 # Explicitly convert vectors to 2D arrays to keep the logic of the internal
2938 # _multi_dot_* functions as simple as possible.
2939 if arrays[0].ndim == 1:
2940 arrays[0] = atleast_2d(arrays[0])
2941 if arrays[-1].ndim == 1:
2942 arrays[-1] = atleast_2d(arrays[-1]).T
2943 _assert_2d(*arrays)
2945 # _multi_dot_three is much faster than _multi_dot_matrix_chain_order
2946 if n == 3:
2947 result = _multi_dot_three(arrays[0], arrays[1], arrays[2], out=out)
2948 else:
2949 order = _multi_dot_matrix_chain_order(arrays)
2950 result = _multi_dot(arrays, order, 0, n - 1, out=out)
2952 # return proper shape
2953 if ndim_first == 1 and ndim_last == 1:
2954 return result[0, 0] # scalar
2955 elif ndim_first == 1 or ndim_last == 1:
2956 return result.ravel() # 1-D
2957 else:
2958 return result
2961def _multi_dot_three(A, B, C, out=None):
2962 """
2963 Find the best order for three arrays and do the multiplication.
2965 For three arguments `_multi_dot_three` is approximately 15 times faster
2966 than `_multi_dot_matrix_chain_order`
2968 """
2969 a0, a1b0 = A.shape
2970 b1c0, c1 = C.shape
2971 # cost1 = cost((AB)C) = a0*a1b0*b1c0 + a0*b1c0*c1
2972 cost1 = a0 * b1c0 * (a1b0 + c1)
2973 # cost2 = cost(A(BC)) = a1b0*b1c0*c1 + a0*a1b0*c1
2974 cost2 = a1b0 * c1 * (a0 + b1c0)
2976 if cost1 < cost2:
2977 return dot(dot(A, B), C, out=out)
2978 else:
2979 return dot(A, dot(B, C), out=out)
2982def _multi_dot_matrix_chain_order(arrays, return_costs=False):
2983 """
2984 Return a np.array that encodes the optimal order of multiplications.
2986 The optimal order array is then used by `_multi_dot()` to do the
2987 multiplication.
2989 Also return the cost matrix if `return_costs` is `True`
2991 The implementation CLOSELY follows Cormen, "Introduction to Algorithms",
2992 Chapter 15.2, p. 370-378. Note that Cormen uses 1-based indices.
2994 cost[i, j] = min([
2995 cost[prefix] + cost[suffix] + cost_mult(prefix, suffix)
2996 for k in range(i, j)])
2998 """
2999 n = len(arrays)
3000 # p stores the dimensions of the matrices
3001 # Example for p: A_{10x100}, B_{100x5}, C_{5x50} --> p = [10, 100, 5, 50]
3002 p = [a.shape[0] for a in arrays] + [arrays[-1].shape[1]]
3003 # m is a matrix of costs of the subproblems
3004 # m[i,j]: min number of scalar multiplications needed to compute A_{i..j}
3005 m = zeros((n, n), dtype=double)
3006 # s is the actual ordering
3007 # s[i, j] is the value of k at which we split the product A_i..A_j
3008 s = empty((n, n), dtype=intp)
3010 for l in range(1, n):
3011 for i in range(n - l):
3012 j = i + l
3013 m[i, j] = inf
3014 for k in range(i, j):
3015 q = m[i, k] + m[k+1, j] + p[i]*p[k+1]*p[j+1]
3016 if q < m[i, j]:
3017 m[i, j] = q
3018 s[i, j] = k # Note that Cormen uses 1-based index
3020 return (s, m) if return_costs else s
3023def _multi_dot(arrays, order, i, j, out=None):
3024 """Actually do the multiplication with the given order."""
3025 if i == j:
3026 # the initial call with non-None out should never get here
3027 assert out is None
3029 return arrays[i]
3030 else:
3031 return dot(_multi_dot(arrays, order, i, order[i, j]),
3032 _multi_dot(arrays, order, order[i, j] + 1, j),
3033 out=out)
3036# diagonal
3038def _diagonal_dispatcher(x, /, *, offset=None):
3039 return (x,)
3042@array_function_dispatch(_diagonal_dispatcher)
3043def diagonal(x, /, *, offset=0):
3044 """
3045 Returns specified diagonals of a matrix (or a stack of matrices) ``x``.
3047 This function is Array API compatible, contrary to
3048 :py:func:`numpy.diagonal`, the matrix is assumed
3049 to be defined by the last two dimensions.
3051 Parameters
3052 ----------
3053 x : (...,M,N) array_like
3054 Input array having shape (..., M, N) and whose innermost two
3055 dimensions form MxN matrices.
3056 offset : int, optional
3057 Offset specifying the off-diagonal relative to the main diagonal,
3058 where::
3060 * offset = 0: the main diagonal.
3061 * offset > 0: off-diagonal above the main diagonal.
3062 * offset < 0: off-diagonal below the main diagonal.
3064 Returns
3065 -------
3066 out : (...,min(N,M)) ndarray
3067 An array containing the diagonals and whose shape is determined by
3068 removing the last two dimensions and appending a dimension equal to
3069 the size of the resulting diagonals. The returned array must have
3070 the same data type as ``x``.
3072 See Also
3073 --------
3074 numpy.diagonal
3076 Examples
3077 --------
3078 >>> a = np.arange(4).reshape(2, 2); a
3079 array([[0, 1],
3080 [2, 3]])
3081 >>> np.linalg.diagonal(a)
3082 array([0, 3])
3084 A 3-D example:
3086 >>> a = np.arange(8).reshape(2, 2, 2); a
3087 array([[[0, 1],
3088 [2, 3]],
3089 [[4, 5],
3090 [6, 7]]])
3091 >>> np.linalg.diagonal(a)
3092 array([[0, 3],
3093 [4, 7]])
3095 Diagonals adjacent to the main diagonal can be obtained by using the
3096 `offset` argument:
3098 >>> a = np.arange(9).reshape(3, 3)
3099 >>> a
3100 array([[0, 1, 2],
3101 [3, 4, 5],
3102 [6, 7, 8]])
3103 >>> np.linalg.diagonal(a, offset=1) # First superdiagonal
3104 array([1, 5])
3105 >>> np.linalg.diagonal(a, offset=2) # Second superdiagonal
3106 array([2])
3107 >>> np.linalg.diagonal(a, offset=-1) # First subdiagonal
3108 array([3, 7])
3109 >>> np.linalg.diagonal(a, offset=-2) # Second subdiagonal
3110 array([6])
3112 The anti-diagonal can be obtained by reversing the order of elements
3113 using either `numpy.flipud` or `numpy.fliplr`.
3115 >>> a = np.arange(9).reshape(3, 3)
3116 >>> a
3117 array([[0, 1, 2],
3118 [3, 4, 5],
3119 [6, 7, 8]])
3120 >>> np.linalg.diagonal(np.fliplr(a)) # Horizontal flip
3121 array([2, 4, 6])
3122 >>> np.linalg.diagonal(np.flipud(a)) # Vertical flip
3123 array([6, 4, 2])
3125 Note that the order in which the diagonal is retrieved varies depending
3126 on the flip function.
3128 """
3129 return _core_diagonal(x, offset, axis1=-2, axis2=-1)
3132# trace
3134def _trace_dispatcher(x, /, *, offset=None, dtype=None):
3135 return (x,)
3138@array_function_dispatch(_trace_dispatcher)
3139def trace(x, /, *, offset=0, dtype=None):
3140 """
3141 Returns the sum along the specified diagonals of a matrix
3142 (or a stack of matrices) ``x``.
3144 This function is Array API compatible, contrary to
3145 :py:func:`numpy.trace`.
3147 Parameters
3148 ----------
3149 x : (...,M,N) array_like
3150 Input array having shape (..., M, N) and whose innermost two
3151 dimensions form MxN matrices.
3152 offset : int, optional
3153 Offset specifying the off-diagonal relative to the main diagonal,
3154 where::
3156 * offset = 0: the main diagonal.
3157 * offset > 0: off-diagonal above the main diagonal.
3158 * offset < 0: off-diagonal below the main diagonal.
3160 dtype : dtype, optional
3161 Data type of the returned array.
3163 Returns
3164 -------
3165 out : ndarray
3166 An array containing the traces and whose shape is determined by
3167 removing the last two dimensions and storing the traces in the last
3168 array dimension. For example, if x has rank k and shape:
3169 (I, J, K, ..., L, M, N), then an output array has rank k-2 and shape:
3170 (I, J, K, ..., L) where::
3172 out[i, j, k, ..., l] = trace(a[i, j, k, ..., l, :, :])
3174 The returned array must have a data type as described by the dtype
3175 parameter above.
3177 See Also
3178 --------
3179 numpy.trace
3181 Examples
3182 --------
3183 >>> np.linalg.trace(np.eye(3))
3184 3.0
3185 >>> a = np.arange(8).reshape((2, 2, 2))
3186 >>> np.linalg.trace(a)
3187 array([3, 11])
3189 Trace is computed with the last two axes as the 2-d sub-arrays.
3190 This behavior differs from :py:func:`numpy.trace` which uses the first two
3191 axes by default.
3193 >>> a = np.arange(24).reshape((3, 2, 2, 2))
3194 >>> np.linalg.trace(a).shape
3195 (3, 2)
3197 Traces adjacent to the main diagonal can be obtained by using the
3198 `offset` argument:
3200 >>> a = np.arange(9).reshape((3, 3)); a
3201 array([[0, 1, 2],
3202 [3, 4, 5],
3203 [6, 7, 8]])
3204 >>> np.linalg.trace(a, offset=1) # First superdiagonal
3205 6
3206 >>> np.linalg.trace(a, offset=2) # Second superdiagonal
3207 2
3208 >>> np.linalg.trace(a, offset=-1) # First subdiagonal
3209 10
3210 >>> np.linalg.trace(a, offset=-2) # Second subdiagonal
3211 6
3213 """
3214 return _core_trace(x, offset, axis1=-2, axis2=-1, dtype=dtype)
3217# cross
3219def _cross_dispatcher(x1, x2, /, *, axis=None):
3220 return (x1, x2,)
3223@array_function_dispatch(_cross_dispatcher)
3224def cross(x1, x2, /, *, axis=-1):
3225 """
3226 Returns the cross product of 3-element vectors.
3228 If ``x1`` and/or ``x2`` are multi-dimensional arrays, then
3229 the cross-product of each pair of corresponding 3-element vectors
3230 is independently computed.
3232 This function is Array API compatible, contrary to
3233 :func:`numpy.cross`.
3235 Parameters
3236 ----------
3237 x1 : array_like
3238 The first input array.
3239 x2 : array_like
3240 The second input array. Must be compatible with ``x1`` for all
3241 non-compute axes. The size of the axis over which to compute
3242 the cross-product must be the same size as the respective axis
3243 in ``x1``.
3244 axis : int, optional
3245 The axis (dimension) of ``x1`` and ``x2`` containing the vectors for
3246 which to compute the cross-product. Default: ``-1``.
3248 Returns
3249 -------
3250 out : ndarray
3251 An array containing the cross products.
3253 See Also
3254 --------
3255 numpy.cross
3257 Examples
3258 --------
3259 Vector cross-product.
3261 >>> x = np.array([1, 2, 3])
3262 >>> y = np.array([4, 5, 6])
3263 >>> np.linalg.cross(x, y)
3264 array([-3, 6, -3])
3266 Multiple vector cross-products. Note that the direction of the cross
3267 product vector is defined by the *right-hand rule*.
3269 >>> x = np.array([[1,2,3], [4,5,6]])
3270 >>> y = np.array([[4,5,6], [1,2,3]])
3271 >>> np.linalg.cross(x, y)
3272 array([[-3, 6, -3],
3273 [ 3, -6, 3]])
3275 >>> x = np.array([[1, 2], [3, 4], [5, 6]])
3276 >>> y = np.array([[4, 5], [6, 1], [2, 3]])
3277 >>> np.linalg.cross(x, y, axis=0)
3278 array([[-24, 6],
3279 [ 18, 24],
3280 [-6, -18]])
3282 """
3283 x1 = asanyarray(x1)
3284 x2 = asanyarray(x2)
3286 if x1.shape[axis] != 3 or x2.shape[axis] != 3:
3287 raise ValueError(
3288 "Both input arrays must be (arrays of) 3-dimensional vectors, "
3289 f"but they are {x1.shape[axis]} and {x2.shape[axis]} "
3290 "dimensional instead."
3291 )
3293 return _core_cross(x1, x2, axis=axis)
3296# matmul
3298def _matmul_dispatcher(x1, x2, /):
3299 return (x1, x2)
3302@array_function_dispatch(_matmul_dispatcher)
3303def matmul(x1, x2, /):
3304 """
3305 Computes the matrix product.
3307 This function is Array API compatible, contrary to
3308 :func:`numpy.matmul`.
3310 Parameters
3311 ----------
3312 x1 : array_like
3313 The first input array.
3314 x2 : array_like
3315 The second input array.
3317 Returns
3318 -------
3319 out : ndarray
3320 The matrix product of the inputs.
3321 This is a scalar only when both ``x1``, ``x2`` are 1-d vectors.
3323 Raises
3324 ------
3325 ValueError
3326 If the last dimension of ``x1`` is not the same size as
3327 the second-to-last dimension of ``x2``.
3329 If a scalar value is passed in.
3331 See Also
3332 --------
3333 numpy.matmul
3335 Examples
3336 --------
3337 For 2-D arrays it is the matrix product:
3339 >>> a = np.array([[1, 0],
3340 ... [0, 1]])
3341 >>> b = np.array([[4, 1],
3342 ... [2, 2]])
3343 >>> np.linalg.matmul(a, b)
3344 array([[4, 1],
3345 [2, 2]])
3347 For 2-D mixed with 1-D, the result is the usual.
3349 >>> a = np.array([[1, 0],
3350 ... [0, 1]])
3351 >>> b = np.array([1, 2])
3352 >>> np.linalg.matmul(a, b)
3353 array([1, 2])
3354 >>> np.linalg.matmul(b, a)
3355 array([1, 2])
3358 Broadcasting is conventional for stacks of arrays
3360 >>> a = np.arange(2 * 2 * 4).reshape((2, 2, 4))
3361 >>> b = np.arange(2 * 2 * 4).reshape((2, 4, 2))
3362 >>> np.linalg.matmul(a,b).shape
3363 (2, 2, 2)
3364 >>> np.linalg.matmul(a, b)[0, 1, 1]
3365 98
3366 >>> sum(a[0, 1, :] * b[0 , :, 1])
3367 98
3369 Vector, vector returns the scalar inner product, but neither argument
3370 is complex-conjugated:
3372 >>> np.linalg.matmul([2j, 3j], [2j, 3j])
3373 (-13+0j)
3375 Scalar multiplication raises an error.
3377 >>> np.linalg.matmul([1,2], 3)
3378 Traceback (most recent call last):
3379 ...
3380 ValueError: matmul: Input operand 1 does not have enough dimensions ...
3382 """
3383 return _core_matmul(x1, x2)
3386# tensordot
3388def _tensordot_dispatcher(x1, x2, /, *, axes=None):
3389 return (x1, x2)
3392@array_function_dispatch(_tensordot_dispatcher)
3393def tensordot(x1, x2, /, *, axes=2):
3394 return _core_tensordot(x1, x2, axes=axes)
3397tensordot.__doc__ = _core_tensordot.__doc__
3400# matrix_transpose
3402def _matrix_transpose_dispatcher(x):
3403 return (x,)
3405@array_function_dispatch(_matrix_transpose_dispatcher)
3406def matrix_transpose(x, /):
3407 return _core_matrix_transpose(x)
3410matrix_transpose.__doc__ = _core_matrix_transpose.__doc__
3413# matrix_norm
3415def _matrix_norm_dispatcher(x, /, *, keepdims=None, ord=None):
3416 return (x,)
3418@array_function_dispatch(_matrix_norm_dispatcher)
3419def matrix_norm(x, /, *, keepdims=False, ord="fro"):
3420 """
3421 Computes the matrix norm of a matrix (or a stack of matrices) ``x``.
3423 This function is Array API compatible.
3425 Parameters
3426 ----------
3427 x : array_like
3428 Input array having shape (..., M, N) and whose two innermost
3429 dimensions form ``MxN`` matrices.
3430 keepdims : bool, optional
3431 If this is set to True, the axes which are normed over are left in
3432 the result as dimensions with size one. Default: False.
3433 ord : {1, -1, 2, -2, inf, -inf, 'fro', 'nuc'}, optional
3434 The order of the norm. For details see the table under ``Notes``
3435 in `numpy.linalg.norm`.
3437 See Also
3438 --------
3439 numpy.linalg.norm : Generic norm function
3441 Examples
3442 --------
3443 >>> from numpy import linalg as LA
3444 >>> a = np.arange(9) - 4
3445 >>> a
3446 array([-4, -3, -2, ..., 2, 3, 4])
3447 >>> b = a.reshape((3, 3))
3448 >>> b
3449 array([[-4, -3, -2],
3450 [-1, 0, 1],
3451 [ 2, 3, 4]])
3453 >>> LA.matrix_norm(b)
3454 7.745966692414834
3455 >>> LA.matrix_norm(b, ord='fro')
3456 7.745966692414834
3457 >>> LA.matrix_norm(b, ord=np.inf)
3458 9.0
3459 >>> LA.matrix_norm(b, ord=-np.inf)
3460 2.0
3462 >>> LA.matrix_norm(b, ord=1)
3463 7.0
3464 >>> LA.matrix_norm(b, ord=-1)
3465 6.0
3466 >>> LA.matrix_norm(b, ord=2)
3467 7.3484692283495345
3468 >>> LA.matrix_norm(b, ord=-2)
3469 1.8570331885190563e-016 # may vary
3471 """
3472 x = asanyarray(x)
3473 return norm(x, axis=(-2, -1), keepdims=keepdims, ord=ord)
3476# vector_norm
3478def _vector_norm_dispatcher(x, /, *, axis=None, keepdims=None, ord=None):
3479 return (x,)
3481@array_function_dispatch(_vector_norm_dispatcher)
3482def vector_norm(x, /, *, axis=None, keepdims=False, ord=2):
3483 """
3484 Computes the vector norm of a vector (or batch of vectors) ``x``.
3486 This function is Array API compatible.
3488 Parameters
3489 ----------
3490 x : array_like
3491 Input array.
3492 axis : {None, int, 2-tuple of ints}, optional
3493 If an integer, ``axis`` specifies the axis (dimension) along which
3494 to compute vector norms. If an n-tuple, ``axis`` specifies the axes
3495 (dimensions) along which to compute batched vector norms. If ``None``,
3496 the vector norm must be computed over all array values (i.e.,
3497 equivalent to computing the vector norm of a flattened array).
3498 Default: ``None``.
3499 keepdims : bool, optional
3500 If this is set to True, the axes which are normed over are left in
3501 the result as dimensions with size one. Default: False.
3502 ord : {int, float, inf, -inf}, optional
3503 The order of the norm. For details see the table under ``Notes``
3504 in `numpy.linalg.norm`.
3506 See Also
3507 --------
3508 numpy.linalg.norm : Generic norm function
3510 Examples
3511 --------
3512 >>> from numpy import linalg as LA
3513 >>> a = np.arange(9) + 1
3514 >>> a
3515 array([1, 2, 3, 4, 5, 6, 7, 8, 9])
3516 >>> b = a.reshape((3, 3))
3517 >>> b
3518 array([[1, 2, 3],
3519 [4, 5, 6],
3520 [7, 8, 9]])
3522 >>> LA.vector_norm(b)
3523 16.881943016134134
3524 >>> LA.vector_norm(b, ord=np.inf)
3525 9.0
3526 >>> LA.vector_norm(b, ord=-np.inf)
3527 1.0
3529 >>> LA.vector_norm(b, ord=0)
3530 9.0
3531 >>> LA.vector_norm(b, ord=1)
3532 45.0
3533 >>> LA.vector_norm(b, ord=-1)
3534 0.3534857623790153
3535 >>> LA.vector_norm(b, ord=2)
3536 16.881943016134134
3537 >>> LA.vector_norm(b, ord=-2)
3538 0.8058837395885292
3540 """
3541 x = asanyarray(x)
3542 shape = list(x.shape)
3543 if axis is None:
3544 # Note: np.linalg.norm() doesn't handle 0-D arrays
3545 x = x.ravel()
3546 _axis = 0
3547 elif isinstance(axis, tuple):
3548 # Note: The axis argument supports any number of axes, whereas
3549 # np.linalg.norm() only supports a single axis for vector norm.
3550 normalized_axis = normalize_axis_tuple(axis, x.ndim)
3551 rest = tuple(i for i in range(x.ndim) if i not in normalized_axis)
3552 newshape = axis + rest
3553 x = _core_transpose(x, newshape).reshape(
3554 (
3555 prod([x.shape[i] for i in axis], dtype=int),
3556 *[x.shape[i] for i in rest]
3557 )
3558 )
3559 _axis = 0
3560 else:
3561 _axis = axis
3563 res = norm(x, axis=_axis, ord=ord)
3565 if keepdims:
3566 # We can't reuse np.linalg.norm(keepdims) because of the reshape hacks
3567 # above to avoid matrix norm logic.
3568 _axis = normalize_axis_tuple(
3569 range(len(shape)) if axis is None else axis, len(shape)
3570 )
3571 for i in _axis:
3572 shape[i] = 1
3573 res = res.reshape(tuple(shape))
3575 return res
3578# vecdot
3580def _vecdot_dispatcher(x1, x2, /, *, axis=None):
3581 return (x1, x2)
3583@array_function_dispatch(_vecdot_dispatcher)
3584def vecdot(x1, x2, /, *, axis=-1):
3585 """
3586 Computes the vector dot product.
3588 This function is restricted to arguments compatible with the Array API,
3589 contrary to :func:`numpy.vecdot`.
3591 Let :math:`\\mathbf{a}` be a vector in ``x1`` and :math:`\\mathbf{b}` be
3592 a corresponding vector in ``x2``. The dot product is defined as:
3594 .. math::
3595 \\mathbf{a} \\cdot \\mathbf{b} = \\sum_{i=0}^{n-1} \\overline{a_i}b_i
3597 over the dimension specified by ``axis`` and where :math:`\\overline{a_i}`
3598 denotes the complex conjugate if :math:`a_i` is complex and the identity
3599 otherwise.
3601 Parameters
3602 ----------
3603 x1 : array_like
3604 First input array.
3605 x2 : array_like
3606 Second input array.
3607 axis : int, optional
3608 Axis over which to compute the dot product. Default: ``-1``.
3610 Returns
3611 -------
3612 output : ndarray
3613 The vector dot product of the input.
3615 See Also
3616 --------
3617 numpy.vecdot
3619 Examples
3620 --------
3621 Get the projected size along a given normal for an array of vectors.
3623 >>> v = np.array([[0., 5., 0.], [0., 0., 10.], [0., 6., 8.]])
3624 >>> n = np.array([0., 0.6, 0.8])
3625 >>> np.linalg.vecdot(v, n)
3626 array([ 3., 8., 10.])
3628 """
3629 return _core_vecdot(x1, x2, axis=axis)