Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.11/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 Any, NamedTuple
24from numpy._core import (
25 abs,
26 add,
27 all,
28 amax,
29 amin,
30 argsort,
31 array,
32 asanyarray,
33 asarray,
34 atleast_2d,
35 cdouble,
36 complexfloating,
37 count_nonzero,
38 cross as _core_cross,
39 csingle,
40 diagonal as _core_diagonal,
41 divide,
42 dot,
43 double,
44 empty,
45 empty_like,
46 errstate,
47 finfo,
48 inexact,
49 inf,
50 intc,
51 intp,
52 isfinite,
53 isnan,
54 matmul as _core_matmul,
55 matrix_transpose as _core_matrix_transpose,
56 moveaxis,
57 multiply,
58 newaxis,
59 object_,
60 outer as _core_outer,
61 overrides,
62 prod,
63 reciprocal,
64 single,
65 sort,
66 sqrt,
67 sum,
68 swapaxes,
69 tensordot as _core_tensordot,
70 trace as _core_trace,
71 transpose as _core_transpose,
72 vecdot as _core_vecdot,
73 zeros,
74)
75from numpy._globals import _NoValue
76from numpy._typing import NDArray
77from numpy._utils import set_module
78from numpy.lib._twodim_base_impl import eye, triu
79from numpy.lib.array_utils import normalize_axis_index, normalize_axis_tuple
80from numpy.linalg import _umath_linalg
83class EigResult(NamedTuple):
84 eigenvalues: NDArray[Any]
85 eigenvectors: NDArray[Any]
87class EighResult(NamedTuple):
88 eigenvalues: NDArray[Any]
89 eigenvectors: NDArray[Any]
91class QRResult(NamedTuple):
92 Q: NDArray[Any]
93 R: NDArray[Any]
95class SlogdetResult(NamedTuple):
96 sign: NDArray[Any]
97 logabsdet: NDArray[Any]
99class SVDResult(NamedTuple):
100 U: NDArray[Any]
101 S: NDArray[Any]
102 Vh: NDArray[Any]
105array_function_dispatch = functools.partial(
106 overrides.array_function_dispatch, module='numpy.linalg'
107)
110fortran_int = intc
113@set_module('numpy.linalg')
114class LinAlgError(ValueError):
115 """
116 Generic Python-exception-derived object raised by linalg functions.
118 General purpose exception class, derived from Python's ValueError
119 class, programmatically raised in linalg functions when a Linear
120 Algebra-related condition would prevent further correct execution of the
121 function.
123 Parameters
124 ----------
125 None
127 Examples
128 --------
129 >>> from numpy import linalg as LA
130 >>> LA.inv(np.zeros((2,2)))
131 Traceback (most recent call last):
132 File "<stdin>", line 1, in <module>
133 File "...linalg.py", line 350,
134 in inv return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))
135 File "...linalg.py", line 249,
136 in solve
137 raise LinAlgError('Singular matrix')
138 numpy.linalg.LinAlgError: Singular matrix
140 """
143def _raise_linalgerror_singular(err, flag):
144 raise LinAlgError("Singular matrix")
146def _raise_linalgerror_nonposdef(err, flag):
147 raise LinAlgError("Matrix is not positive definite")
149def _raise_linalgerror_eigenvalues_nonconvergence(err, flag):
150 raise LinAlgError("Eigenvalues did not converge")
152def _raise_linalgerror_svd_nonconvergence(err, flag):
153 raise LinAlgError("SVD did not converge")
155def _raise_linalgerror_lstsq(err, flag):
156 raise LinAlgError("SVD did not converge in Linear Least Squares")
158def _raise_linalgerror_qr(err, flag):
159 raise LinAlgError("Incorrect argument found while performing "
160 "QR factorization")
163def _makearray(a):
164 new = asarray(a)
165 wrap = getattr(a, "__array_wrap__", new.__array_wrap__)
166 return new, wrap
168def isComplexType(t):
169 return issubclass(t, complexfloating)
172_real_types_map = {single: single,
173 double: double,
174 csingle: single,
175 cdouble: double}
177_complex_types_map = {single: csingle,
178 double: cdouble,
179 csingle: csingle,
180 cdouble: cdouble}
182def _realType(t, default=double):
183 return _real_types_map.get(t, default)
185def _complexType(t, default=cdouble):
186 return _complex_types_map.get(t, default)
188def _commonType(*arrays):
189 # in lite version, use higher precision (always double or cdouble)
190 result_type = single
191 is_complex = False
192 for a in arrays:
193 type_ = a.dtype.type
194 if issubclass(type_, inexact):
195 if isComplexType(type_):
196 is_complex = True
197 rt = _realType(type_, default=None)
198 if rt is double:
199 result_type = double
200 elif rt is None:
201 # unsupported inexact scalar
202 raise TypeError(f"array type {a.dtype.name} is unsupported in linalg")
203 else:
204 result_type = double
205 if is_complex:
206 result_type = _complex_types_map[result_type]
207 return cdouble, result_type
208 else:
209 return double, result_type
212def _to_native_byte_order(*arrays):
213 ret = []
214 for arr in arrays:
215 if arr.dtype.byteorder not in ('=', '|'):
216 ret.append(asarray(arr, dtype=arr.dtype.newbyteorder('=')))
217 else:
218 ret.append(arr)
219 if len(ret) == 1:
220 return ret[0]
221 else:
222 return ret
225def _assert_2d(*arrays):
226 for a in arrays:
227 if a.ndim != 2:
228 raise LinAlgError('%d-dimensional array given. Array must be '
229 'two-dimensional' % a.ndim)
231def _assert_stacked_2d(*arrays):
232 for a in arrays:
233 if a.ndim < 2:
234 raise LinAlgError('%d-dimensional array given. Array must be '
235 'at least two-dimensional' % a.ndim)
237def _assert_stacked_square(*arrays):
238 for a in arrays:
239 try:
240 m, n = a.shape[-2:]
241 except ValueError:
242 raise LinAlgError('%d-dimensional array given. Array must be '
243 'at least two-dimensional' % a.ndim)
244 if m != n:
245 raise LinAlgError('Last 2 dimensions of the array must be square')
247def _assert_finite(*arrays):
248 for a in arrays:
249 if not isfinite(a).all():
250 raise LinAlgError("Array must not contain infs or NaNs")
252def _is_empty_2d(arr):
253 # check size first for efficiency
254 return arr.size == 0 and prod(arr.shape[-2:]) == 0
257def transpose(a):
258 """
259 Transpose each matrix in a stack of matrices.
261 Unlike np.transpose, this only swaps the last two axes, rather than all of
262 them
264 Parameters
265 ----------
266 a : (...,M,N) array_like
268 Returns
269 -------
270 aT : (...,N,M) ndarray
271 """
272 return swapaxes(a, -1, -2)
274# Linear equations
276def _tensorsolve_dispatcher(a, b, axes=None):
277 return (a, b)
280@array_function_dispatch(_tensorsolve_dispatcher)
281def tensorsolve(a, b, axes=None):
282 """
283 Solve the tensor equation ``a x = b`` for x.
285 It is assumed that all indices of `x` are summed over in the product,
286 together with the rightmost indices of `a`, as is done in, for example,
287 ``tensordot(a, x, axes=x.ndim)``.
289 Parameters
290 ----------
291 a : array_like
292 Coefficient tensor, of shape ``b.shape + Q``. `Q`, a tuple, equals
293 the shape of that sub-tensor of `a` consisting of the appropriate
294 number of its rightmost indices, and must be such that
295 ``prod(Q) == prod(b.shape)`` (in which sense `a` is said to be
296 'square').
297 b : array_like
298 Right-hand tensor, which can be of any shape.
299 axes : tuple of ints, optional
300 Axes in `a` to reorder to the right, before inversion.
301 If None (default), no reordering is done.
303 Returns
304 -------
305 x : ndarray, shape Q
307 Raises
308 ------
309 LinAlgError
310 If `a` is singular or not 'square' (in the above sense).
312 See Also
313 --------
314 numpy.tensordot, tensorinv, numpy.einsum
316 Examples
317 --------
318 >>> import numpy as np
319 >>> a = np.eye(2*3*4).reshape((2*3, 4, 2, 3, 4))
320 >>> rng = np.random.default_rng()
321 >>> b = rng.normal(size=(2*3, 4))
322 >>> x = np.linalg.tensorsolve(a, b)
323 >>> x.shape
324 (2, 3, 4)
325 >>> np.allclose(np.tensordot(a, x, axes=3), b)
326 True
328 """
329 a, wrap = _makearray(a)
330 b = asarray(b)
331 an = a.ndim
333 if axes is not None:
334 allaxes = list(range(an))
335 for k in axes:
336 allaxes.remove(k)
337 allaxes.insert(an, k)
338 a = a.transpose(allaxes)
340 oldshape = a.shape[-(an - b.ndim):]
341 prod = 1
342 for k in oldshape:
343 prod *= k
345 if a.size != prod ** 2:
346 raise LinAlgError(
347 "Input arrays must satisfy the requirement \
348 prod(a.shape[b.ndim:]) == prod(a.shape[:b.ndim])"
349 )
351 a = a.reshape(prod, prod)
352 b = b.ravel()
353 res = wrap(solve(a, b))
354 res.shape = oldshape
355 return res
358def _solve_dispatcher(a, b):
359 return (a, b)
362@array_function_dispatch(_solve_dispatcher)
363def solve(a, b):
364 """
365 Solve a linear matrix equation, or system of linear scalar equations.
367 Computes the "exact" solution, `x`, of the well-determined, i.e., full
368 rank, linear matrix equation `ax = b`.
370 Parameters
371 ----------
372 a : (..., M, M) array_like
373 Coefficient matrix.
374 b : {(M,), (..., M, K)}, array_like
375 Ordinate or "dependent variable" values.
377 Returns
378 -------
379 x : {(..., M,), (..., M, K)} ndarray
380 Solution to the system a x = b. Returned shape is (..., M) if b is
381 shape (M,) and (..., M, K) if b is (..., M, K), where the "..." part is
382 broadcasted between a and b.
384 Raises
385 ------
386 LinAlgError
387 If `a` is singular or not square.
389 See Also
390 --------
391 scipy.linalg.solve : Similar function in SciPy.
393 Notes
394 -----
395 Broadcasting rules apply, see the `numpy.linalg` documentation for
396 details.
398 The solutions are computed using LAPACK routine ``_gesv``.
400 `a` must be square and of full-rank, i.e., all rows (or, equivalently,
401 columns) must be linearly independent; if either is not true, use
402 `lstsq` for the least-squares best "solution" of the
403 system/equation.
405 .. versionchanged:: 2.0
407 The b array is only treated as a shape (M,) column vector if it is
408 exactly 1-dimensional. In all other instances it is treated as a stack
409 of (M, K) matrices. Previously b would be treated as a stack of (M,)
410 vectors if b.ndim was equal to a.ndim - 1.
412 References
413 ----------
414 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
415 FL, Academic Press, Inc., 1980, pg. 22.
417 Examples
418 --------
419 Solve the system of equations:
420 ``x0 + 2 * x1 = 1`` and
421 ``3 * x0 + 5 * x1 = 2``:
423 >>> import numpy as np
424 >>> a = np.array([[1, 2], [3, 5]])
425 >>> b = np.array([1, 2])
426 >>> x = np.linalg.solve(a, b)
427 >>> x
428 array([-1., 1.])
430 Check that the solution is correct:
432 >>> np.allclose(np.dot(a, x), b)
433 True
435 """
436 a, _ = _makearray(a)
437 _assert_stacked_square(a)
438 b, wrap = _makearray(b)
439 t, result_t = _commonType(a, b)
441 # We use the b = (..., M,) logic, only if the number of extra dimensions
442 # match exactly
443 if b.ndim == 1:
444 gufunc = _umath_linalg.solve1
445 else:
446 gufunc = _umath_linalg.solve
448 signature = 'DD->D' if isComplexType(t) else 'dd->d'
449 with errstate(call=_raise_linalgerror_singular, invalid='call',
450 over='ignore', divide='ignore', under='ignore'):
451 r = gufunc(a, b, signature=signature)
453 return wrap(r.astype(result_t, copy=False))
456def _tensorinv_dispatcher(a, ind=None):
457 return (a,)
460@array_function_dispatch(_tensorinv_dispatcher)
461def tensorinv(a, ind=2):
462 """
463 Compute the 'inverse' of an N-dimensional array.
465 The result is an inverse for `a` relative to the tensordot operation
466 ``tensordot(a, b, ind)``, i. e., up to floating-point accuracy,
467 ``tensordot(tensorinv(a), a, ind)`` is the "identity" tensor for the
468 tensordot operation.
470 Parameters
471 ----------
472 a : array_like
473 Tensor to 'invert'. Its shape must be 'square', i. e.,
474 ``prod(a.shape[:ind]) == prod(a.shape[ind:])``.
475 ind : int, optional
476 Number of first indices that are involved in the inverse sum.
477 Must be a positive integer, default is 2.
479 Returns
480 -------
481 b : ndarray
482 `a`'s tensordot inverse, shape ``a.shape[ind:] + a.shape[:ind]``.
484 Raises
485 ------
486 LinAlgError
487 If `a` is singular or not 'square' (in the above sense).
489 See Also
490 --------
491 numpy.tensordot, tensorsolve
493 Examples
494 --------
495 >>> import numpy as np
496 >>> a = np.eye(4*6).reshape((4, 6, 8, 3))
497 >>> ainv = np.linalg.tensorinv(a, ind=2)
498 >>> ainv.shape
499 (8, 3, 4, 6)
500 >>> rng = np.random.default_rng()
501 >>> b = rng.normal(size=(4, 6))
502 >>> np.allclose(np.tensordot(ainv, b), np.linalg.tensorsolve(a, b))
503 True
505 >>> a = np.eye(4*6).reshape((24, 8, 3))
506 >>> ainv = np.linalg.tensorinv(a, ind=1)
507 >>> ainv.shape
508 (8, 3, 24)
509 >>> rng = np.random.default_rng()
510 >>> b = rng.normal(size=24)
511 >>> np.allclose(np.tensordot(ainv, b, 1), np.linalg.tensorsolve(a, b))
512 True
514 """
515 a = asarray(a)
516 oldshape = a.shape
517 prod = 1
518 if ind > 0:
519 invshape = oldshape[ind:] + oldshape[:ind]
520 for k in oldshape[ind:]:
521 prod *= k
522 else:
523 raise ValueError("Invalid ind argument.")
524 a = a.reshape(prod, -1)
525 ia = inv(a)
526 return ia.reshape(*invshape)
529# Matrix inversion
531def _unary_dispatcher(a):
532 return (a,)
535@array_function_dispatch(_unary_dispatcher)
536def inv(a):
537 """
538 Compute the inverse of a matrix.
540 Given a square matrix `a`, return the matrix `ainv` satisfying
541 ``a @ ainv = ainv @ a = eye(a.shape[0])``.
543 Parameters
544 ----------
545 a : (..., M, M) array_like
546 Matrix to be inverted.
548 Returns
549 -------
550 ainv : (..., M, M) ndarray or matrix
551 Inverse of the matrix `a`.
553 Raises
554 ------
555 LinAlgError
556 If `a` is not square or inversion fails.
558 See Also
559 --------
560 scipy.linalg.inv : Similar function in SciPy.
561 numpy.linalg.cond : Compute the condition number of a matrix.
562 numpy.linalg.svd : Compute the singular value decomposition of a matrix.
564 Notes
565 -----
566 Broadcasting rules apply, see the `numpy.linalg` documentation for
567 details.
569 If `a` is detected to be singular, a `LinAlgError` is raised. If `a` is
570 ill-conditioned, a `LinAlgError` may or may not be raised, and results may
571 be inaccurate due to floating-point errors.
573 References
574 ----------
575 .. [1] Wikipedia, "Condition number",
576 https://en.wikipedia.org/wiki/Condition_number
578 Examples
579 --------
580 >>> import numpy as np
581 >>> from numpy.linalg import inv
582 >>> a = np.array([[1., 2.], [3., 4.]])
583 >>> ainv = inv(a)
584 >>> np.allclose(a @ ainv, np.eye(2))
585 True
586 >>> np.allclose(ainv @ a, np.eye(2))
587 True
589 If a is a matrix object, then the return value is a matrix as well:
591 >>> ainv = inv(np.matrix(a))
592 >>> ainv
593 matrix([[-2. , 1. ],
594 [ 1.5, -0.5]])
596 Inverses of several matrices can be computed at once:
598 >>> a = np.array([[[1., 2.], [3., 4.]], [[1, 3], [3, 5]]])
599 >>> inv(a)
600 array([[[-2. , 1. ],
601 [ 1.5 , -0.5 ]],
602 [[-1.25, 0.75],
603 [ 0.75, -0.25]]])
605 If a matrix is close to singular, the computed inverse may not satisfy
606 ``a @ ainv = ainv @ a = eye(a.shape[0])`` even if a `LinAlgError`
607 is not raised:
609 >>> a = np.array([[2,4,6],[2,0,2],[6,8,14]])
610 >>> inv(a) # No errors raised
611 array([[-1.12589991e+15, -5.62949953e+14, 5.62949953e+14],
612 [-1.12589991e+15, -5.62949953e+14, 5.62949953e+14],
613 [ 1.12589991e+15, 5.62949953e+14, -5.62949953e+14]])
614 >>> a @ inv(a)
615 array([[ 0. , -0.5 , 0. ], # may vary
616 [-0.5 , 0.625, 0.25 ],
617 [ 0. , 0. , 1. ]])
619 To detect ill-conditioned matrices, you can use `numpy.linalg.cond` to
620 compute its *condition number* [1]_. The larger the condition number, the
621 more ill-conditioned the matrix is. As a rule of thumb, if the condition
622 number ``cond(a) = 10**k``, then you may lose up to ``k`` digits of
623 accuracy on top of what would be lost to the numerical method due to loss
624 of precision from arithmetic methods.
626 >>> from numpy.linalg import cond
627 >>> cond(a)
628 np.float64(8.659885634118668e+17) # may vary
630 It is also possible to detect ill-conditioning by inspecting the matrix's
631 singular values directly. The ratio between the largest and the smallest
632 singular value is the condition number:
634 >>> from numpy.linalg import svd
635 >>> sigma = svd(a, compute_uv=False) # Do not compute singular vectors
636 >>> sigma.max()/sigma.min()
637 8.659885634118668e+17 # may vary
639 """
640 a, wrap = _makearray(a)
641 _assert_stacked_square(a)
642 t, result_t = _commonType(a)
644 signature = 'D->D' if isComplexType(t) else 'd->d'
645 with errstate(call=_raise_linalgerror_singular, invalid='call',
646 over='ignore', divide='ignore', under='ignore'):
647 ainv = _umath_linalg.inv(a, signature=signature)
648 return wrap(ainv.astype(result_t, copy=False))
651def _matrix_power_dispatcher(a, n):
652 return (a,)
655@array_function_dispatch(_matrix_power_dispatcher)
656def matrix_power(a, n):
657 """
658 Raise a square matrix to the (integer) power `n`.
660 For positive integers `n`, the power is computed by repeated matrix
661 squarings and matrix multiplications. If ``n == 0``, the identity matrix
662 of the same shape as M is returned. If ``n < 0``, the inverse
663 is computed and then raised to the ``abs(n)``.
665 .. note:: Stacks of object matrices are not currently supported.
667 Parameters
668 ----------
669 a : (..., M, M) array_like
670 Matrix to be "powered".
671 n : int
672 The exponent can be any integer or long integer, positive,
673 negative, or zero.
675 Returns
676 -------
677 a**n : (..., M, M) ndarray or matrix object
678 The return value is the same shape and type as `M`;
679 if the exponent is positive or zero then the type of the
680 elements is the same as those of `M`. If the exponent is
681 negative the elements are floating-point.
683 Raises
684 ------
685 LinAlgError
686 For matrices that are not square or that (for negative powers) cannot
687 be inverted numerically.
689 Examples
690 --------
691 >>> import numpy as np
692 >>> from numpy.linalg import matrix_power
693 >>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit
694 >>> matrix_power(i, 3) # should = -i
695 array([[ 0, -1],
696 [ 1, 0]])
697 >>> matrix_power(i, 0)
698 array([[1, 0],
699 [0, 1]])
700 >>> matrix_power(i, -3) # should = 1/(-i) = i, but w/ f.p. elements
701 array([[ 0., 1.],
702 [-1., 0.]])
704 Somewhat more sophisticated example
706 >>> q = np.zeros((4, 4))
707 >>> q[0:2, 0:2] = -i
708 >>> q[2:4, 2:4] = i
709 >>> q # one of the three quaternion units not equal to 1
710 array([[ 0., -1., 0., 0.],
711 [ 1., 0., 0., 0.],
712 [ 0., 0., 0., 1.],
713 [ 0., 0., -1., 0.]])
714 >>> matrix_power(q, 2) # = -np.eye(4)
715 array([[-1., 0., 0., 0.],
716 [ 0., -1., 0., 0.],
717 [ 0., 0., -1., 0.],
718 [ 0., 0., 0., -1.]])
720 """
721 a = asanyarray(a)
722 _assert_stacked_square(a)
724 try:
725 n = operator.index(n)
726 except TypeError as e:
727 raise TypeError("exponent must be an integer") from e
729 # Fall back on dot for object arrays. Object arrays are not supported by
730 # the current implementation of matmul using einsum
731 if a.dtype != object:
732 fmatmul = matmul
733 elif a.ndim == 2:
734 fmatmul = dot
735 else:
736 raise NotImplementedError(
737 "matrix_power not supported for stacks of object arrays")
739 if n == 0:
740 a = empty_like(a)
741 a[...] = eye(a.shape[-2], dtype=a.dtype)
742 return a
744 elif n < 0:
745 a = inv(a)
746 n = abs(n)
748 # short-cuts.
749 if n == 1:
750 return a
752 elif n == 2:
753 return fmatmul(a, a)
755 elif n == 3:
756 return fmatmul(fmatmul(a, a), a)
758 # Use binary decomposition to reduce the number of matrix multiplications.
759 # Here, we iterate over the bits of n, from LSB to MSB, raise `a` to
760 # increasing powers of 2, and multiply into the result as needed.
761 z = result = None
762 while n > 0:
763 z = a if z is None else fmatmul(z, z)
764 n, bit = divmod(n, 2)
765 if bit:
766 result = z if result is None else fmatmul(result, z)
768 return result
771# Cholesky decomposition
773def _cholesky_dispatcher(a, /, *, upper=None):
774 return (a,)
777@array_function_dispatch(_cholesky_dispatcher)
778def cholesky(a, /, *, upper=False):
779 """
780 Cholesky decomposition.
782 Return the lower or upper Cholesky decomposition, ``L * L.H`` or
783 ``U.H * U``, of the square matrix ``a``, where ``L`` is lower-triangular,
784 ``U`` is upper-triangular, and ``.H`` is the conjugate transpose operator
785 (which is the ordinary transpose if ``a`` is real-valued). ``a`` must be
786 Hermitian (symmetric if real-valued) and positive-definite. No checking is
787 performed to verify whether ``a`` is Hermitian or not. In addition, only
788 the lower or upper-triangular and diagonal elements of ``a`` are used.
789 Only ``L`` or ``U`` is actually returned.
791 Parameters
792 ----------
793 a : (..., M, M) array_like
794 Hermitian (symmetric if all elements are real), positive-definite
795 input matrix.
796 upper : bool
797 If ``True``, the result must be the upper-triangular Cholesky factor.
798 If ``False``, the result must be the lower-triangular Cholesky factor.
799 Default: ``False``.
801 Returns
802 -------
803 L : (..., M, M) array_like
804 Lower or upper-triangular Cholesky factor of `a`. Returns a matrix
805 object if `a` is a matrix object.
807 Raises
808 ------
809 LinAlgError
810 If the decomposition fails, for example, if `a` is not
811 positive-definite.
813 See Also
814 --------
815 scipy.linalg.cholesky : Similar function in SciPy.
816 scipy.linalg.cholesky_banded : Cholesky decompose a banded Hermitian
817 positive-definite matrix.
818 scipy.linalg.cho_factor : Cholesky decomposition of a matrix, to use in
819 `scipy.linalg.cho_solve`.
821 Notes
822 -----
823 Broadcasting rules apply, see the `numpy.linalg` documentation for
824 details.
826 The Cholesky decomposition is often used as a fast way of solving
828 .. math:: A \\mathbf{x} = \\mathbf{b}
830 (when `A` is both Hermitian/symmetric and positive-definite).
832 First, we solve for :math:`\\mathbf{y}` in
834 .. math:: L \\mathbf{y} = \\mathbf{b},
836 and then for :math:`\\mathbf{x}` in
838 .. math:: L^{H} \\mathbf{x} = \\mathbf{y}.
840 Examples
841 --------
842 >>> import numpy as np
843 >>> A = np.array([[1,-2j],[2j,5]])
844 >>> A
845 array([[ 1.+0.j, -0.-2.j],
846 [ 0.+2.j, 5.+0.j]])
847 >>> L = np.linalg.cholesky(A)
848 >>> L
849 array([[1.+0.j, 0.+0.j],
850 [0.+2.j, 1.+0.j]])
851 >>> np.dot(L, L.T.conj()) # verify that L * L.H = A
852 array([[1.+0.j, 0.-2.j],
853 [0.+2.j, 5.+0.j]])
854 >>> A = [[1,-2j],[2j,5]] # what happens if A is only array_like?
855 >>> np.linalg.cholesky(A) # an ndarray object is returned
856 array([[1.+0.j, 0.+0.j],
857 [0.+2.j, 1.+0.j]])
858 >>> # But a matrix object is returned if A is a matrix object
859 >>> np.linalg.cholesky(np.matrix(A))
860 matrix([[ 1.+0.j, 0.+0.j],
861 [ 0.+2.j, 1.+0.j]])
862 >>> # The upper-triangular Cholesky factor can also be obtained.
863 >>> np.linalg.cholesky(A, upper=True)
864 array([[1.-0.j, 0.-2.j],
865 [0.-0.j, 1.-0.j]])
867 """
868 gufunc = _umath_linalg.cholesky_up if upper else _umath_linalg.cholesky_lo
869 a, wrap = _makearray(a)
870 _assert_stacked_square(a)
871 t, result_t = _commonType(a)
872 signature = 'D->D' if isComplexType(t) else 'd->d'
873 with errstate(call=_raise_linalgerror_nonposdef, invalid='call',
874 over='ignore', divide='ignore', under='ignore'):
875 r = gufunc(a, signature=signature)
876 return wrap(r.astype(result_t, copy=False))
879# outer product
882def _outer_dispatcher(x1, x2):
883 return (x1, x2)
886@array_function_dispatch(_outer_dispatcher)
887def outer(x1, x2, /):
888 """
889 Compute the outer product of two vectors.
891 This function is Array API compatible. Compared to ``np.outer``
892 it accepts 1-dimensional inputs only.
894 Parameters
895 ----------
896 x1 : (M,) array_like
897 One-dimensional input array of size ``N``.
898 Must have a numeric data type.
899 x2 : (N,) array_like
900 One-dimensional input array of size ``M``.
901 Must have a numeric data type.
903 Returns
904 -------
905 out : (M, N) ndarray
906 ``out[i, j] = a[i] * b[j]``
908 See also
909 --------
910 outer
912 Examples
913 --------
914 Make a (*very* coarse) grid for computing a Mandelbrot set:
916 >>> rl = np.linalg.outer(np.ones((5,)), np.linspace(-2, 2, 5))
917 >>> rl
918 array([[-2., -1., 0., 1., 2.],
919 [-2., -1., 0., 1., 2.],
920 [-2., -1., 0., 1., 2.],
921 [-2., -1., 0., 1., 2.],
922 [-2., -1., 0., 1., 2.]])
923 >>> im = np.linalg.outer(1j*np.linspace(2, -2, 5), np.ones((5,)))
924 >>> im
925 array([[0.+2.j, 0.+2.j, 0.+2.j, 0.+2.j, 0.+2.j],
926 [0.+1.j, 0.+1.j, 0.+1.j, 0.+1.j, 0.+1.j],
927 [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
928 [0.-1.j, 0.-1.j, 0.-1.j, 0.-1.j, 0.-1.j],
929 [0.-2.j, 0.-2.j, 0.-2.j, 0.-2.j, 0.-2.j]])
930 >>> grid = rl + im
931 >>> grid
932 array([[-2.+2.j, -1.+2.j, 0.+2.j, 1.+2.j, 2.+2.j],
933 [-2.+1.j, -1.+1.j, 0.+1.j, 1.+1.j, 2.+1.j],
934 [-2.+0.j, -1.+0.j, 0.+0.j, 1.+0.j, 2.+0.j],
935 [-2.-1.j, -1.-1.j, 0.-1.j, 1.-1.j, 2.-1.j],
936 [-2.-2.j, -1.-2.j, 0.-2.j, 1.-2.j, 2.-2.j]])
938 An example using a "vector" of letters:
940 >>> x = np.array(['a', 'b', 'c'], dtype=object)
941 >>> np.linalg.outer(x, [1, 2, 3])
942 array([['a', 'aa', 'aaa'],
943 ['b', 'bb', 'bbb'],
944 ['c', 'cc', 'ccc']], dtype=object)
946 """
947 x1 = asanyarray(x1)
948 x2 = asanyarray(x2)
949 if x1.ndim != 1 or x2.ndim != 1:
950 raise ValueError(
951 "Input arrays must be one-dimensional, but they are "
952 f"{x1.ndim=} and {x2.ndim=}."
953 )
954 return _core_outer(x1, x2, out=None)
957# QR decomposition
960def _qr_dispatcher(a, mode=None):
961 return (a,)
964@array_function_dispatch(_qr_dispatcher)
965def qr(a, mode='reduced'):
966 """
967 Compute the qr factorization of a matrix.
969 Factor the matrix `a` as *qr*, where `q` is orthonormal and `r` is
970 upper-triangular.
972 Parameters
973 ----------
974 a : array_like, shape (..., M, N)
975 An array-like object with the dimensionality of at least 2.
976 mode : {'reduced', 'complete', 'r', 'raw'}, optional, default: 'reduced'
977 If K = min(M, N), then
979 * 'reduced' : returns Q, R with dimensions (..., M, K), (..., K, N)
980 * 'complete' : returns Q, R with dimensions (..., M, M), (..., M, N)
981 * 'r' : returns R only with dimensions (..., K, N)
982 * 'raw' : returns h, tau with dimensions (..., N, M), (..., K,)
984 The options 'reduced', 'complete, and 'raw' are new in numpy 1.8,
985 see the notes for more information. The default is 'reduced', and to
986 maintain backward compatibility with earlier versions of numpy both
987 it and the old default 'full' can be omitted. Note that array h
988 returned in 'raw' mode is transposed for calling Fortran. The
989 'economic' mode is deprecated. The modes 'full' and 'economic' may
990 be passed using only the first letter for backwards compatibility,
991 but all others must be spelled out. See the Notes for more
992 explanation.
995 Returns
996 -------
997 Q : ndarray of float or complex, optional
998 A matrix with orthonormal columns. When mode = 'complete' the
999 result is an orthogonal/unitary matrix depending on whether or not
1000 a is real/complex. The determinant may be either +/- 1 in that
1001 case. In case the number of dimensions in the input array is
1002 greater than 2 then a stack of the matrices with above properties
1003 is returned.
1004 R : ndarray of float or complex, optional
1005 The upper-triangular matrix or a stack of upper-triangular
1006 matrices if the number of dimensions in the input array is greater
1007 than 2.
1008 (h, tau) : ndarrays of np.double or np.cdouble, optional
1009 The array h contains the Householder reflectors that generate q
1010 along with r. The tau array contains scaling factors for the
1011 reflectors. In the deprecated 'economic' mode only h is returned.
1013 Raises
1014 ------
1015 LinAlgError
1016 If factoring fails.
1018 See Also
1019 --------
1020 scipy.linalg.qr : Similar function in SciPy.
1021 scipy.linalg.rq : Compute RQ decomposition of a matrix.
1023 Notes
1024 -----
1025 When mode is 'reduced' or 'complete', the result will be a namedtuple with
1026 the attributes ``Q`` and ``R``.
1028 This is an interface to the LAPACK routines ``dgeqrf``, ``zgeqrf``,
1029 ``dorgqr``, and ``zungqr``.
1031 For more information on the qr factorization, see for example:
1032 https://en.wikipedia.org/wiki/QR_factorization
1034 Subclasses of `ndarray` are preserved except for the 'raw' mode. So if
1035 `a` is of type `matrix`, all the return values will be matrices too.
1037 New 'reduced', 'complete', and 'raw' options for mode were added in
1038 NumPy 1.8.0 and the old option 'full' was made an alias of 'reduced'. In
1039 addition the options 'full' and 'economic' were deprecated. Because
1040 'full' was the previous default and 'reduced' is the new default,
1041 backward compatibility can be maintained by letting `mode` default.
1042 The 'raw' option was added so that LAPACK routines that can multiply
1043 arrays by q using the Householder reflectors can be used. Note that in
1044 this case the returned arrays are of type np.double or np.cdouble and
1045 the h array is transposed to be FORTRAN compatible. No routines using
1046 the 'raw' return are currently exposed by numpy, but some are available
1047 in lapack_lite and just await the necessary work.
1049 Examples
1050 --------
1051 >>> import numpy as np
1052 >>> rng = np.random.default_rng()
1053 >>> a = rng.normal(size=(9, 6))
1054 >>> Q, R = np.linalg.qr(a)
1055 >>> np.allclose(a, np.dot(Q, R)) # a does equal QR
1056 True
1057 >>> R2 = np.linalg.qr(a, mode='r')
1058 >>> np.allclose(R, R2) # mode='r' returns the same R as mode='full'
1059 True
1060 >>> a = np.random.normal(size=(3, 2, 2)) # Stack of 2 x 2 matrices as input
1061 >>> Q, R = np.linalg.qr(a)
1062 >>> Q.shape
1063 (3, 2, 2)
1064 >>> R.shape
1065 (3, 2, 2)
1066 >>> np.allclose(a, np.matmul(Q, R))
1067 True
1069 Example illustrating a common use of `qr`: solving of least squares
1070 problems
1072 What are the least-squares-best `m` and `y0` in ``y = y0 + mx`` for
1073 the following data: {(0,1), (1,0), (1,2), (2,1)}. (Graph the points
1074 and you'll see that it should be y0 = 0, m = 1.) The answer is provided
1075 by solving the over-determined matrix equation ``Ax = b``, where::
1077 A = array([[0, 1], [1, 1], [1, 1], [2, 1]])
1078 x = array([[y0], [m]])
1079 b = array([[1], [0], [2], [1]])
1081 If A = QR such that Q is orthonormal (which is always possible via
1082 Gram-Schmidt), then ``x = inv(R) * (Q.T) * b``. (In numpy practice,
1083 however, we simply use `lstsq`.)
1085 >>> A = np.array([[0, 1], [1, 1], [1, 1], [2, 1]])
1086 >>> A
1087 array([[0, 1],
1088 [1, 1],
1089 [1, 1],
1090 [2, 1]])
1091 >>> b = np.array([1, 2, 2, 3])
1092 >>> Q, R = np.linalg.qr(A)
1093 >>> p = np.dot(Q.T, b)
1094 >>> np.dot(np.linalg.inv(R), p)
1095 array([ 1., 1.])
1097 """
1098 if mode not in ('reduced', 'complete', 'r', 'raw'):
1099 if mode in ('f', 'full'):
1100 # 2013-04-01, 1.8
1101 msg = (
1102 "The 'full' option is deprecated in favor of 'reduced'.\n"
1103 "For backward compatibility let mode default."
1104 )
1105 warnings.warn(msg, DeprecationWarning, stacklevel=2)
1106 mode = 'reduced'
1107 elif mode in ('e', 'economic'):
1108 # 2013-04-01, 1.8
1109 msg = "The 'economic' option is deprecated."
1110 warnings.warn(msg, DeprecationWarning, stacklevel=2)
1111 mode = 'economic'
1112 else:
1113 raise ValueError(f"Unrecognized mode '{mode}'")
1115 a, wrap = _makearray(a)
1116 _assert_stacked_2d(a)
1117 m, n = a.shape[-2:]
1118 t, result_t = _commonType(a)
1119 a = a.astype(t, copy=True)
1120 a = _to_native_byte_order(a)
1121 mn = min(m, n)
1123 signature = 'D->D' if isComplexType(t) else 'd->d'
1124 with errstate(call=_raise_linalgerror_qr, invalid='call',
1125 over='ignore', divide='ignore', under='ignore'):
1126 tau = _umath_linalg.qr_r_raw(a, signature=signature)
1128 # handle modes that don't return q
1129 if mode == 'r':
1130 r = triu(a[..., :mn, :])
1131 r = r.astype(result_t, copy=False)
1132 return wrap(r)
1134 if mode == 'raw':
1135 q = transpose(a)
1136 q = q.astype(result_t, copy=False)
1137 tau = tau.astype(result_t, copy=False)
1138 return wrap(q), tau
1140 if mode == 'economic':
1141 a = a.astype(result_t, copy=False)
1142 return wrap(a)
1144 # mc is the number of columns in the resulting q
1145 # matrix. If the mode is complete then it is
1146 # same as number of rows, and if the mode is reduced,
1147 # then it is the minimum of number of rows and columns.
1148 if mode == 'complete' and m > n:
1149 mc = m
1150 gufunc = _umath_linalg.qr_complete
1151 else:
1152 mc = mn
1153 gufunc = _umath_linalg.qr_reduced
1155 signature = 'DD->D' if isComplexType(t) else 'dd->d'
1156 with errstate(call=_raise_linalgerror_qr, invalid='call',
1157 over='ignore', divide='ignore', under='ignore'):
1158 q = gufunc(a, tau, signature=signature)
1159 r = triu(a[..., :mc, :])
1161 q = q.astype(result_t, copy=False)
1162 r = r.astype(result_t, copy=False)
1164 return QRResult(wrap(q), wrap(r))
1166# Eigenvalues
1169@array_function_dispatch(_unary_dispatcher)
1170def eigvals(a):
1171 """
1172 Compute the eigenvalues of a general matrix.
1174 Main difference between `eigvals` and `eig`: the eigenvectors aren't
1175 returned.
1177 Parameters
1178 ----------
1179 a : (..., M, M) array_like
1180 A complex- or real-valued matrix whose eigenvalues will be computed.
1182 Returns
1183 -------
1184 w : (..., M,) ndarray
1185 The eigenvalues, each repeated according to its multiplicity.
1186 They are not necessarily ordered, nor are they necessarily
1187 real for real matrices.
1189 Raises
1190 ------
1191 LinAlgError
1192 If the eigenvalue computation does not converge.
1194 See Also
1195 --------
1196 eig : eigenvalues and right eigenvectors of general arrays
1197 eigvalsh : eigenvalues of real symmetric or complex Hermitian
1198 (conjugate symmetric) arrays.
1199 eigh : eigenvalues and eigenvectors of real symmetric or complex
1200 Hermitian (conjugate symmetric) arrays.
1201 scipy.linalg.eigvals : Similar function in SciPy.
1203 Notes
1204 -----
1205 Broadcasting rules apply, see the `numpy.linalg` documentation for
1206 details.
1208 This is implemented using the ``_geev`` LAPACK routines which compute
1209 the eigenvalues and eigenvectors of general square arrays.
1211 Examples
1212 --------
1213 Illustration, using the fact that the eigenvalues of a diagonal matrix
1214 are its diagonal elements, that multiplying a matrix on the left
1215 by an orthogonal matrix, `Q`, and on the right by `Q.T` (the transpose
1216 of `Q`), preserves the eigenvalues of the "middle" matrix. In other words,
1217 if `Q` is orthogonal, then ``Q * A * Q.T`` has the same eigenvalues as
1218 ``A``:
1220 >>> import numpy as np
1221 >>> from numpy import linalg as LA
1222 >>> x = np.random.random()
1223 >>> Q = np.array([[np.cos(x), -np.sin(x)], [np.sin(x), np.cos(x)]])
1224 >>> LA.norm(Q[0, :]), LA.norm(Q[1, :]), np.dot(Q[0, :],Q[1, :])
1225 (1.0, 1.0, 0.0)
1227 Now multiply a diagonal matrix by ``Q`` on one side and
1228 by ``Q.T`` on the other:
1230 >>> D = np.diag((-1,1))
1231 >>> LA.eigvals(D)
1232 array([-1., 1.])
1233 >>> A = np.dot(Q, D)
1234 >>> A = np.dot(A, Q.T)
1235 >>> LA.eigvals(A)
1236 array([ 1., -1.]) # random
1238 """
1239 a, wrap = _makearray(a)
1240 _assert_stacked_square(a)
1241 _assert_finite(a)
1242 t, result_t = _commonType(a)
1244 signature = 'D->D' if isComplexType(t) else 'd->D'
1245 with errstate(call=_raise_linalgerror_eigenvalues_nonconvergence,
1246 invalid='call', over='ignore', divide='ignore',
1247 under='ignore'):
1248 w = _umath_linalg.eigvals(a, signature=signature)
1250 if not isComplexType(t):
1251 if all(w.imag == 0):
1252 w = w.real
1253 result_t = _realType(result_t)
1254 else:
1255 result_t = _complexType(result_t)
1257 return w.astype(result_t, copy=False)
1260def _eigvalsh_dispatcher(a, UPLO=None):
1261 return (a,)
1264@array_function_dispatch(_eigvalsh_dispatcher)
1265def eigvalsh(a, UPLO='L'):
1266 """
1267 Compute the eigenvalues of a complex Hermitian or real symmetric matrix.
1269 Main difference from eigh: the eigenvectors are not computed.
1271 Parameters
1272 ----------
1273 a : (..., M, M) array_like
1274 A complex- or real-valued matrix whose eigenvalues are to be
1275 computed.
1276 UPLO : {'L', 'U'}, optional
1277 Specifies whether the calculation is done with the lower triangular
1278 part of `a` ('L', default) or the upper triangular part ('U').
1279 Irrespective of this value only the real parts of the diagonal will
1280 be considered in the computation to preserve the notion of a Hermitian
1281 matrix. It therefore follows that the imaginary part of the diagonal
1282 will always be treated as zero.
1284 Returns
1285 -------
1286 w : (..., M,) ndarray
1287 The eigenvalues in ascending order, each repeated according to
1288 its multiplicity.
1290 Raises
1291 ------
1292 LinAlgError
1293 If the eigenvalue computation does not converge.
1295 See Also
1296 --------
1297 eigh : eigenvalues and eigenvectors of real symmetric or complex Hermitian
1298 (conjugate symmetric) arrays.
1299 eigvals : eigenvalues of general real or complex arrays.
1300 eig : eigenvalues and right eigenvectors of general real or complex
1301 arrays.
1302 scipy.linalg.eigvalsh : Similar function in SciPy.
1304 Notes
1305 -----
1306 Broadcasting rules apply, see the `numpy.linalg` documentation for
1307 details.
1309 The eigenvalues are computed using LAPACK routines ``_syevd``, ``_heevd``.
1311 Examples
1312 --------
1313 >>> import numpy as np
1314 >>> from numpy import linalg as LA
1315 >>> a = np.array([[1, -2j], [2j, 5]])
1316 >>> LA.eigvalsh(a)
1317 array([ 0.17157288, 5.82842712]) # may vary
1319 >>> # demonstrate the treatment of the imaginary part of the diagonal
1320 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1321 >>> a
1322 array([[5.+2.j, 9.-2.j],
1323 [0.+2.j, 2.-1.j]])
1324 >>> # with UPLO='L' this is numerically equivalent to using LA.eigvals()
1325 >>> # with:
1326 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1327 >>> b
1328 array([[5.+0.j, 0.-2.j],
1329 [0.+2.j, 2.+0.j]])
1330 >>> wa = LA.eigvalsh(a)
1331 >>> wb = LA.eigvals(b)
1332 >>> wa
1333 array([1., 6.])
1334 >>> wb
1335 array([6.+0.j, 1.+0.j])
1337 """
1338 UPLO = UPLO.upper()
1339 if UPLO not in ('L', 'U'):
1340 raise ValueError("UPLO argument must be 'L' or 'U'")
1342 if UPLO == 'L':
1343 gufunc = _umath_linalg.eigvalsh_lo
1344 else:
1345 gufunc = _umath_linalg.eigvalsh_up
1347 a, wrap = _makearray(a)
1348 _assert_stacked_square(a)
1349 t, result_t = _commonType(a)
1350 signature = 'D->d' if isComplexType(t) else 'd->d'
1351 with errstate(call=_raise_linalgerror_eigenvalues_nonconvergence,
1352 invalid='call', over='ignore', divide='ignore',
1353 under='ignore'):
1354 w = gufunc(a, signature=signature)
1355 return w.astype(_realType(result_t), copy=False)
1358# Eigenvectors
1361@array_function_dispatch(_unary_dispatcher)
1362def eig(a):
1363 """
1364 Compute the eigenvalues and right eigenvectors of a square array.
1366 Parameters
1367 ----------
1368 a : (..., M, M) array
1369 Matrices for which the eigenvalues and right eigenvectors will
1370 be computed
1372 Returns
1373 -------
1374 A namedtuple with the following attributes:
1376 eigenvalues : (..., M) array
1377 The eigenvalues, each repeated according to its multiplicity.
1378 The eigenvalues are not necessarily ordered. The resulting
1379 array will be of complex type, unless the imaginary part is
1380 zero in which case it will be cast to a real type. When `a`
1381 is real the resulting eigenvalues will be real (0 imaginary
1382 part) or occur in conjugate pairs
1384 eigenvectors : (..., M, M) array
1385 The normalized (unit "length") eigenvectors, such that the
1386 column ``eigenvectors[:,i]`` is the eigenvector corresponding to the
1387 eigenvalue ``eigenvalues[i]``.
1389 Raises
1390 ------
1391 LinAlgError
1392 If the eigenvalue computation does not converge.
1394 See Also
1395 --------
1396 eigvals : eigenvalues of a non-symmetric array.
1397 eigh : eigenvalues and eigenvectors of a real symmetric or complex
1398 Hermitian (conjugate symmetric) array.
1399 eigvalsh : eigenvalues of a real symmetric or complex Hermitian
1400 (conjugate symmetric) array.
1401 scipy.linalg.eig : Similar function in SciPy that also solves the
1402 generalized eigenvalue problem.
1403 scipy.linalg.schur : Best choice for unitary and other non-Hermitian
1404 normal matrices.
1406 Notes
1407 -----
1408 Broadcasting rules apply, see the `numpy.linalg` documentation for
1409 details.
1411 This is implemented using the ``_geev`` LAPACK routines which compute
1412 the eigenvalues and eigenvectors of general square arrays.
1414 The number `w` is an eigenvalue of `a` if there exists a vector `v` such
1415 that ``a @ v = w * v``. Thus, the arrays `a`, `eigenvalues`, and
1416 `eigenvectors` satisfy the equations ``a @ eigenvectors[:,i] =
1417 eigenvalues[i] * eigenvectors[:,i]`` for :math:`i \\in \\{0,...,M-1\\}`.
1419 The array `eigenvectors` may not be of maximum rank, that is, some of the
1420 columns may be linearly dependent, although round-off error may obscure
1421 that fact. If the eigenvalues are all different, then theoretically the
1422 eigenvectors are linearly independent and `a` can be diagonalized by a
1423 similarity transformation using `eigenvectors`, i.e, ``inv(eigenvectors) @
1424 a @ eigenvectors`` is diagonal.
1426 For non-Hermitian normal matrices the SciPy function `scipy.linalg.schur`
1427 is preferred because the matrix `eigenvectors` is guaranteed to be
1428 unitary, which is not the case when using `eig`. The Schur factorization
1429 produces an upper triangular matrix rather than a diagonal matrix, but for
1430 normal matrices only the diagonal of the upper triangular matrix is
1431 needed, the rest is roundoff error.
1433 Finally, it is emphasized that `eigenvectors` consists of the *right* (as
1434 in right-hand side) eigenvectors of `a`. A vector `y` satisfying ``y.T @ a
1435 = z * y.T`` for some number `z` is called a *left* eigenvector of `a`,
1436 and, in general, the left and right eigenvectors of a matrix are not
1437 necessarily the (perhaps conjugate) transposes of each other.
1439 References
1440 ----------
1441 G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, FL,
1442 Academic Press, Inc., 1980, Various pp.
1444 Examples
1445 --------
1446 >>> import numpy as np
1447 >>> from numpy import linalg as LA
1449 (Almost) trivial example with real eigenvalues and eigenvectors.
1451 >>> eigenvalues, eigenvectors = LA.eig(np.diag((1, 2, 3)))
1452 >>> eigenvalues
1453 array([1., 2., 3.])
1454 >>> eigenvectors
1455 array([[1., 0., 0.],
1456 [0., 1., 0.],
1457 [0., 0., 1.]])
1459 Real matrix possessing complex eigenvalues and eigenvectors;
1460 note that the eigenvalues are complex conjugates of each other.
1462 >>> eigenvalues, eigenvectors = LA.eig(np.array([[1, -1], [1, 1]]))
1463 >>> eigenvalues
1464 array([1.+1.j, 1.-1.j])
1465 >>> eigenvectors
1466 array([[0.70710678+0.j , 0.70710678-0.j ],
1467 [0. -0.70710678j, 0. +0.70710678j]])
1469 Complex-valued matrix with real eigenvalues (but complex-valued
1470 eigenvectors); note that ``a.conj().T == a``, i.e., `a` is Hermitian.
1472 >>> a = np.array([[1, 1j], [-1j, 1]])
1473 >>> eigenvalues, eigenvectors = LA.eig(a)
1474 >>> eigenvalues
1475 array([2.+0.j, 0.+0.j])
1476 >>> eigenvectors
1477 array([[ 0. +0.70710678j, 0.70710678+0.j ], # may vary
1478 [ 0.70710678+0.j , -0. +0.70710678j]])
1480 Be careful about round-off error!
1482 >>> a = np.array([[1 + 1e-9, 0], [0, 1 - 1e-9]])
1483 >>> # Theor. eigenvalues are 1 +/- 1e-9
1484 >>> eigenvalues, eigenvectors = LA.eig(a)
1485 >>> eigenvalues
1486 array([1., 1.])
1487 >>> eigenvectors
1488 array([[1., 0.],
1489 [0., 1.]])
1491 """
1492 a, wrap = _makearray(a)
1493 _assert_stacked_square(a)
1494 _assert_finite(a)
1495 t, result_t = _commonType(a)
1497 signature = 'D->DD' if isComplexType(t) else 'd->DD'
1498 with errstate(call=_raise_linalgerror_eigenvalues_nonconvergence,
1499 invalid='call', over='ignore', divide='ignore',
1500 under='ignore'):
1501 w, vt = _umath_linalg.eig(a, signature=signature)
1503 if not isComplexType(t) and all(w.imag == 0.0):
1504 w = w.real
1505 vt = vt.real
1506 result_t = _realType(result_t)
1507 else:
1508 result_t = _complexType(result_t)
1510 vt = vt.astype(result_t, copy=False)
1511 return EigResult(w.astype(result_t, copy=False), wrap(vt))
1514@array_function_dispatch(_eigvalsh_dispatcher)
1515def eigh(a, UPLO='L'):
1516 """
1517 Return the eigenvalues and eigenvectors of a complex Hermitian
1518 (conjugate symmetric) or a real symmetric matrix.
1520 Returns two objects, a 1-D array containing the eigenvalues of `a`, and
1521 a 2-D square array or matrix (depending on the input type) of the
1522 corresponding eigenvectors (in columns).
1524 Parameters
1525 ----------
1526 a : (..., M, M) array
1527 Hermitian or real symmetric matrices whose eigenvalues and
1528 eigenvectors are to be computed.
1529 UPLO : {'L', 'U'}, optional
1530 Specifies whether the calculation is done with the lower triangular
1531 part of `a` ('L', default) or the upper triangular part ('U').
1532 Irrespective of this value only the real parts of the diagonal will
1533 be considered in the computation to preserve the notion of a Hermitian
1534 matrix. It therefore follows that the imaginary part of the diagonal
1535 will always be treated as zero.
1537 Returns
1538 -------
1539 A namedtuple with the following attributes:
1541 eigenvalues : (..., M) ndarray
1542 The eigenvalues in ascending order, each repeated according to
1543 its multiplicity.
1544 eigenvectors : {(..., M, M) ndarray, (..., M, M) matrix}
1545 The column ``eigenvectors[:, i]`` is the normalized eigenvector
1546 corresponding to the eigenvalue ``eigenvalues[i]``. Will return a
1547 matrix object if `a` is a matrix object.
1549 Raises
1550 ------
1551 LinAlgError
1552 If the eigenvalue computation does not converge.
1554 See Also
1555 --------
1556 eigvalsh : eigenvalues of real symmetric or complex Hermitian
1557 (conjugate symmetric) arrays.
1558 eig : eigenvalues and right eigenvectors for non-symmetric arrays.
1559 eigvals : eigenvalues of non-symmetric arrays.
1560 scipy.linalg.eigh : Similar function in SciPy (but also solves the
1561 generalized eigenvalue problem).
1563 Notes
1564 -----
1565 Broadcasting rules apply, see the `numpy.linalg` documentation for
1566 details.
1568 The eigenvalues/eigenvectors are computed using LAPACK routines ``_syevd``,
1569 ``_heevd``.
1571 The eigenvalues of real symmetric or complex Hermitian matrices are always
1572 real. [1]_ The array `eigenvalues` of (column) eigenvectors is unitary and
1573 `a`, `eigenvalues`, and `eigenvectors` satisfy the equations ``dot(a,
1574 eigenvectors[:, i]) = eigenvalues[i] * eigenvectors[:, i]``.
1576 References
1577 ----------
1578 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
1579 FL, Academic Press, Inc., 1980, pg. 222.
1581 Examples
1582 --------
1583 >>> import numpy as np
1584 >>> from numpy import linalg as LA
1585 >>> a = np.array([[1, -2j], [2j, 5]])
1586 >>> a
1587 array([[ 1.+0.j, -0.-2.j],
1588 [ 0.+2.j, 5.+0.j]])
1589 >>> eigenvalues, eigenvectors = LA.eigh(a)
1590 >>> eigenvalues
1591 array([0.17157288, 5.82842712])
1592 >>> eigenvectors
1593 array([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1594 [ 0. +0.38268343j, 0. -0.92387953j]])
1596 >>> (np.dot(a, eigenvectors[:, 0]) -
1597 ... eigenvalues[0] * eigenvectors[:, 0]) # verify 1st eigenval/vec pair
1598 array([5.55111512e-17+0.0000000e+00j, 0.00000000e+00+1.2490009e-16j])
1599 >>> (np.dot(a, eigenvectors[:, 1]) -
1600 ... eigenvalues[1] * eigenvectors[:, 1]) # verify 2nd eigenval/vec pair
1601 array([0.+0.j, 0.+0.j])
1603 >>> A = np.matrix(a) # what happens if input is a matrix object
1604 >>> A
1605 matrix([[ 1.+0.j, -0.-2.j],
1606 [ 0.+2.j, 5.+0.j]])
1607 >>> eigenvalues, eigenvectors = LA.eigh(A)
1608 >>> eigenvalues
1609 array([0.17157288, 5.82842712])
1610 >>> eigenvectors
1611 matrix([[-0.92387953+0.j , -0.38268343+0.j ], # may vary
1612 [ 0. +0.38268343j, 0. -0.92387953j]])
1614 >>> # demonstrate the treatment of the imaginary part of the diagonal
1615 >>> a = np.array([[5+2j, 9-2j], [0+2j, 2-1j]])
1616 >>> a
1617 array([[5.+2.j, 9.-2.j],
1618 [0.+2.j, 2.-1.j]])
1619 >>> # with UPLO='L' this is numerically equivalent to using LA.eig() with:
1620 >>> b = np.array([[5.+0.j, 0.-2.j], [0.+2.j, 2.-0.j]])
1621 >>> b
1622 array([[5.+0.j, 0.-2.j],
1623 [0.+2.j, 2.+0.j]])
1624 >>> wa, va = LA.eigh(a)
1625 >>> wb, vb = LA.eig(b)
1626 >>> wa
1627 array([1., 6.])
1628 >>> wb
1629 array([6.+0.j, 1.+0.j])
1630 >>> va
1631 array([[-0.4472136 +0.j , -0.89442719+0.j ], # may vary
1632 [ 0. +0.89442719j, 0. -0.4472136j ]])
1633 >>> vb
1634 array([[ 0.89442719+0.j , -0. +0.4472136j],
1635 [-0. +0.4472136j, 0.89442719+0.j ]])
1637 """
1638 UPLO = UPLO.upper()
1639 if UPLO not in ('L', 'U'):
1640 raise ValueError("UPLO argument must be 'L' or 'U'")
1642 a, wrap = _makearray(a)
1643 _assert_stacked_square(a)
1644 t, result_t = _commonType(a)
1646 if UPLO == 'L':
1647 gufunc = _umath_linalg.eigh_lo
1648 else:
1649 gufunc = _umath_linalg.eigh_up
1651 signature = 'D->dD' if isComplexType(t) else 'd->dd'
1652 with errstate(call=_raise_linalgerror_eigenvalues_nonconvergence,
1653 invalid='call', over='ignore', divide='ignore',
1654 under='ignore'):
1655 w, vt = gufunc(a, signature=signature)
1656 w = w.astype(_realType(result_t), copy=False)
1657 vt = vt.astype(result_t, copy=False)
1658 return EighResult(w, wrap(vt))
1661# Singular value decomposition
1663def _svd_dispatcher(a, full_matrices=None, compute_uv=None, hermitian=None):
1664 return (a,)
1667@array_function_dispatch(_svd_dispatcher)
1668def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
1669 """
1670 Singular Value Decomposition.
1672 When `a` is a 2D array, and ``full_matrices=False``, then it is
1673 factorized as ``u @ np.diag(s) @ vh = (u * s) @ vh``, where
1674 `u` and the Hermitian transpose of `vh` are 2D arrays with
1675 orthonormal columns and `s` is a 1D array of `a`'s singular
1676 values. When `a` is higher-dimensional, SVD is applied in
1677 stacked mode as explained below.
1679 Parameters
1680 ----------
1681 a : (..., M, N) array_like
1682 A real or complex array with ``a.ndim >= 2``.
1683 full_matrices : bool, optional
1684 If True (default), `u` and `vh` have the shapes ``(..., M, M)`` and
1685 ``(..., N, N)``, respectively. Otherwise, the shapes are
1686 ``(..., M, K)`` and ``(..., K, N)``, respectively, where
1687 ``K = min(M, N)``.
1688 compute_uv : bool, optional
1689 Whether or not to compute `u` and `vh` in addition to `s`. True
1690 by default.
1691 hermitian : bool, optional
1692 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
1693 enabling a more efficient method for finding singular values.
1694 Defaults to False.
1696 Returns
1697 -------
1698 U : { (..., M, M), (..., M, K) } array
1699 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1700 size as those of the input `a`. The size of the last two dimensions
1701 depends on the value of `full_matrices`. Only returned when
1702 `compute_uv` is True.
1703 S : (..., K) array
1704 Vector(s) with the singular values, within each vector sorted in
1705 descending order. The first ``a.ndim - 2`` dimensions have the same
1706 size as those of the input `a`.
1707 Vh : { (..., N, N), (..., K, N) } array
1708 Unitary array(s). The first ``a.ndim - 2`` dimensions have the same
1709 size as those of the input `a`. The size of the last two dimensions
1710 depends on the value of `full_matrices`. Only returned when
1711 `compute_uv` is True.
1713 Raises
1714 ------
1715 LinAlgError
1716 If SVD computation does not converge.
1718 See Also
1719 --------
1720 scipy.linalg.svd : Similar function in SciPy.
1721 scipy.linalg.svdvals : Compute singular values of a matrix.
1723 Notes
1724 -----
1725 When `compute_uv` is True, the result is a namedtuple with the following
1726 attribute names: `U`, `S`, and `Vh`.
1728 The decomposition is performed using LAPACK routine ``_gesdd``.
1730 SVD is usually described for the factorization of a 2D matrix :math:`A`.
1731 The higher-dimensional case will be discussed below. In the 2D case, SVD is
1732 written as :math:`A = U S V^H`, where :math:`A = a`, :math:`U= u`,
1733 :math:`S= \\mathtt{np.diag}(s)` and :math:`V^H = vh`. The 1D array `s`
1734 contains the singular values of `a` and `u` and `vh` are unitary. The rows
1735 of `vh` are the eigenvectors of :math:`A^H A` and the columns of `u` are
1736 the eigenvectors of :math:`A A^H`. In both cases the corresponding
1737 (possibly non-zero) eigenvalues are given by ``s**2``.
1739 If `a` has more than two dimensions, then broadcasting rules apply, as
1740 explained in :ref:`routines.linalg-broadcasting`. This means that SVD is
1741 working in "stacked" mode: it iterates over all indices of the first
1742 ``a.ndim - 2`` dimensions and for each combination SVD is applied to the
1743 last two indices. The matrix `a` can be reconstructed from the
1744 decomposition with either ``(u * s[..., None, :]) @ vh`` or
1745 ``u @ (s[..., None] * vh)``. (The ``@`` operator can be replaced by the
1746 function ``np.matmul`` for python versions below 3.5.)
1748 If `a` is a ``matrix`` object (as opposed to an ``ndarray``), then so are
1749 all the return values.
1751 Examples
1752 --------
1753 >>> import numpy as np
1754 >>> rng = np.random.default_rng()
1755 >>> a = rng.normal(size=(9, 6)) + 1j*rng.normal(size=(9, 6))
1756 >>> b = rng.normal(size=(2, 7, 8, 3)) + 1j*rng.normal(size=(2, 7, 8, 3))
1759 Reconstruction based on full SVD, 2D case:
1761 >>> U, S, Vh = np.linalg.svd(a, full_matrices=True)
1762 >>> U.shape, S.shape, Vh.shape
1763 ((9, 9), (6,), (6, 6))
1764 >>> np.allclose(a, np.dot(U[:, :6] * S, Vh))
1765 True
1766 >>> smat = np.zeros((9, 6), dtype=complex)
1767 >>> smat[:6, :6] = np.diag(S)
1768 >>> np.allclose(a, np.dot(U, np.dot(smat, Vh)))
1769 True
1771 Reconstruction based on reduced SVD, 2D case:
1773 >>> U, S, Vh = np.linalg.svd(a, full_matrices=False)
1774 >>> U.shape, S.shape, Vh.shape
1775 ((9, 6), (6,), (6, 6))
1776 >>> np.allclose(a, np.dot(U * S, Vh))
1777 True
1778 >>> smat = np.diag(S)
1779 >>> np.allclose(a, np.dot(U, np.dot(smat, Vh)))
1780 True
1782 Reconstruction based on full SVD, 4D case:
1784 >>> U, S, Vh = np.linalg.svd(b, full_matrices=True)
1785 >>> U.shape, S.shape, Vh.shape
1786 ((2, 7, 8, 8), (2, 7, 3), (2, 7, 3, 3))
1787 >>> np.allclose(b, np.matmul(U[..., :3] * S[..., None, :], Vh))
1788 True
1789 >>> np.allclose(b, np.matmul(U[..., :3], S[..., None] * Vh))
1790 True
1792 Reconstruction based on reduced SVD, 4D case:
1794 >>> U, S, Vh = np.linalg.svd(b, full_matrices=False)
1795 >>> U.shape, S.shape, Vh.shape
1796 ((2, 7, 8, 3), (2, 7, 3), (2, 7, 3, 3))
1797 >>> np.allclose(b, np.matmul(U * S[..., None, :], Vh))
1798 True
1799 >>> np.allclose(b, np.matmul(U, S[..., None] * Vh))
1800 True
1802 """
1803 import numpy as np
1804 a, wrap = _makearray(a)
1806 if hermitian:
1807 # note: lapack svd returns eigenvalues with s ** 2 sorted descending,
1808 # but eig returns s sorted ascending, so we re-order the eigenvalues
1809 # and related arrays to have the correct order
1810 if compute_uv:
1811 s, u = eigh(a)
1812 # avoid zero sign
1813 sgn = np.copysign(1.0, s)
1814 s = abs(s)
1815 sidx = argsort(s)[..., ::-1]
1816 sgn = np.take_along_axis(sgn, sidx, axis=-1)
1817 s = np.take_along_axis(s, sidx, axis=-1)
1818 u = np.take_along_axis(u, sidx[..., None, :], axis=-1)
1819 # singular values are unsigned, move the sign into v
1820 vt = transpose(u * sgn[..., None, :]).conjugate()
1821 return SVDResult(wrap(u), s, wrap(vt))
1822 else:
1823 s = eigvalsh(a)
1824 s = abs(s)
1825 return sort(s)[..., ::-1]
1827 _assert_stacked_2d(a)
1828 t, result_t = _commonType(a)
1830 m, n = a.shape[-2:]
1831 if compute_uv:
1832 if full_matrices:
1833 gufunc = _umath_linalg.svd_f
1834 else:
1835 gufunc = _umath_linalg.svd_s
1837 signature = 'D->DdD' if isComplexType(t) else 'd->ddd'
1838 with errstate(call=_raise_linalgerror_svd_nonconvergence,
1839 invalid='call', over='ignore', divide='ignore',
1840 under='ignore'):
1841 u, s, vh = gufunc(a, signature=signature)
1842 u = u.astype(result_t, copy=False)
1843 s = s.astype(_realType(result_t), copy=False)
1844 vh = vh.astype(result_t, copy=False)
1845 return SVDResult(wrap(u), s, wrap(vh))
1846 else:
1847 signature = 'D->d' if isComplexType(t) else 'd->d'
1848 with errstate(call=_raise_linalgerror_svd_nonconvergence,
1849 invalid='call', over='ignore', divide='ignore',
1850 under='ignore'):
1851 s = _umath_linalg.svd(a, signature=signature)
1852 s = s.astype(_realType(result_t), copy=False)
1853 return s
1856def _svdvals_dispatcher(x):
1857 return (x,)
1860@array_function_dispatch(_svdvals_dispatcher)
1861def svdvals(x, /):
1862 """
1863 Returns the singular values of a matrix (or a stack of matrices) ``x``.
1864 When x is a stack of matrices, the function will compute the singular
1865 values for each matrix in the stack.
1867 This function is Array API compatible.
1869 Calling ``np.svdvals(x)`` to get singular values is the same as
1870 ``np.svd(x, compute_uv=False, hermitian=False)``.
1872 Parameters
1873 ----------
1874 x : (..., M, N) array_like
1875 Input array having shape (..., M, N) and whose last two
1876 dimensions form matrices on which to perform singular value
1877 decomposition. Should have a floating-point data type.
1879 Returns
1880 -------
1881 out : ndarray
1882 An array with shape (..., K) that contains the vector(s)
1883 of singular values of length K, where K = min(M, N).
1885 See Also
1886 --------
1887 scipy.linalg.svdvals : Compute singular values of a matrix.
1889 Examples
1890 --------
1892 >>> np.linalg.svdvals([[1, 2, 3, 4, 5],
1893 ... [1, 4, 9, 16, 25],
1894 ... [1, 8, 27, 64, 125]])
1895 array([146.68862757, 5.57510612, 0.60393245])
1897 Determine the rank of a matrix using singular values:
1899 >>> s = np.linalg.svdvals([[1, 2, 3],
1900 ... [2, 4, 6],
1901 ... [-1, 1, -1]]); s
1902 array([8.38434191e+00, 1.64402274e+00, 2.31534378e-16])
1903 >>> np.count_nonzero(s > 1e-10) # Matrix of rank 2
1904 2
1906 """
1907 return svd(x, compute_uv=False, hermitian=False)
1910def _cond_dispatcher(x, p=None):
1911 return (x,)
1914@array_function_dispatch(_cond_dispatcher)
1915def cond(x, p=None):
1916 """
1917 Compute the condition number of a matrix.
1919 This function is capable of returning the condition number using
1920 one of seven different norms, depending on the value of `p` (see
1921 Parameters below).
1923 Parameters
1924 ----------
1925 x : (..., M, N) array_like
1926 The matrix whose condition number is sought.
1927 p : {None, 1, -1, 2, -2, inf, -inf, 'fro'}, optional
1928 Order of the norm used in the condition number computation:
1930 ===== ============================
1931 p norm for matrices
1932 ===== ============================
1933 None 2-norm, computed directly using the ``SVD``
1934 'fro' Frobenius norm
1935 inf max(sum(abs(x), axis=1))
1936 -inf min(sum(abs(x), axis=1))
1937 1 max(sum(abs(x), axis=0))
1938 -1 min(sum(abs(x), axis=0))
1939 2 2-norm (largest sing. value)
1940 -2 smallest singular value
1941 ===== ============================
1943 inf means the `numpy.inf` object, and the Frobenius norm is
1944 the root-of-sum-of-squares norm.
1946 Returns
1947 -------
1948 c : {float, inf}
1949 The condition number of the matrix. May be infinite.
1951 See Also
1952 --------
1953 numpy.linalg.norm
1955 Notes
1956 -----
1957 The condition number of `x` is defined as the norm of `x` times the
1958 norm of the inverse of `x` [1]_; the norm can be the usual L2-norm
1959 (root-of-sum-of-squares) or one of a number of other matrix norms.
1961 References
1962 ----------
1963 .. [1] G. Strang, *Linear Algebra and Its Applications*, Orlando, FL,
1964 Academic Press, Inc., 1980, pg. 285.
1966 Examples
1967 --------
1968 >>> import numpy as np
1969 >>> from numpy import linalg as LA
1970 >>> a = np.array([[1, 0, -1], [0, 1, 0], [1, 0, 1]])
1971 >>> a
1972 array([[ 1, 0, -1],
1973 [ 0, 1, 0],
1974 [ 1, 0, 1]])
1975 >>> LA.cond(a)
1976 1.4142135623730951
1977 >>> LA.cond(a, 'fro')
1978 3.1622776601683795
1979 >>> LA.cond(a, np.inf)
1980 2.0
1981 >>> LA.cond(a, -np.inf)
1982 1.0
1983 >>> LA.cond(a, 1)
1984 2.0
1985 >>> LA.cond(a, -1)
1986 1.0
1987 >>> LA.cond(a, 2)
1988 1.4142135623730951
1989 >>> LA.cond(a, -2)
1990 0.70710678118654746 # may vary
1991 >>> (min(LA.svd(a, compute_uv=False)) *
1992 ... min(LA.svd(LA.inv(a), compute_uv=False)))
1993 0.70710678118654746 # may vary
1995 """
1996 x = asarray(x) # in case we have a matrix
1997 if _is_empty_2d(x):
1998 raise LinAlgError("cond is not defined on empty arrays")
1999 if p is None or p in {2, -2}:
2000 s = svd(x, compute_uv=False)
2001 with errstate(all='ignore'):
2002 if p == -2:
2003 r = s[..., -1] / s[..., 0]
2004 else:
2005 r = s[..., 0] / s[..., -1]
2006 else:
2007 # Call inv(x) ignoring errors. The result array will
2008 # contain nans in the entries where inversion failed.
2009 _assert_stacked_square(x)
2010 t, result_t = _commonType(x)
2011 result_t = _realType(result_t) # condition number is always real
2012 signature = 'D->D' if isComplexType(t) else 'd->d'
2013 with errstate(all='ignore'):
2014 invx = _umath_linalg.inv(x, signature=signature)
2015 r = norm(x, p, axis=(-2, -1)) * norm(invx, p, axis=(-2, -1))
2016 r = r.astype(result_t, copy=False)
2018 # Convert nans to infs unless the original array had nan entries
2019 nan_mask = isnan(r)
2020 if nan_mask.any():
2021 nan_mask &= ~isnan(x).any(axis=(-2, -1))
2022 if r.ndim > 0:
2023 r[nan_mask] = inf
2024 elif nan_mask:
2025 # Convention is to return scalars instead of 0d arrays.
2026 r = r.dtype.type(inf)
2028 return r
2031def _matrix_rank_dispatcher(A, tol=None, hermitian=None, *, rtol=None):
2032 return (A,)
2035@array_function_dispatch(_matrix_rank_dispatcher)
2036def matrix_rank(A, tol=None, hermitian=False, *, rtol=None):
2037 """
2038 Return matrix rank of array using SVD method
2040 Rank of the array is the number of singular values of the array that are
2041 greater than `tol`.
2043 Parameters
2044 ----------
2045 A : {(M,), (..., M, N)} array_like
2046 Input vector or stack of matrices.
2047 tol : (...) array_like, float, optional
2048 Threshold below which SVD values are considered zero. If `tol` is
2049 None, and ``S`` is an array with singular values for `M`, and
2050 ``eps`` is the epsilon value for datatype of ``S``, then `tol` is
2051 set to ``S.max() * max(M, N) * eps``.
2052 hermitian : bool, optional
2053 If True, `A` is assumed to be Hermitian (symmetric if real-valued),
2054 enabling a more efficient method for finding singular values.
2055 Defaults to False.
2056 rtol : (...) array_like, float, optional
2057 Parameter for the relative tolerance component. Only ``tol`` or
2058 ``rtol`` can be set at a time. Defaults to ``max(M, N) * eps``.
2060 .. versionadded:: 2.0.0
2062 Returns
2063 -------
2064 rank : (...) array_like
2065 Rank of A.
2067 Notes
2068 -----
2069 The default threshold to detect rank deficiency is a test on the magnitude
2070 of the singular values of `A`. By default, we identify singular values
2071 less than ``S.max() * max(M, N) * eps`` as indicating rank deficiency
2072 (with the symbols defined above). This is the algorithm MATLAB uses [1]_.
2073 It also appears in *Numerical recipes* in the discussion of SVD solutions
2074 for linear least squares [2]_.
2076 This default threshold is designed to detect rank deficiency accounting
2077 for the numerical errors of the SVD computation. Imagine that there
2078 is a column in `A` that is an exact (in floating point) linear combination
2079 of other columns in `A`. Computing the SVD on `A` will not produce
2080 a singular value exactly equal to 0 in general: any difference of
2081 the smallest SVD value from 0 will be caused by numerical imprecision
2082 in the calculation of the SVD. Our threshold for small SVD values takes
2083 this numerical imprecision into account, and the default threshold will
2084 detect such numerical rank deficiency. The threshold may declare a matrix
2085 `A` rank deficient even if the linear combination of some columns of `A`
2086 is not exactly equal to another column of `A` but only numerically very
2087 close to another column of `A`.
2089 We chose our default threshold because it is in wide use. Other thresholds
2090 are possible. For example, elsewhere in the 2007 edition of *Numerical
2091 recipes* there is an alternative threshold of ``S.max() *
2092 np.finfo(A.dtype).eps / 2. * np.sqrt(m + n + 1.)``. The authors describe
2093 this threshold as being based on "expected roundoff error" (p 71).
2095 The thresholds above deal with floating point roundoff error in the
2096 calculation of the SVD. However, you may have more information about
2097 the sources of error in `A` that would make you consider other tolerance
2098 values to detect *effective* rank deficiency. The most useful measure
2099 of the tolerance depends on the operations you intend to use on your
2100 matrix. For example, if your data come from uncertain measurements with
2101 uncertainties greater than floating point epsilon, choosing a tolerance
2102 near that uncertainty may be preferable. The tolerance may be absolute
2103 if the uncertainties are absolute rather than relative.
2105 References
2106 ----------
2107 .. [1] MATLAB reference documentation, "Rank"
2108 https://www.mathworks.com/help/techdoc/ref/rank.html
2109 .. [2] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery,
2110 "Numerical Recipes (3rd edition)", Cambridge University Press, 2007,
2111 page 795.
2113 Examples
2114 --------
2115 >>> import numpy as np
2116 >>> from numpy.linalg import matrix_rank
2117 >>> matrix_rank(np.eye(4)) # Full rank matrix
2118 4
2119 >>> I=np.eye(4); I[-1,-1] = 0. # rank deficient matrix
2120 >>> matrix_rank(I)
2121 3
2122 >>> matrix_rank(np.ones((4,))) # 1 dimension - rank 1 unless all 0
2123 1
2124 >>> matrix_rank(np.zeros((4,)))
2125 0
2126 """
2127 if rtol is not None and tol is not None:
2128 raise ValueError("`tol` and `rtol` can't be both set.")
2130 A = asarray(A)
2131 if A.ndim < 2:
2132 return int(not all(A == 0))
2134 S = svd(A, compute_uv=False, hermitian=hermitian)
2136 if tol is None:
2137 if rtol is None:
2138 rtol = max(A.shape[-2:]) * finfo(S.dtype).eps
2139 else:
2140 rtol = asarray(rtol)[..., newaxis]
2141 tol = S.max(axis=-1, keepdims=True, initial=0) * rtol
2142 else:
2143 tol = asarray(tol)[..., newaxis]
2145 return count_nonzero(S > tol, axis=-1)
2148# Generalized inverse
2150def _pinv_dispatcher(a, rcond=None, hermitian=None, *, rtol=None):
2151 return (a,)
2154@array_function_dispatch(_pinv_dispatcher)
2155def pinv(a, rcond=None, hermitian=False, *, rtol=_NoValue):
2156 """
2157 Compute the (Moore-Penrose) pseudo-inverse of a matrix.
2159 Calculate the generalized inverse of a matrix using its
2160 singular-value decomposition (SVD) and including all
2161 *large* singular values.
2163 Parameters
2164 ----------
2165 a : (..., M, N) array_like
2166 Matrix or stack of matrices to be pseudo-inverted.
2167 rcond : (...) array_like of float, optional
2168 Cutoff for small singular values.
2169 Singular values less than or equal to
2170 ``rcond * largest_singular_value`` are set to zero.
2171 Broadcasts against the stack of matrices. Default: ``1e-15``.
2172 hermitian : bool, optional
2173 If True, `a` is assumed to be Hermitian (symmetric if real-valued),
2174 enabling a more efficient method for finding singular values.
2175 Defaults to False.
2176 rtol : (...) array_like of float, optional
2177 Same as `rcond`, but it's an Array API compatible parameter name.
2178 Only `rcond` or `rtol` can be set at a time. If none of them are
2179 provided then NumPy's ``1e-15`` default is used. If ``rtol=None``
2180 is passed then the API standard default is used.
2182 .. versionadded:: 2.0.0
2184 Returns
2185 -------
2186 B : (..., N, M) ndarray
2187 The pseudo-inverse of `a`. If `a` is a `matrix` instance, then so
2188 is `B`.
2190 Raises
2191 ------
2192 LinAlgError
2193 If the SVD computation does not converge.
2195 See Also
2196 --------
2197 scipy.linalg.pinv : Similar function in SciPy.
2198 scipy.linalg.pinvh : Compute the (Moore-Penrose) pseudo-inverse of a
2199 Hermitian matrix.
2201 Notes
2202 -----
2203 The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
2204 defined as: "the matrix that 'solves' [the least-squares problem]
2205 :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then
2206 :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`.
2208 It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular
2209 value decomposition of A, then
2210 :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are
2211 orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting
2212 of A's so-called singular values, (followed, typically, by
2213 zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix
2214 consisting of the reciprocals of A's singular values
2215 (again, followed by zeros). [1]_
2217 References
2218 ----------
2219 .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
2220 FL, Academic Press, Inc., 1980, pp. 139-142.
2222 Examples
2223 --------
2224 The following example checks that ``a * a+ * a == a`` and
2225 ``a+ * a * a+ == a+``:
2227 >>> import numpy as np
2228 >>> rng = np.random.default_rng()
2229 >>> a = rng.normal(size=(9, 6))
2230 >>> B = np.linalg.pinv(a)
2231 >>> np.allclose(a, np.dot(a, np.dot(B, a)))
2232 True
2233 >>> np.allclose(B, np.dot(B, np.dot(a, B)))
2234 True
2236 """
2237 a, wrap = _makearray(a)
2238 if rcond is None:
2239 if rtol is _NoValue:
2240 rcond = 1e-15
2241 elif rtol is None:
2242 rcond = max(a.shape[-2:]) * finfo(a.dtype).eps
2243 else:
2244 rcond = rtol
2245 elif rtol is not _NoValue:
2246 raise ValueError("`rtol` and `rcond` can't be both set.")
2247 else:
2248 # NOTE: Deprecate `rcond` in a few versions.
2249 pass
2251 rcond = asarray(rcond)
2252 if _is_empty_2d(a):
2253 m, n = a.shape[-2:]
2254 res = empty(a.shape[:-2] + (n, m), dtype=a.dtype)
2255 return wrap(res)
2256 a = a.conjugate()
2257 u, s, vt = svd(a, full_matrices=False, hermitian=hermitian)
2259 # discard small singular values
2260 cutoff = rcond[..., newaxis] * amax(s, axis=-1, keepdims=True)
2261 large = s > cutoff
2262 s = divide(1, s, where=large, out=s)
2263 s[~large] = 0
2265 res = matmul(transpose(vt), multiply(s[..., newaxis], transpose(u)))
2266 return wrap(res)
2269# Determinant
2272@array_function_dispatch(_unary_dispatcher)
2273def slogdet(a):
2274 """
2275 Compute the sign and (natural) logarithm of the determinant of an array.
2277 If an array has a very small or very large determinant, then a call to
2278 `det` may overflow or underflow. This routine is more robust against such
2279 issues, because it computes the logarithm of the determinant rather than
2280 the determinant itself.
2282 Parameters
2283 ----------
2284 a : (..., M, M) array_like
2285 Input array, has to be a square 2-D array.
2287 Returns
2288 -------
2289 A namedtuple with the following attributes:
2291 sign : (...) array_like
2292 A number representing the sign of the determinant. For a real matrix,
2293 this is 1, 0, or -1. For a complex matrix, this is a complex number
2294 with absolute value 1 (i.e., it is on the unit circle), or else 0.
2295 logabsdet : (...) array_like
2296 The natural log of the absolute value of the determinant.
2298 If the determinant is zero, then `sign` will be 0 and `logabsdet`
2299 will be -inf. In all cases, the determinant is equal to
2300 ``sign * np.exp(logabsdet)``.
2302 See Also
2303 --------
2304 det
2306 Notes
2307 -----
2308 Broadcasting rules apply, see the `numpy.linalg` documentation for
2309 details.
2311 The determinant is computed via LU factorization using the LAPACK
2312 routine ``z/dgetrf``.
2314 Examples
2315 --------
2316 The determinant of a 2-D array ``[[a, b], [c, d]]`` is ``ad - bc``:
2318 >>> import numpy as np
2319 >>> a = np.array([[1, 2], [3, 4]])
2320 >>> (sign, logabsdet) = np.linalg.slogdet(a)
2321 >>> (sign, logabsdet)
2322 (-1, 0.69314718055994529) # may vary
2323 >>> sign * np.exp(logabsdet)
2324 -2.0
2326 Computing log-determinants for a stack of matrices:
2328 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2329 >>> a.shape
2330 (3, 2, 2)
2331 >>> sign, logabsdet = np.linalg.slogdet(a)
2332 >>> (sign, logabsdet)
2333 (array([-1., -1., -1.]), array([ 0.69314718, 1.09861229, 2.07944154]))
2334 >>> sign * np.exp(logabsdet)
2335 array([-2., -3., -8.])
2337 This routine succeeds where ordinary `det` does not:
2339 >>> np.linalg.det(np.eye(500) * 0.1)
2340 0.0
2341 >>> np.linalg.slogdet(np.eye(500) * 0.1)
2342 (1, -1151.2925464970228)
2344 """
2345 a = asarray(a)
2346 _assert_stacked_square(a)
2347 t, result_t = _commonType(a)
2348 real_t = _realType(result_t)
2349 signature = 'D->Dd' if isComplexType(t) else 'd->dd'
2350 sign, logdet = _umath_linalg.slogdet(a, signature=signature)
2351 sign = sign.astype(result_t, copy=False)
2352 logdet = logdet.astype(real_t, copy=False)
2353 return SlogdetResult(sign, logdet)
2356@array_function_dispatch(_unary_dispatcher)
2357def det(a):
2358 """
2359 Compute the determinant of an array.
2361 Parameters
2362 ----------
2363 a : (..., M, M) array_like
2364 Input array to compute determinants for.
2366 Returns
2367 -------
2368 det : (...) array_like
2369 Determinant of `a`.
2371 See Also
2372 --------
2373 slogdet : Another way to represent the determinant, more suitable
2374 for large matrices where underflow/overflow may occur.
2375 scipy.linalg.det : Similar function in SciPy.
2377 Notes
2378 -----
2379 Broadcasting rules apply, see the `numpy.linalg` documentation for
2380 details.
2382 The determinant is computed via LU factorization using the LAPACK
2383 routine ``z/dgetrf``.
2385 Examples
2386 --------
2387 The determinant of a 2-D array [[a, b], [c, d]] is ad - bc:
2389 >>> import numpy as np
2390 >>> a = np.array([[1, 2], [3, 4]])
2391 >>> np.linalg.det(a)
2392 -2.0 # may vary
2394 Computing determinants for a stack of matrices:
2396 >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
2397 >>> a.shape
2398 (3, 2, 2)
2399 >>> np.linalg.det(a)
2400 array([-2., -3., -8.])
2402 """
2403 a = asarray(a)
2404 _assert_stacked_square(a)
2405 t, result_t = _commonType(a)
2406 signature = 'D->D' if isComplexType(t) else 'd->d'
2407 r = _umath_linalg.det(a, signature=signature)
2408 r = r.astype(result_t, copy=False)
2409 return r
2412# Linear Least Squares
2414def _lstsq_dispatcher(a, b, rcond=None):
2415 return (a, b)
2418@array_function_dispatch(_lstsq_dispatcher)
2419def lstsq(a, b, rcond=None):
2420 r"""
2421 Return the least-squares solution to a linear matrix equation.
2423 Computes the vector `x` that approximately solves the equation
2424 ``a @ x = b``. The equation may be under-, well-, or over-determined
2425 (i.e., the number of linearly independent rows of `a` can be less than,
2426 equal to, or greater than its number of linearly independent columns).
2427 If `a` is square and of full rank, then `x` (but for round-off error)
2428 is the "exact" solution of the equation. Else, `x` minimizes the
2429 Euclidean 2-norm :math:`||b - ax||`. If there are multiple minimizing
2430 solutions, the one with the smallest 2-norm :math:`||x||` is returned.
2432 Parameters
2433 ----------
2434 a : (M, N) array_like
2435 "Coefficient" matrix.
2436 b : {(M,), (M, K)} array_like
2437 Ordinate or "dependent variable" values. If `b` is two-dimensional,
2438 the least-squares solution is calculated for each of the `K` columns
2439 of `b`.
2440 rcond : float, optional
2441 Cut-off ratio for small singular values of `a`.
2442 For the purposes of rank determination, singular values are treated
2443 as zero if they are smaller than `rcond` times the largest singular
2444 value of `a`.
2445 The default uses the machine precision times ``max(M, N)``. Passing
2446 ``-1`` will use machine precision.
2448 .. versionchanged:: 2.0
2449 Previously, the default was ``-1``, but a warning was given that
2450 this would change.
2452 Returns
2453 -------
2454 x : {(N,), (N, K)} ndarray
2455 Least-squares solution. If `b` is two-dimensional,
2456 the solutions are in the `K` columns of `x`.
2457 residuals : {(1,), (K,), (0,)} ndarray
2458 Sums of squared residuals: Squared Euclidean 2-norm for each column in
2459 ``b - a @ x``.
2460 If the rank of `a` is < N or M <= N, this is an empty array.
2461 If `b` is 1-dimensional, this is a (1,) shape array.
2462 Otherwise the shape is (K,).
2463 rank : int
2464 Rank of matrix `a`.
2465 s : (min(M, N),) ndarray
2466 Singular values of `a`.
2468 Raises
2469 ------
2470 LinAlgError
2471 If computation does not converge.
2473 See Also
2474 --------
2475 scipy.linalg.lstsq : Similar function in SciPy.
2477 Notes
2478 -----
2479 If `b` is a matrix, then all array results are returned as matrices.
2481 Examples
2482 --------
2483 Fit a line, ``y = mx + c``, through some noisy data-points:
2485 >>> import numpy as np
2486 >>> x = np.array([0, 1, 2, 3])
2487 >>> y = np.array([-1, 0.2, 0.9, 2.1])
2489 By examining the coefficients, we see that the line should have a
2490 gradient of roughly 1 and cut the y-axis at, more or less, -1.
2492 We can rewrite the line equation as ``y = Ap``, where ``A = [[x 1]]``
2493 and ``p = [[m], [c]]``. Now use `lstsq` to solve for `p`:
2495 >>> A = np.vstack([x, np.ones(len(x))]).T
2496 >>> A
2497 array([[ 0., 1.],
2498 [ 1., 1.],
2499 [ 2., 1.],
2500 [ 3., 1.]])
2502 >>> m, c = np.linalg.lstsq(A, y)[0]
2503 >>> m, c
2504 (1.0 -0.95) # may vary
2506 Plot the data along with the fitted line:
2508 >>> import matplotlib.pyplot as plt
2509 >>> _ = plt.plot(x, y, 'o', label='Original data', markersize=10)
2510 >>> _ = plt.plot(x, m*x + c, 'r', label='Fitted line')
2511 >>> _ = plt.legend()
2512 >>> plt.show()
2514 """
2515 a, _ = _makearray(a)
2516 b, wrap = _makearray(b)
2517 is_1d = b.ndim == 1
2518 if is_1d:
2519 b = b[:, newaxis]
2520 _assert_2d(a, b)
2521 m, n = a.shape[-2:]
2522 m2, n_rhs = b.shape[-2:]
2523 if m != m2:
2524 raise LinAlgError('Incompatible dimensions')
2526 t, result_t = _commonType(a, b)
2527 result_real_t = _realType(result_t)
2529 if rcond is None:
2530 rcond = finfo(t).eps * max(n, m)
2532 signature = 'DDd->Ddid' if isComplexType(t) else 'ddd->ddid'
2533 if n_rhs == 0:
2534 # lapack can't handle n_rhs = 0 - so allocate
2535 # the array one larger in that axis
2536 b = zeros(b.shape[:-2] + (m, n_rhs + 1), dtype=b.dtype)
2538 with errstate(call=_raise_linalgerror_lstsq, invalid='call',
2539 over='ignore', divide='ignore', under='ignore'):
2540 x, resids, rank, s = _umath_linalg.lstsq(a, b, rcond,
2541 signature=signature)
2542 if m == 0:
2543 x[...] = 0
2544 if n_rhs == 0:
2545 # remove the item we added
2546 x = x[..., :n_rhs]
2547 resids = resids[..., :n_rhs]
2549 # remove the axis we added
2550 if is_1d:
2551 x = x.squeeze(axis=-1)
2552 # we probably should squeeze resids too, but we can't
2553 # without breaking compatibility.
2555 # as documented
2556 if rank != n or m <= n:
2557 resids = array([], result_real_t)
2559 # coerce output arrays
2560 s = s.astype(result_real_t, copy=False)
2561 resids = resids.astype(result_real_t, copy=False)
2562 # Copying lets the memory in r_parts be freed
2563 x = x.astype(result_t, copy=True)
2564 return wrap(x), wrap(resids), rank, s
2567def _multi_svd_norm(x, row_axis, col_axis, op, initial=None):
2568 """Compute a function of the singular values of the 2-D matrices in `x`.
2570 This is a private utility function used by `numpy.linalg.norm()`.
2572 Parameters
2573 ----------
2574 x : ndarray
2575 row_axis, col_axis : int
2576 The axes of `x` that hold the 2-D matrices.
2577 op : callable
2578 This should be either numpy.amin or `numpy.amax` or `numpy.sum`.
2580 Returns
2581 -------
2582 result : float or ndarray
2583 If `x` is 2-D, the return values is a float.
2584 Otherwise, it is an array with ``x.ndim - 2`` dimensions.
2585 The return values are either the minimum or maximum or sum of the
2586 singular values of the matrices, depending on whether `op`
2587 is `numpy.amin` or `numpy.amax` or `numpy.sum`.
2589 """
2590 y = moveaxis(x, (row_axis, col_axis), (-2, -1))
2591 result = op(svd(y, compute_uv=False), axis=-1, initial=initial)
2592 return result
2595def _norm_dispatcher(x, ord=None, axis=None, keepdims=None):
2596 return (x,)
2599@array_function_dispatch(_norm_dispatcher)
2600def norm(x, ord=None, axis=None, keepdims=False):
2601 """
2602 Matrix or vector norm.
2604 This function is able to return one of eight different matrix norms,
2605 or one of an infinite number of vector norms (described below), depending
2606 on the value of the ``ord`` parameter.
2608 Parameters
2609 ----------
2610 x : array_like
2611 Input array. If `axis` is None, `x` must be 1-D or 2-D, unless `ord`
2612 is None. If both `axis` and `ord` are None, the 2-norm of
2613 ``x.ravel`` will be returned.
2614 ord : {int, float, inf, -inf, 'fro', 'nuc'}, optional
2615 Order of the norm (see table under ``Notes`` for what values are
2616 supported for matrices and vectors respectively). inf means numpy's
2617 `inf` object. The default is None.
2618 axis : {None, int, 2-tuple of ints}, optional.
2619 If `axis` is an integer, it specifies the axis of `x` along which to
2620 compute the vector norms. If `axis` is a 2-tuple, it specifies the
2621 axes that hold 2-D matrices, and the matrix norms of these matrices
2622 are computed. If `axis` is None then either a vector norm (when `x`
2623 is 1-D) or a matrix norm (when `x` is 2-D) is returned. The default
2624 is None.
2626 keepdims : bool, optional
2627 If this is set to True, the axes which are normed over are left in the
2628 result as dimensions with size one. With this option the result will
2629 broadcast correctly against the original `x`.
2631 Returns
2632 -------
2633 n : float or ndarray
2634 Norm of the matrix or vector(s).
2636 See Also
2637 --------
2638 scipy.linalg.norm : Similar function in SciPy.
2640 Notes
2641 -----
2642 For values of ``ord < 1``, the result is, strictly speaking, not a
2643 mathematical 'norm', but it may still be useful for various numerical
2644 purposes.
2646 The following norms can be calculated:
2648 ===== ============================ ==========================
2649 ord norm for matrices norm for vectors
2650 ===== ============================ ==========================
2651 None Frobenius norm 2-norm
2652 'fro' Frobenius norm --
2653 'nuc' nuclear norm --
2654 inf max(sum(abs(x), axis=1)) max(abs(x))
2655 -inf min(sum(abs(x), axis=1)) min(abs(x))
2656 0 -- sum(x != 0)
2657 1 max(sum(abs(x), axis=0)) as below
2658 -1 min(sum(abs(x), axis=0)) as below
2659 2 2-norm (largest sing. value) as below
2660 -2 smallest singular value as below
2661 other -- sum(abs(x)**ord)**(1./ord)
2662 ===== ============================ ==========================
2664 The Frobenius norm is given by [1]_:
2666 :math:`||A||_F = [\\sum_{i,j} abs(a_{i,j})^2]^{1/2}`
2668 The nuclear norm is the sum of the singular values.
2670 Both the Frobenius and nuclear norm orders are only defined for
2671 matrices and raise a ValueError when ``x.ndim != 2``.
2673 References
2674 ----------
2675 .. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*,
2676 Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15
2678 Examples
2679 --------
2681 >>> import numpy as np
2682 >>> from numpy import linalg as LA
2683 >>> a = np.arange(9) - 4
2684 >>> a
2685 array([-4, -3, -2, ..., 2, 3, 4])
2686 >>> b = a.reshape((3, 3))
2687 >>> b
2688 array([[-4, -3, -2],
2689 [-1, 0, 1],
2690 [ 2, 3, 4]])
2692 >>> LA.norm(a)
2693 7.745966692414834
2694 >>> LA.norm(b)
2695 7.745966692414834
2696 >>> LA.norm(b, 'fro')
2697 7.745966692414834
2698 >>> LA.norm(a, np.inf)
2699 4.0
2700 >>> LA.norm(b, np.inf)
2701 9.0
2702 >>> LA.norm(a, -np.inf)
2703 0.0
2704 >>> LA.norm(b, -np.inf)
2705 2.0
2707 >>> LA.norm(a, 1)
2708 20.0
2709 >>> LA.norm(b, 1)
2710 7.0
2711 >>> LA.norm(a, -1)
2712 -4.6566128774142013e-010
2713 >>> LA.norm(b, -1)
2714 6.0
2715 >>> LA.norm(a, 2)
2716 7.745966692414834
2717 >>> LA.norm(b, 2)
2718 7.3484692283495345
2720 >>> LA.norm(a, -2)
2721 0.0
2722 >>> LA.norm(b, -2)
2723 1.8570331885190563e-016 # may vary
2724 >>> LA.norm(a, 3)
2725 5.8480354764257312 # may vary
2726 >>> LA.norm(a, -3)
2727 0.0
2729 Using the `axis` argument to compute vector norms:
2731 >>> c = np.array([[ 1, 2, 3],
2732 ... [-1, 1, 4]])
2733 >>> LA.norm(c, axis=0)
2734 array([ 1.41421356, 2.23606798, 5. ])
2735 >>> LA.norm(c, axis=1)
2736 array([ 3.74165739, 4.24264069])
2737 >>> LA.norm(c, ord=1, axis=1)
2738 array([ 6., 6.])
2740 Using the `axis` argument to compute matrix norms:
2742 >>> m = np.arange(8).reshape(2,2,2)
2743 >>> LA.norm(m, axis=(1,2))
2744 array([ 3.74165739, 11.22497216])
2745 >>> LA.norm(m[0, :, :]), LA.norm(m[1, :, :])
2746 (3.7416573867739413, 11.224972160321824)
2748 """
2749 x = asarray(x)
2751 if not issubclass(x.dtype.type, (inexact, object_)):
2752 x = x.astype(float)
2754 # Immediately handle some default, simple, fast, and common cases.
2755 if axis is None:
2756 ndim = x.ndim
2757 if (
2758 (ord is None) or
2759 (ord in ('f', 'fro') and ndim == 2) or
2760 (ord == 2 and ndim == 1)
2761 ):
2762 x = x.ravel(order='K')
2763 if isComplexType(x.dtype.type):
2764 x_real = x.real
2765 x_imag = x.imag
2766 sqnorm = x_real.dot(x_real) + x_imag.dot(x_imag)
2767 else:
2768 sqnorm = x.dot(x)
2769 ret = sqrt(sqnorm)
2770 if keepdims:
2771 ret = ret.reshape(ndim * [1])
2772 return ret
2774 # Normalize the `axis` argument to a tuple.
2775 nd = x.ndim
2776 if axis is None:
2777 axis = tuple(range(nd))
2778 elif not isinstance(axis, tuple):
2779 try:
2780 axis = int(axis)
2781 except Exception as e:
2782 raise TypeError(
2783 "'axis' must be None, an integer or a tuple of integers"
2784 ) from e
2785 axis = (axis,)
2787 if len(axis) == 1:
2788 if ord == inf:
2789 return abs(x).max(axis=axis, keepdims=keepdims, initial=0)
2790 elif ord == -inf:
2791 return abs(x).min(axis=axis, keepdims=keepdims)
2792 elif ord == 0:
2793 # Zero norm
2794 return (
2795 (x != 0)
2796 .astype(x.real.dtype)
2797 .sum(axis=axis, keepdims=keepdims)
2798 )
2799 elif ord == 1:
2800 # special case for speedup
2801 return add.reduce(abs(x), axis=axis, keepdims=keepdims)
2802 elif ord is None or ord == 2:
2803 # special case for speedup
2804 s = (x.conj() * x).real
2805 return sqrt(add.reduce(s, axis=axis, keepdims=keepdims))
2806 # None of the str-type keywords for ord ('fro', 'nuc')
2807 # are valid for vectors
2808 elif isinstance(ord, str):
2809 raise ValueError(f"Invalid norm order '{ord}' for vectors")
2810 else:
2811 absx = abs(x)
2812 absx **= ord
2813 ret = add.reduce(absx, axis=axis, keepdims=keepdims)
2814 ret **= reciprocal(ord, dtype=ret.dtype)
2815 return ret
2816 elif len(axis) == 2:
2817 row_axis, col_axis = axis
2818 row_axis = normalize_axis_index(row_axis, nd)
2819 col_axis = normalize_axis_index(col_axis, nd)
2820 if row_axis == col_axis:
2821 raise ValueError('Duplicate axes given.')
2822 if ord == 2:
2823 ret = _multi_svd_norm(x, row_axis, col_axis, amax, 0)
2824 elif ord == -2:
2825 ret = _multi_svd_norm(x, row_axis, col_axis, amin)
2826 elif ord == 1:
2827 if col_axis > row_axis:
2828 col_axis -= 1
2829 ret = add.reduce(abs(x), axis=row_axis).max(axis=col_axis, initial=0)
2830 elif ord == inf:
2831 if row_axis > col_axis:
2832 row_axis -= 1
2833 ret = add.reduce(abs(x), axis=col_axis).max(axis=row_axis, initial=0)
2834 elif ord == -1:
2835 if col_axis > row_axis:
2836 col_axis -= 1
2837 ret = add.reduce(abs(x), axis=row_axis).min(axis=col_axis)
2838 elif ord == -inf:
2839 if row_axis > col_axis:
2840 row_axis -= 1
2841 ret = add.reduce(abs(x), axis=col_axis).min(axis=row_axis)
2842 elif ord in [None, 'fro', 'f']:
2843 ret = sqrt(add.reduce((x.conj() * x).real, axis=axis))
2844 elif ord == 'nuc':
2845 ret = _multi_svd_norm(x, row_axis, col_axis, sum, 0)
2846 else:
2847 raise ValueError("Invalid norm order for matrices.")
2848 if keepdims:
2849 ret_shape = list(x.shape)
2850 ret_shape[axis[0]] = 1
2851 ret_shape[axis[1]] = 1
2852 ret = ret.reshape(ret_shape)
2853 return ret
2854 else:
2855 raise ValueError("Improper number of dimensions to norm.")
2858# multi_dot
2860def _multidot_dispatcher(arrays, *, out=None):
2861 yield from arrays
2862 yield out
2865@array_function_dispatch(_multidot_dispatcher)
2866def multi_dot(arrays, *, out=None):
2867 """
2868 Compute the dot product of two or more arrays in a single function call,
2869 while automatically selecting the fastest evaluation order.
2871 `multi_dot` chains `numpy.dot` and uses optimal parenthesization
2872 of the matrices [1]_ [2]_. Depending on the shapes of the matrices,
2873 this can speed up the multiplication a lot.
2875 If the first argument is 1-D it is treated as a row vector.
2876 If the last argument is 1-D it is treated as a column vector.
2877 The other arguments must be 2-D.
2879 Think of `multi_dot` as::
2881 def multi_dot(arrays): return functools.reduce(np.dot, arrays)
2884 Parameters
2885 ----------
2886 arrays : sequence of array_like
2887 If the first argument is 1-D it is treated as row vector.
2888 If the last argument is 1-D it is treated as column vector.
2889 The other arguments must be 2-D.
2890 out : ndarray, optional
2891 Output argument. This must have the exact kind that would be returned
2892 if it was not used. In particular, it must have the right type, must be
2893 C-contiguous, and its dtype must be the dtype that would be returned
2894 for `dot(a, b)`. This is a performance feature. Therefore, if these
2895 conditions are not met, an exception is raised, instead of attempting
2896 to be flexible.
2898 Returns
2899 -------
2900 output : ndarray
2901 Returns the dot product of the supplied arrays.
2903 See Also
2904 --------
2905 numpy.dot : dot multiplication with two arguments.
2907 References
2908 ----------
2910 .. [1] Cormen, "Introduction to Algorithms", Chapter 15.2, p. 370-378
2911 .. [2] https://en.wikipedia.org/wiki/Matrix_chain_multiplication
2913 Examples
2914 --------
2915 `multi_dot` allows you to write::
2917 >>> import numpy as np
2918 >>> from numpy.linalg import multi_dot
2919 >>> # Prepare some data
2920 >>> A = np.random.random((10000, 100))
2921 >>> B = np.random.random((100, 1000))
2922 >>> C = np.random.random((1000, 5))
2923 >>> D = np.random.random((5, 333))
2924 >>> # the actual dot multiplication
2925 >>> _ = multi_dot([A, B, C, D])
2927 instead of::
2929 >>> _ = np.dot(np.dot(np.dot(A, B), C), D)
2930 >>> # or
2931 >>> _ = A.dot(B).dot(C).dot(D)
2933 Notes
2934 -----
2935 The cost for a matrix multiplication can be calculated with the
2936 following function::
2938 def cost(A, B):
2939 return A.shape[0] * A.shape[1] * B.shape[1]
2941 Assume we have three matrices
2942 :math:`A_{10 \\times 100}, B_{100 \\times 5}, C_{5 \\times 50}`.
2944 The costs for the two different parenthesizations are as follows::
2946 cost((AB)C) = 10*100*5 + 10*5*50 = 5000 + 2500 = 7500
2947 cost(A(BC)) = 10*100*50 + 100*5*50 = 50000 + 25000 = 75000
2949 """
2950 n = len(arrays)
2951 # optimization only makes sense for len(arrays) > 2
2952 if n < 2:
2953 raise ValueError("Expecting at least two arrays.")
2954 elif n == 2:
2955 return dot(arrays[0], arrays[1], out=out)
2957 arrays = [asanyarray(a) for a in arrays]
2959 # save original ndim to reshape the result array into the proper form later
2960 ndim_first, ndim_last = arrays[0].ndim, arrays[-1].ndim
2961 # Explicitly convert vectors to 2D arrays to keep the logic of the internal
2962 # _multi_dot_* functions as simple as possible.
2963 if arrays[0].ndim == 1:
2964 arrays[0] = atleast_2d(arrays[0])
2965 if arrays[-1].ndim == 1:
2966 arrays[-1] = atleast_2d(arrays[-1]).T
2967 _assert_2d(*arrays)
2969 # _multi_dot_three is much faster than _multi_dot_matrix_chain_order
2970 if n == 3:
2971 result = _multi_dot_three(arrays[0], arrays[1], arrays[2], out=out)
2972 else:
2973 order = _multi_dot_matrix_chain_order(arrays)
2974 result = _multi_dot(arrays, order, 0, n - 1, out=out)
2976 # return proper shape
2977 if ndim_first == 1 and ndim_last == 1:
2978 return result[0, 0] # scalar
2979 elif ndim_first == 1 or ndim_last == 1:
2980 return result.ravel() # 1-D
2981 else:
2982 return result
2985def _multi_dot_three(A, B, C, out=None):
2986 """
2987 Find the best order for three arrays and do the multiplication.
2989 For three arguments `_multi_dot_three` is approximately 15 times faster
2990 than `_multi_dot_matrix_chain_order`
2992 """
2993 a0, a1b0 = A.shape
2994 b1c0, c1 = C.shape
2995 # cost1 = cost((AB)C) = a0*a1b0*b1c0 + a0*b1c0*c1
2996 cost1 = a0 * b1c0 * (a1b0 + c1)
2997 # cost2 = cost(A(BC)) = a1b0*b1c0*c1 + a0*a1b0*c1
2998 cost2 = a1b0 * c1 * (a0 + b1c0)
3000 if cost1 < cost2:
3001 return dot(dot(A, B), C, out=out)
3002 else:
3003 return dot(A, dot(B, C), out=out)
3006def _multi_dot_matrix_chain_order(arrays, return_costs=False):
3007 """
3008 Return a np.array that encodes the optimal order of multiplications.
3010 The optimal order array is then used by `_multi_dot()` to do the
3011 multiplication.
3013 Also return the cost matrix if `return_costs` is `True`
3015 The implementation CLOSELY follows Cormen, "Introduction to Algorithms",
3016 Chapter 15.2, p. 370-378. Note that Cormen uses 1-based indices.
3018 cost[i, j] = min([
3019 cost[prefix] + cost[suffix] + cost_mult(prefix, suffix)
3020 for k in range(i, j)])
3022 """
3023 n = len(arrays)
3024 # p stores the dimensions of the matrices
3025 # Example for p: A_{10x100}, B_{100x5}, C_{5x50} --> p = [10, 100, 5, 50]
3026 p = [a.shape[0] for a in arrays] + [arrays[-1].shape[1]]
3027 # m is a matrix of costs of the subproblems
3028 # m[i,j]: min number of scalar multiplications needed to compute A_{i..j}
3029 m = zeros((n, n), dtype=double)
3030 # s is the actual ordering
3031 # s[i, j] is the value of k at which we split the product A_i..A_j
3032 s = empty((n, n), dtype=intp)
3034 for l in range(1, n):
3035 for i in range(n - l):
3036 j = i + l
3037 m[i, j] = inf
3038 for k in range(i, j):
3039 q = m[i, k] + m[k + 1, j] + p[i] * p[k + 1] * p[j + 1]
3040 if q < m[i, j]:
3041 m[i, j] = q
3042 s[i, j] = k # Note that Cormen uses 1-based index
3044 return (s, m) if return_costs else s
3047def _multi_dot(arrays, order, i, j, out=None):
3048 """Actually do the multiplication with the given order."""
3049 if i == j:
3050 # the initial call with non-None out should never get here
3051 assert out is None
3053 return arrays[i]
3054 else:
3055 return dot(_multi_dot(arrays, order, i, order[i, j]),
3056 _multi_dot(arrays, order, order[i, j] + 1, j),
3057 out=out)
3060# diagonal
3062def _diagonal_dispatcher(x, /, *, offset=None):
3063 return (x,)
3066@array_function_dispatch(_diagonal_dispatcher)
3067def diagonal(x, /, *, offset=0):
3068 """
3069 Returns specified diagonals of a matrix (or a stack of matrices) ``x``.
3071 This function is Array API compatible, contrary to
3072 :py:func:`numpy.diagonal`, the matrix is assumed
3073 to be defined by the last two dimensions.
3075 Parameters
3076 ----------
3077 x : (...,M,N) array_like
3078 Input array having shape (..., M, N) and whose innermost two
3079 dimensions form MxN matrices.
3080 offset : int, optional
3081 Offset specifying the off-diagonal relative to the main diagonal,
3082 where::
3084 * offset = 0: the main diagonal.
3085 * offset > 0: off-diagonal above the main diagonal.
3086 * offset < 0: off-diagonal below the main diagonal.
3088 Returns
3089 -------
3090 out : (...,min(N,M)) ndarray
3091 An array containing the diagonals and whose shape is determined by
3092 removing the last two dimensions and appending a dimension equal to
3093 the size of the resulting diagonals. The returned array must have
3094 the same data type as ``x``.
3096 See Also
3097 --------
3098 numpy.diagonal
3100 Examples
3101 --------
3102 >>> a = np.arange(4).reshape(2, 2); a
3103 array([[0, 1],
3104 [2, 3]])
3105 >>> np.linalg.diagonal(a)
3106 array([0, 3])
3108 A 3-D example:
3110 >>> a = np.arange(8).reshape(2, 2, 2); a
3111 array([[[0, 1],
3112 [2, 3]],
3113 [[4, 5],
3114 [6, 7]]])
3115 >>> np.linalg.diagonal(a)
3116 array([[0, 3],
3117 [4, 7]])
3119 Diagonals adjacent to the main diagonal can be obtained by using the
3120 `offset` argument:
3122 >>> a = np.arange(9).reshape(3, 3)
3123 >>> a
3124 array([[0, 1, 2],
3125 [3, 4, 5],
3126 [6, 7, 8]])
3127 >>> np.linalg.diagonal(a, offset=1) # First superdiagonal
3128 array([1, 5])
3129 >>> np.linalg.diagonal(a, offset=2) # Second superdiagonal
3130 array([2])
3131 >>> np.linalg.diagonal(a, offset=-1) # First subdiagonal
3132 array([3, 7])
3133 >>> np.linalg.diagonal(a, offset=-2) # Second subdiagonal
3134 array([6])
3136 The anti-diagonal can be obtained by reversing the order of elements
3137 using either `numpy.flipud` or `numpy.fliplr`.
3139 >>> a = np.arange(9).reshape(3, 3)
3140 >>> a
3141 array([[0, 1, 2],
3142 [3, 4, 5],
3143 [6, 7, 8]])
3144 >>> np.linalg.diagonal(np.fliplr(a)) # Horizontal flip
3145 array([2, 4, 6])
3146 >>> np.linalg.diagonal(np.flipud(a)) # Vertical flip
3147 array([6, 4, 2])
3149 Note that the order in which the diagonal is retrieved varies depending
3150 on the flip function.
3152 """
3153 return _core_diagonal(x, offset, axis1=-2, axis2=-1)
3156# trace
3158def _trace_dispatcher(x, /, *, offset=None, dtype=None):
3159 return (x,)
3162@array_function_dispatch(_trace_dispatcher)
3163def trace(x, /, *, offset=0, dtype=None):
3164 """
3165 Returns the sum along the specified diagonals of a matrix
3166 (or a stack of matrices) ``x``.
3168 This function is Array API compatible, contrary to
3169 :py:func:`numpy.trace`.
3171 Parameters
3172 ----------
3173 x : (...,M,N) array_like
3174 Input array having shape (..., M, N) and whose innermost two
3175 dimensions form MxN matrices.
3176 offset : int, optional
3177 Offset specifying the off-diagonal relative to the main diagonal,
3178 where::
3180 * offset = 0: the main diagonal.
3181 * offset > 0: off-diagonal above the main diagonal.
3182 * offset < 0: off-diagonal below the main diagonal.
3184 dtype : dtype, optional
3185 Data type of the returned array.
3187 Returns
3188 -------
3189 out : ndarray
3190 An array containing the traces and whose shape is determined by
3191 removing the last two dimensions and storing the traces in the last
3192 array dimension. For example, if x has rank k and shape:
3193 (I, J, K, ..., L, M, N), then an output array has rank k-2 and shape:
3194 (I, J, K, ..., L) where::
3196 out[i, j, k, ..., l] = trace(a[i, j, k, ..., l, :, :])
3198 The returned array must have a data type as described by the dtype
3199 parameter above.
3201 See Also
3202 --------
3203 numpy.trace
3205 Examples
3206 --------
3207 >>> np.linalg.trace(np.eye(3))
3208 3.0
3209 >>> a = np.arange(8).reshape((2, 2, 2))
3210 >>> np.linalg.trace(a)
3211 array([3, 11])
3213 Trace is computed with the last two axes as the 2-d sub-arrays.
3214 This behavior differs from :py:func:`numpy.trace` which uses the first two
3215 axes by default.
3217 >>> a = np.arange(24).reshape((3, 2, 2, 2))
3218 >>> np.linalg.trace(a).shape
3219 (3, 2)
3221 Traces adjacent to the main diagonal can be obtained by using the
3222 `offset` argument:
3224 >>> a = np.arange(9).reshape((3, 3)); a
3225 array([[0, 1, 2],
3226 [3, 4, 5],
3227 [6, 7, 8]])
3228 >>> np.linalg.trace(a, offset=1) # First superdiagonal
3229 6
3230 >>> np.linalg.trace(a, offset=2) # Second superdiagonal
3231 2
3232 >>> np.linalg.trace(a, offset=-1) # First subdiagonal
3233 10
3234 >>> np.linalg.trace(a, offset=-2) # Second subdiagonal
3235 6
3237 """
3238 return _core_trace(x, offset, axis1=-2, axis2=-1, dtype=dtype)
3241# cross
3243def _cross_dispatcher(x1, x2, /, *, axis=None):
3244 return (x1, x2,)
3247@array_function_dispatch(_cross_dispatcher)
3248def cross(x1, x2, /, *, axis=-1):
3249 """
3250 Returns the cross product of 3-element vectors.
3252 If ``x1`` and/or ``x2`` are multi-dimensional arrays, then
3253 the cross-product of each pair of corresponding 3-element vectors
3254 is independently computed.
3256 This function is Array API compatible, contrary to
3257 :func:`numpy.cross`.
3259 Parameters
3260 ----------
3261 x1 : array_like
3262 The first input array.
3263 x2 : array_like
3264 The second input array. Must be compatible with ``x1`` for all
3265 non-compute axes. The size of the axis over which to compute
3266 the cross-product must be the same size as the respective axis
3267 in ``x1``.
3268 axis : int, optional
3269 The axis (dimension) of ``x1`` and ``x2`` containing the vectors for
3270 which to compute the cross-product. Default: ``-1``.
3272 Returns
3273 -------
3274 out : ndarray
3275 An array containing the cross products.
3277 See Also
3278 --------
3279 numpy.cross
3281 Examples
3282 --------
3283 Vector cross-product.
3285 >>> x = np.array([1, 2, 3])
3286 >>> y = np.array([4, 5, 6])
3287 >>> np.linalg.cross(x, y)
3288 array([-3, 6, -3])
3290 Multiple vector cross-products. Note that the direction of the cross
3291 product vector is defined by the *right-hand rule*.
3293 >>> x = np.array([[1,2,3], [4,5,6]])
3294 >>> y = np.array([[4,5,6], [1,2,3]])
3295 >>> np.linalg.cross(x, y)
3296 array([[-3, 6, -3],
3297 [ 3, -6, 3]])
3299 >>> x = np.array([[1, 2], [3, 4], [5, 6]])
3300 >>> y = np.array([[4, 5], [6, 1], [2, 3]])
3301 >>> np.linalg.cross(x, y, axis=0)
3302 array([[-24, 6],
3303 [ 18, 24],
3304 [-6, -18]])
3306 """
3307 x1 = asanyarray(x1)
3308 x2 = asanyarray(x2)
3310 if x1.shape[axis] != 3 or x2.shape[axis] != 3:
3311 raise ValueError(
3312 "Both input arrays must be (arrays of) 3-dimensional vectors, "
3313 f"but they are {x1.shape[axis]} and {x2.shape[axis]} "
3314 "dimensional instead."
3315 )
3317 return _core_cross(x1, x2, axis=axis)
3320# matmul
3322def _matmul_dispatcher(x1, x2, /):
3323 return (x1, x2)
3326@array_function_dispatch(_matmul_dispatcher)
3327def matmul(x1, x2, /):
3328 """
3329 Computes the matrix product.
3331 This function is Array API compatible, contrary to
3332 :func:`numpy.matmul`.
3334 Parameters
3335 ----------
3336 x1 : array_like
3337 The first input array.
3338 x2 : array_like
3339 The second input array.
3341 Returns
3342 -------
3343 out : ndarray
3344 The matrix product of the inputs.
3345 This is a scalar only when both ``x1``, ``x2`` are 1-d vectors.
3347 Raises
3348 ------
3349 ValueError
3350 If the last dimension of ``x1`` is not the same size as
3351 the second-to-last dimension of ``x2``.
3353 If a scalar value is passed in.
3355 See Also
3356 --------
3357 numpy.matmul
3359 Examples
3360 --------
3361 For 2-D arrays it is the matrix product:
3363 >>> a = np.array([[1, 0],
3364 ... [0, 1]])
3365 >>> b = np.array([[4, 1],
3366 ... [2, 2]])
3367 >>> np.linalg.matmul(a, b)
3368 array([[4, 1],
3369 [2, 2]])
3371 For 2-D mixed with 1-D, the result is the usual.
3373 >>> a = np.array([[1, 0],
3374 ... [0, 1]])
3375 >>> b = np.array([1, 2])
3376 >>> np.linalg.matmul(a, b)
3377 array([1, 2])
3378 >>> np.linalg.matmul(b, a)
3379 array([1, 2])
3382 Broadcasting is conventional for stacks of arrays
3384 >>> a = np.arange(2 * 2 * 4).reshape((2, 2, 4))
3385 >>> b = np.arange(2 * 2 * 4).reshape((2, 4, 2))
3386 >>> np.linalg.matmul(a,b).shape
3387 (2, 2, 2)
3388 >>> np.linalg.matmul(a, b)[0, 1, 1]
3389 98
3390 >>> sum(a[0, 1, :] * b[0 , :, 1])
3391 98
3393 Vector, vector returns the scalar inner product, but neither argument
3394 is complex-conjugated:
3396 >>> np.linalg.matmul([2j, 3j], [2j, 3j])
3397 (-13+0j)
3399 Scalar multiplication raises an error.
3401 >>> np.linalg.matmul([1,2], 3)
3402 Traceback (most recent call last):
3403 ...
3404 ValueError: matmul: Input operand 1 does not have enough dimensions ...
3406 """
3407 return _core_matmul(x1, x2)
3410# tensordot
3412def _tensordot_dispatcher(x1, x2, /, *, axes=None):
3413 return (x1, x2)
3416@array_function_dispatch(_tensordot_dispatcher)
3417def tensordot(x1, x2, /, *, axes=2):
3418 return _core_tensordot(x1, x2, axes=axes)
3421tensordot.__doc__ = _core_tensordot.__doc__
3424# matrix_transpose
3426def _matrix_transpose_dispatcher(x):
3427 return (x,)
3429@array_function_dispatch(_matrix_transpose_dispatcher)
3430def matrix_transpose(x, /):
3431 return _core_matrix_transpose(x)
3434matrix_transpose.__doc__ = f"""{_core_matrix_transpose.__doc__}
3436 Notes
3437 -----
3438 This function is an alias of `numpy.matrix_transpose`.
3439"""
3442# matrix_norm
3444def _matrix_norm_dispatcher(x, /, *, keepdims=None, ord=None):
3445 return (x,)
3447@array_function_dispatch(_matrix_norm_dispatcher)
3448def matrix_norm(x, /, *, keepdims=False, ord="fro"):
3449 """
3450 Computes the matrix norm of a matrix (or a stack of matrices) ``x``.
3452 This function is Array API compatible.
3454 Parameters
3455 ----------
3456 x : array_like
3457 Input array having shape (..., M, N) and whose two innermost
3458 dimensions form ``MxN`` matrices.
3459 keepdims : bool, optional
3460 If this is set to True, the axes which are normed over are left in
3461 the result as dimensions with size one. Default: False.
3462 ord : {1, -1, 2, -2, inf, -inf, 'fro', 'nuc'}, optional
3463 The order of the norm. For details see the table under ``Notes``
3464 in `numpy.linalg.norm`.
3466 See Also
3467 --------
3468 numpy.linalg.norm : Generic norm function
3470 Examples
3471 --------
3472 >>> from numpy import linalg as LA
3473 >>> a = np.arange(9) - 4
3474 >>> a
3475 array([-4, -3, -2, ..., 2, 3, 4])
3476 >>> b = a.reshape((3, 3))
3477 >>> b
3478 array([[-4, -3, -2],
3479 [-1, 0, 1],
3480 [ 2, 3, 4]])
3482 >>> LA.matrix_norm(b)
3483 7.745966692414834
3484 >>> LA.matrix_norm(b, ord='fro')
3485 7.745966692414834
3486 >>> LA.matrix_norm(b, ord=np.inf)
3487 9.0
3488 >>> LA.matrix_norm(b, ord=-np.inf)
3489 2.0
3491 >>> LA.matrix_norm(b, ord=1)
3492 7.0
3493 >>> LA.matrix_norm(b, ord=-1)
3494 6.0
3495 >>> LA.matrix_norm(b, ord=2)
3496 7.3484692283495345
3497 >>> LA.matrix_norm(b, ord=-2)
3498 1.8570331885190563e-016 # may vary
3500 """
3501 x = asanyarray(x)
3502 return norm(x, axis=(-2, -1), keepdims=keepdims, ord=ord)
3505# vector_norm
3507def _vector_norm_dispatcher(x, /, *, axis=None, keepdims=None, ord=None):
3508 return (x,)
3510@array_function_dispatch(_vector_norm_dispatcher)
3511def vector_norm(x, /, *, axis=None, keepdims=False, ord=2):
3512 """
3513 Computes the vector norm of a vector (or batch of vectors) ``x``.
3515 This function is Array API compatible.
3517 Parameters
3518 ----------
3519 x : array_like
3520 Input array.
3521 axis : {None, int, 2-tuple of ints}, optional
3522 If an integer, ``axis`` specifies the axis (dimension) along which
3523 to compute vector norms. If an n-tuple, ``axis`` specifies the axes
3524 (dimensions) along which to compute batched vector norms. If ``None``,
3525 the vector norm must be computed over all array values (i.e.,
3526 equivalent to computing the vector norm of a flattened array).
3527 Default: ``None``.
3528 keepdims : bool, optional
3529 If this is set to True, the axes which are normed over are left in
3530 the result as dimensions with size one. Default: False.
3531 ord : {int, float, inf, -inf}, optional
3532 The order of the norm. For details see the table under ``Notes``
3533 in `numpy.linalg.norm`.
3535 See Also
3536 --------
3537 numpy.linalg.norm : Generic norm function
3539 Examples
3540 --------
3541 >>> from numpy import linalg as LA
3542 >>> a = np.arange(9) + 1
3543 >>> a
3544 array([1, 2, 3, 4, 5, 6, 7, 8, 9])
3545 >>> b = a.reshape((3, 3))
3546 >>> b
3547 array([[1, 2, 3],
3548 [4, 5, 6],
3549 [7, 8, 9]])
3551 >>> LA.vector_norm(b)
3552 16.881943016134134
3553 >>> LA.vector_norm(b, ord=np.inf)
3554 9.0
3555 >>> LA.vector_norm(b, ord=-np.inf)
3556 1.0
3558 >>> LA.vector_norm(b, ord=0)
3559 9.0
3560 >>> LA.vector_norm(b, ord=1)
3561 45.0
3562 >>> LA.vector_norm(b, ord=-1)
3563 0.3534857623790153
3564 >>> LA.vector_norm(b, ord=2)
3565 16.881943016134134
3566 >>> LA.vector_norm(b, ord=-2)
3567 0.8058837395885292
3569 """
3570 x = asanyarray(x)
3571 shape = list(x.shape)
3572 if axis is None:
3573 # Note: np.linalg.norm() doesn't handle 0-D arrays
3574 x = x.ravel()
3575 _axis = 0
3576 elif isinstance(axis, tuple):
3577 # Note: The axis argument supports any number of axes, whereas
3578 # np.linalg.norm() only supports a single axis for vector norm.
3579 normalized_axis = normalize_axis_tuple(axis, x.ndim)
3580 rest = tuple(i for i in range(x.ndim) if i not in normalized_axis)
3581 newshape = axis + rest
3582 x = _core_transpose(x, newshape).reshape(
3583 (
3584 prod([x.shape[i] for i in axis], dtype=int),
3585 *[x.shape[i] for i in rest]
3586 )
3587 )
3588 _axis = 0
3589 else:
3590 _axis = axis
3592 res = norm(x, axis=_axis, ord=ord)
3594 if keepdims:
3595 # We can't reuse np.linalg.norm(keepdims) because of the reshape hacks
3596 # above to avoid matrix norm logic.
3597 _axis = normalize_axis_tuple(
3598 range(len(shape)) if axis is None else axis, len(shape)
3599 )
3600 for i in _axis:
3601 shape[i] = 1
3602 res = res.reshape(tuple(shape))
3604 return res
3607# vecdot
3609def _vecdot_dispatcher(x1, x2, /, *, axis=None):
3610 return (x1, x2)
3612@array_function_dispatch(_vecdot_dispatcher)
3613def vecdot(x1, x2, /, *, axis=-1):
3614 """
3615 Computes the vector dot product.
3617 This function is restricted to arguments compatible with the Array API,
3618 contrary to :func:`numpy.vecdot`.
3620 Let :math:`\\mathbf{a}` be a vector in ``x1`` and :math:`\\mathbf{b}` be
3621 a corresponding vector in ``x2``. The dot product is defined as:
3623 .. math::
3624 \\mathbf{a} \\cdot \\mathbf{b} = \\sum_{i=0}^{n-1} \\overline{a_i}b_i
3626 over the dimension specified by ``axis`` and where :math:`\\overline{a_i}`
3627 denotes the complex conjugate if :math:`a_i` is complex and the identity
3628 otherwise.
3630 Parameters
3631 ----------
3632 x1 : array_like
3633 First input array.
3634 x2 : array_like
3635 Second input array.
3636 axis : int, optional
3637 Axis over which to compute the dot product. Default: ``-1``.
3639 Returns
3640 -------
3641 output : ndarray
3642 The vector dot product of the input.
3644 See Also
3645 --------
3646 numpy.vecdot
3648 Examples
3649 --------
3650 Get the projected size along a given normal for an array of vectors.
3652 >>> v = np.array([[0., 5., 0.], [0., 0., 10.], [0., 6., 8.]])
3653 >>> n = np.array([0., 0.6, 0.8])
3654 >>> np.linalg.vecdot(v, n)
3655 array([ 3., 8., 10.])
3657 """
3658 return _core_vecdot(x1, x2, axis=axis)