"""Lite version of scipy.linalg.

Notes

-----

This module is a lite version of the linalg.py module in SciPy which

contains high-level Python interface to the LAPACK library. The lite

version only accesses the following LAPACK functions: dgesv, zgesv,

dgeev, zgeev, dgesdd, zgesdd, dgelsd, zgelsd, dsyevd, zheevd, dgetrf,

zgetrf, dpotrf, zpotrf, dgeqrf, zgeqrf, zungqr, dorgqr.

"""

__all__ = ['matrix_power', 'solve', 'tensorsolve', 'tensorinv', 'inv',

'cholesky', 'eigvals', 'eigvalsh', 'pinv', 'slogdet', 'det',

'svd', 'svdvals', 'eig', 'eigh', 'lstsq', 'norm', 'qr', 'cond',

'matrix_rank', 'LinAlgError', 'multi_dot', 'trace', 'diagonal',

'cross', 'outer', 'tensordot', 'matmul', 'matrix_transpose',

'matrix_norm', 'vector_norm', 'vecdot']

import functools

import operator

import warnings

from typing import Any, NamedTuple

from numpy._core import (

abs,

add,

all,

amax,

amin,

argsort,

array,

asanyarray,

asarray,

atleast_2d,

cdouble,

complexfloating,

count_nonzero,

cross as _core_cross,

csingle,

diagonal as _core_diagonal,

divide,

dot,

double,

empty,

empty_like,

errstate,

finfo,

inexact,

inf,

intc,

intp,

isfinite,

isnan,

matmul as _core_matmul,

matrix_transpose as _core_matrix_transpose,

moveaxis,

multiply,

newaxis,

object_,

outer as _core_outer,

overrides,

prod,

reciprocal,

sign,

single,

sort,

sqrt,

sum,

swapaxes,

tensordot as _core_tensordot,

trace as _core_trace,

transpose as _core_transpose,

vecdot as _core_vecdot,

zeros,

)

from numpy._globals import _NoValue

from numpy._typing import NDArray

from numpy._utils import set_module

from numpy.lib._twodim_base_impl import eye, triu

from numpy.lib.array_utils import normalize_axis_index, normalize_axis_tuple

from numpy.linalg import _umath_linalg

class EigResult(NamedTuple):

eigenvalues: NDArray[Any]

eigenvectors: NDArray[Any]

class EighResult(NamedTuple):

eigenvalues: NDArray[Any]

eigenvectors: NDArray[Any]

class QRResult(NamedTuple):

Q: NDArray[Any]

R: NDArray[Any]

class SlogdetResult(NamedTuple):

sign: NDArray[Any]

logabsdet: NDArray[Any]

class SVDResult(NamedTuple):

U: NDArray[Any]

S: NDArray[Any]

Vh: NDArray[Any]

array_function_dispatch = functools.partial(

overrides.array_function_dispatch, module='numpy.linalg'

)

fortran_int = intc

@set_module('numpy.linalg')

class LinAlgError(ValueError):

"""

Generic Python-exception-derived object raised by linalg functions.

General purpose exception class, derived from Python's ValueError

class, programmatically raised in linalg functions when a Linear

Algebra-related condition would prevent further correct execution of the

function.

Parameters

----------

None

Examples

--------

>>> from numpy import linalg as LA

>>> LA.inv(np.zeros((2,2)))

Traceback (most recent call last):

File "<stdin>", line 1, in <module>

File "...linalg.py", line 350,

in inv return wrap(solve(a, identity(a.shape[0], dtype=a.dtype)))

File "...linalg.py", line 249,

in solve

raise LinAlgError('Singular matrix')

numpy.linalg.LinAlgError: Singular matrix

"""

def _raise_linalgerror_singular(err, flag):

raise LinAlgError("Singular matrix")

def _raise_linalgerror_nonposdef(err, flag):

raise LinAlgError("Matrix is not positive definite")

def _raise_linalgerror_eigenvalues_nonconvergence(err, flag):

raise LinAlgError("Eigenvalues did not converge")

def _raise_linalgerror_svd_nonconvergence(err, flag):

raise LinAlgError("SVD did not converge")

def _raise_linalgerror_lstsq(err, flag):

raise LinAlgError("SVD did not converge in Linear Least Squares")

def _raise_linalgerror_qr(err, flag):

raise LinAlgError("Incorrect argument found while performing "

"QR factorization")

def _makearray(a):

new = asarray(a)

wrap = getattr(a, "__array_wrap__", new.__array_wrap__)

return new, wrap

def isComplexType(t):

return issubclass(t, complexfloating)

_real_types_map = {single: single,

double: double,

csingle: single,

cdouble: double}

_complex_types_map = {single: csingle,

double: cdouble,

csingle: csingle,

cdouble: cdouble}

def _realType(t, default=double):

return _real_types_map.get(t, default)

def _complexType(t, default=cdouble):

return _complex_types_map.get(t, default)

def _to_real_if_imag_zero(w, t):

"""Backwards compat helper: force w to be real if t.dtype is real and w.imag == 0

"""

result_t = t.dtype.type

if not isComplexType(result_t) and all(w.imag == 0.0):

w = w.real

result_t = _realType(result_t)

else:

result_t = _complexType(result_t)

return w.astype(result_t, copy=False)

def _commonType(*arrays):

# in lite version, use higher precision (always double or cdouble)

result_type = single

is_complex = False

for a in arrays:

type_ = a.dtype.type

if issubclass(type_, inexact):

if isComplexType(type_):

is_complex = True

rt = _realType(type_, default=None)

if rt is double:

result_type = double

elif rt is None:

# unsupported inexact scalar

raise TypeError(f"array type {a.dtype.name} is unsupported in linalg")

else:

result_type = double

if is_complex:

result_type = _complex_types_map[result_type]

return cdouble, result_type

else:

return double, result_type

def _to_native_byte_order(*arrays):

ret = []

for arr in arrays:

if arr.dtype.byteorder not in ('=', '|'):

ret.append(asarray(arr, dtype=arr.dtype.newbyteorder('=')))

else:

ret.append(arr)

if len(ret) == 1:

return ret[0]

else:

return ret

def _assert_2d(*arrays):

for a in arrays:

if a.ndim != 2:

raise LinAlgError(f'{a.ndim}-dimensional array given. Array must be '

'two-dimensional')

def _assert_stacked_2d(*arrays):

for a in arrays:

if a.ndim < 2:

raise LinAlgError(f'{a.ndim}-dimensional array given. Array must be '

'at least two-dimensional')

def _assert_stacked_square(*arrays):

for a in arrays:

try:

m, n = a.shape[-2:]

except ValueError:

raise LinAlgError(f'{a.ndim}-dimensional array given. Array must be '

'at least two-dimensional')

if m != n:

raise LinAlgError('Last 2 dimensions of the array must be square')

def _assert_finite(*arrays):

for a in arrays:

if not isfinite(a).all():

raise LinAlgError("Array must not contain infs or NaNs")

def _is_empty_2d(arr):

# check size first for efficiency

return arr.size == 0 and prod(arr.shape[-2:]) == 0

def transpose(a):

"""

Transpose each matrix in a stack of matrices.

Unlike np.transpose, this only swaps the last two axes, rather than all of

them

Parameters

----------

a : (...,M,N) array_like

Returns

-------

aT : (...,N,M) ndarray

"""

return swapaxes(a, -1, -2)

# Linear equations

def _tensorsolve_dispatcher(a, b, axes=None):

return (a, b)

@array_function_dispatch(_tensorsolve_dispatcher)

def tensorsolve(a, b, axes=None):

"""

Solve the tensor equation ``a x = b`` for x.

It is assumed that all indices of `x` are summed over in the product,

together with the rightmost indices of `a`, as is done in, for example,

``tensordot(a, x, axes=x.ndim)``.

Parameters

----------

a : array_like

Coefficient tensor, of shape ``b.shape + Q``. `Q`, a tuple, equals

the shape of that sub-tensor of `a` consisting of the appropriate

number of its rightmost indices, and must be such that

``prod(Q) == prod(b.shape)`` (in which sense `a` is said to be

'square').

b : array_like

Right-hand tensor, which can be of any shape.

axes : tuple of ints, optional

Axes in `a` to reorder to the right, before inversion.

If None (default), no reordering is done.

Returns

-------

x : ndarray, shape Q

Raises

------

LinAlgError

If `a` is singular or not 'square' (in the above sense).

See Also

--------

numpy.tensordot, tensorinv, numpy.einsum

Examples

--------

>>> import numpy as np

>>> a = np.eye(2*3*4).reshape((2*3, 4, 2, 3, 4))

>>> rng = np.random.default_rng()

>>> b = rng.normal(size=(2*3, 4))

>>> x = np.linalg.tensorsolve(a, b)

>>> x.shape

(2, 3, 4)

>>> np.allclose(np.tensordot(a, x, axes=3), b)

True

"""

a, wrap = _makearray(a)

b = asarray(b)

an = a.ndim

if axes is not None:

allaxes = list(range(an))

for k in axes:

allaxes.remove(k)

allaxes.insert(an, k)

a = a.transpose(allaxes)

oldshape = a.shape[-(an - b.ndim):]

prod = 1

for k in oldshape:

prod *= k

if a.size != prod ** 2:

raise LinAlgError(

"Input arrays must satisfy the requirement \

prod(a.shape[b.ndim:]) == prod(a.shape[:b.ndim])"

)

a = a.reshape(prod, prod)

b = b.ravel()

res = wrap(solve(a, b))

return res.reshape(oldshape)

def _solve_dispatcher(a, b):

return (a, b)

@array_function_dispatch(_solve_dispatcher)

def solve(a, b):

"""

Solve a linear matrix equation, or system of linear scalar equations.

Computes the "exact" solution, `x`, of the well-determined, i.e., full

rank, linear matrix equation `ax = b`.

Parameters

----------

a : (..., M, M) array_like

Coefficient matrix.

b : {(M,), (..., M, K)}, array_like

Ordinate or "dependent variable" values.

Returns

-------

x : {(..., M,), (..., M, K)} ndarray

Solution to the system a x = b. Returned shape is (..., M) if b is

shape (M,) and (..., M, K) if b is (..., M, K), where the "..." part is

broadcasted between a and b.

Raises

------

LinAlgError

If `a` is singular or not square.

See Also

--------

scipy.linalg.solve : Similar function in SciPy.

Notes

-----

Broadcasting rules apply, see the `numpy.linalg` documentation for

details.

The solutions are computed using LAPACK routine ``_gesv``.

`a` must be square and of full-rank, i.e., all rows (or, equivalently,

columns) must be linearly independent; if either is not true, use

`lstsq` for the least-squares best "solution" of the

system/equation.

.. versionchanged:: 2.0

The b array is only treated as a shape (M,) column vector if it is

exactly 1-dimensional. In all other instances it is treated as a stack

of (M, K) matrices. Previously b would be treated as a stack of (M,)

vectors if b.ndim was equal to a.ndim - 1.

References

----------

.. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,

FL, Academic Press, Inc., 1980, pg. 22.

Examples

--------

Solve the system of equations:

``x0 + 2 * x1 = 1`` and

``3 * x0 + 5 * x1 = 2``:

>>> import numpy as np

>>> a = np.array([[1, 2], [3, 5]])

>>> b = np.array([1, 2])

>>> x = np.linalg.solve(a, b)

>>> x

array([-1., 1.])

Check that the solution is correct:

>>> np.allclose(np.dot(a, x), b)

True

"""

a, _ = _makearray(a)

_assert_stacked_square(a)

b, wrap = _makearray(b)

t, result_t = _commonType(a, b)

# We use the b = (..., M,) logic, only if the number of extra dimensions

# match exactly

if b.ndim == 1:

gufunc = _umath_linalg.solve1

else:

gufunc = _umath_linalg.solve

signature = 'DD->D' if isComplexType(t) else 'dd->d'

with errstate(call=_raise_linalgerror_singular, invalid='call',

over='ignore', divide='ignore', under='ignore'):

r = gufunc(a, b, signature=signature)

return wrap(r.astype(result_t, copy=False))

def _tensorinv_dispatcher(a, ind=None):

return (a,)

@array_function_dispatch(_tensorinv_dispatcher)

def tensorinv(a, ind=2):

"""

Compute the 'inverse' of an N-dimensional array.

The result is an inverse for `a` relative to the tensordot operation

``tensordot(a, b, ind)``, i. e., up to floating-point accuracy,

``tensordot(tensorinv(a), a, ind)`` is the "identity" tensor for the

tensordot operation.

Parameters

----------

a : array_like

Tensor to 'invert'. Its shape must be 'square', i. e.,

``prod(a.shape[:ind]) == prod(a.shape[ind:])``.

ind : int, optional

Number of first indices that are involved in the inverse sum.

Must be a positive integer, default is 2.

Returns

-------

b : ndarray

`a`'s tensordot inverse, shape ``a.shape[ind:] + a.shape[:ind]``.

Raises

------

LinAlgError

If `a` is singular or not 'square' (in the above sense).

See Also

--------

numpy.tensordot, tensorsolve

Examples

--------

>>> import numpy as np

>>> a = np.eye(4*6).reshape((4, 6, 8, 3))

>>> ainv = np.linalg.tensorinv(a, ind=2)

>>> ainv.shape

(8, 3, 4, 6)

>>> rng = np.random.default_rng()

>>> b = rng.normal(size=(4, 6))

>>> np.allclose(np.tensordot(ainv, b), np.linalg.tensorsolve(a, b))

True

>>> a = np.eye(4*6).reshape((24, 8, 3))

>>> ainv = np.linalg.tensorinv(a, ind=1)

>>> ainv.shape

(8, 3, 24)

>>> rng = np.random.default_rng()

>>> b = rng.normal(size=24)

>>> np.allclose(np.tensordot(ainv, b, 1), np.linalg.tensorsolve(a, b))

True

"""

a = asarray(a)

oldshape = a.shape

prod = 1

if ind > 0:

invshape = oldshape[ind:] + oldshape[:ind]

for k in oldshape[ind:]:

prod *= k

else:

raise ValueError("Invalid ind argument.")

a = a.reshape(prod, -1)

ia = inv(a)

return ia.reshape(*invshape)

# Matrix inversion

def _unary_dispatcher(a):

return (a,)

@array_function_dispatch(_unary_dispatcher)

def inv(a):

"""

Compute the inverse of a matrix.

Given a square matrix `a`, return the matrix `ainv` satisfying

``a @ ainv = ainv @ a = eye(a.shape[0])``.

Parameters

----------

a : (..., M, M) array_like

Matrix to be inverted.

Returns

-------

ainv : (..., M, M) ndarray or matrix

Inverse of the matrix `a`.

Raises

------

LinAlgError

If `a` is not square or inversion fails.

See Also

--------

scipy.linalg.inv : Similar function in SciPy.

numpy.linalg.cond : Compute the condition number of a matrix.

numpy.linalg.svd : Compute the singular value decomposition of a matrix.

Notes

-----

Broadcasting rules apply, see the `numpy.linalg` documentation for

details.

If `a` is detected to be singular, a `LinAlgError` is raised. If `a` is

ill-conditioned, a `LinAlgError` may or may not be raised, and results may

be inaccurate due to floating-point errors.

References

----------

.. [1] Wikipedia, "Condition number",

https://en.wikipedia.org/wiki/Condition_number

Examples

--------

>>> import numpy as np

>>> from numpy.linalg import inv

>>> a = np.array([[1., 2.], [3., 4.]])

>>> ainv = inv(a)

>>> np.allclose(a @ ainv, np.eye(2))

True

>>> np.allclose(ainv @ a, np.eye(2))

True

If a is a matrix object, then the return value is a matrix as well:

>>> ainv = inv(np.matrix(a))

>>> ainv

matrix([[-2. , 1. ],

[ 1.5, -0.5]])

Inverses of several matrices can be computed at once:

>>> a = np.array([[[1., 2.], [3., 4.]], [[1, 3], [3, 5]]])

>>> inv(a)

array([[[-2. , 1. ],

[ 1.5 , -0.5 ]],

[[-1.25, 0.75],

[ 0.75, -0.25]]])

If a matrix is close to singular, the computed inverse may not satisfy

``a @ ainv = ainv @ a = eye(a.shape[0])`` even if a `LinAlgError`

is not raised:

>>> a = np.array([[2,4,6],[2,0,2],[6,8,14]])

>>> inv(a) # No errors raised

array([[-1.12589991e+15, -5.62949953e+14, 5.62949953e+14],

[-1.12589991e+15, -5.62949953e+14, 5.62949953e+14],

[ 1.12589991e+15, 5.62949953e+14, -5.62949953e+14]])

>>> a @ inv(a)

array([[ 0. , -0.5 , 0. ], # may vary

[-0.5 , 0.625, 0.25 ],

[ 0. , 0. , 1. ]])

To detect ill-conditioned matrices, you can use `numpy.linalg.cond` to

compute its *condition number* [1]_. The larger the condition number, the

more ill-conditioned the matrix is. As a rule of thumb, if the condition

number ``cond(a) = 10**k``, then you may lose up to ``k`` digits of

accuracy on top of what would be lost to the numerical method due to loss

of precision from arithmetic methods.

>>> from numpy.linalg import cond

>>> cond(a)

np.float64(8.659885634118668e+17) # may vary

It is also possible to detect ill-conditioning by inspecting the matrix's

singular values directly. The ratio between the largest and the smallest

singular value is the condition number:

>>> from numpy.linalg import svd

>>> sigma = svd(a, compute_uv=False) # Do not compute singular vectors

>>> sigma.max()/sigma.min()

8.659885634118668e+17 # may vary

"""

a, wrap = _makearray(a)

_assert_stacked_square(a)

t, result_t = _commonType(a)

signature = 'D->D' if isComplexType(t) else 'd->d'

with errstate(call=_raise_linalgerror_singular, invalid='call',

over='ignore', divide='ignore', under='ignore'):

ainv = _umath_linalg.inv(a, signature=signature)

return wrap(ainv.astype(result_t, copy=False))

def _matrix_power_dispatcher(a, n):

return (a,)

@array_function_dispatch(_matrix_power_dispatcher)

def matrix_power(a, n):

"""

Raise a square matrix to the (integer) power `n`.

For positive integers `n`, the power is computed by repeated matrix

squarings and matrix multiplications. If ``n == 0``, the identity matrix

of the same shape as M is returned. If ``n < 0``, the inverse

is computed and then raised to the ``abs(n)``.

.. note:: Stacks of object matrices are not currently supported.

Parameters

----------

a : (..., M, M) array_like

Matrix to be "powered".

n : int

The exponent can be any integer or long integer, positive,

negative, or zero.

Returns

-------

a**n : (..., M, M) ndarray or matrix object

The return value is the same shape and type as `M`;

if the exponent is positive or zero then the type of the

elements is the same as those of `M`. If the exponent is

negative the elements are floating-point.

Raises

------

LinAlgError

For matrices that are not square or that (for negative powers) cannot

be inverted numerically.

Examples

--------

>>> import numpy as np

>>> from numpy.linalg import matrix_power

>>> i = np.array([[0, 1], [-1, 0]]) # matrix equiv. of the imaginary unit

>>> matrix_power(i, 3) # should = -i

array([[ 0, -1],

[ 1, 0]])

>>> matrix_power(i, 0)

array([[1, 0],

[0, 1]])

>>> matrix_power(i, -3) # should = 1/(-i) = i, but w/ f.p. elements

array([[ 0., 1.],

[-1., 0.]])

Somewhat more sophisticated example

>>> q = np.zeros((4, 4))

>>> q[0:2, 0:2] = -i

>>> q[2:4, 2:4] = i

>>> q # one of the three quaternion units not equal to 1

array([[ 0., -1., 0., 0.],

[ 1., 0., 0., 0.],

[ 0., 0., 0., 1.],

[ 0., 0., -1., 0.]])

>>> matrix_power(q, 2) # = -np.eye(4)

array([[-1., 0., 0., 0.],

[ 0., -1., 0., 0.],

[ 0., 0., -1., 0.],

[ 0., 0., 0., -1.]])

"""

a = asanyarray(a)

_assert_stacked_square(a)

try:

n = operator.index(n)

except TypeError as e:

raise TypeError("exponent must be an integer") from e

# Fall back on dot for object arrays. Object arrays are not supported by

# the current implementation of matmul using einsum

if a.dtype != object:

fmatmul = matmul

elif a.ndim == 2:

fmatmul = dot

else:

raise NotImplementedError(

"matrix_power not supported for stacks of object arrays")

if n == 0:

a = empty_like(a)

a[...] = eye(a.shape[-2], dtype=a.dtype)

return a

elif n < 0:

a = inv(a)

n = abs(n)

# short-cuts.

if n == 1:

return a

elif n == 2:

return fmatmul(a, a)

elif n == 3:

return fmatmul(fmatmul(a, a), a)

# Use binary decomposition to reduce the number of matrix multiplications.

# Here, we iterate over the bits of n, from LSB to MSB, raise `a` to

# increasing powers of 2, and multiply into the result as needed.

z = result = None

while n > 0:

z = a if z is None else fmatmul(z, z)

n, bit = divmod(n, 2)

if bit:

result = z if result is None else fmatmul(result, z)

return result

# Cholesky decomposition

def _cholesky_dispatcher(a, /, *, upper=None):

return (a,)

@array_function_dispatch(_cholesky_dispatcher)

def cholesky(a, /, *, upper=False):

"""

Cholesky decomposition.

Return the lower or upper Cholesky decomposition, ``L * L.H`` or

``U.H * U``, of the square matrix ``a``, where ``L`` is lower-triangular,

``U`` is upper-triangular, and ``.H`` is the conjugate transpose operator

(which is the ordinary transpose if ``a`` is real-valued). ``a`` must be

Hermitian (symmetric if real-valued) and positive-definite. No checking is

performed to verify whether ``a`` is Hermitian or not. In addition, only

the lower or upper-triangular and diagonal elements of ``a`` are used.

Only ``L`` or ``U`` is actually returned.

Parameters

----------

a : (..., M, M) array_like

Hermitian (symmetric if all elements are real), positive-definite

input matrix.

upper : bool

If ``True``, the result must be the upper-triangular Cholesky factor.

If ``False``, the result must be the lower-triangular Cholesky factor.

Default: ``False``.

Returns

-------

L : (..., M, M) array_like

Lower or upper-triangular Cholesky factor of `a`. Returns a matrix

object if `a` is a matrix object.

Raises

------

LinAlgError

If the decomposition fails, for example, if `a` is not

positive-definite.

See Also

--------

scipy.linalg.cholesky : Similar function in SciPy.

scipy.linalg.cholesky_banded : Cholesky decompose a banded Hermitian

positive-definite matrix.

scipy.linalg.cho_factor : Cholesky decomposition of a matrix, to use in

`scipy.linalg.cho_solve`.

Notes

-----

Broadcasting rules apply, see the `numpy.linalg` documentation for

details.

The Cholesky decomposition is often used as a fast way of solving

.. math:: A \\mathbf{x} = \\mathbf{b}

(when `A` is both Hermitian/symmetric and positive-definite).

First, we solve for :math:`\\mathbf{y}` in

.. math:: L \\mathbf{y} = \\mathbf{b},

and then for :math:`\\mathbf{x}` in

.. math:: L^{H} \\mathbf{x} = \\mathbf{y}.

Examples

--------

>>> import numpy as np

>>> A = np.array([[1,-2j],[2j,5]])

>>> A

array([[ 1.+0.j, -0.-2.j],

[ 0.+2.j, 5.+0.j]])

>>> L = np.linalg.cholesky(A)

>>> L

array([[1.+0.j, 0.+0.j],

[0.+2.j, 1.+0.j]])

>>> np.dot(L, L.T.conj()) # verify that L * L.H = A

array([[1.+0.j, 0.-2.j],

[0.+2.j, 5.+0.j]])

>>> A = [[1,-2j],[2j,5]] # what happens if A is only array_like?

>>> np.linalg.cholesky(A) # an ndarray object is returned

array([[1.+0.j, 0.+0.j],

[0.+2.j, 1.+0.j]])

>>> # But a matrix object is returned if A is a matrix object

>>> np.linalg.cholesky(np.matrix(A))

matrix([[ 1.+0.j, 0.+0.j],

[ 0.+2.j, 1.+0.j]])

>>> # The upper-triangular Cholesky factor can also be obtained.

>>> np.linalg.cholesky(A, upper=True)

array([[1.-0.j, 0.-2.j],

[0.-0.j, 1.-0.j]])

"""

gufunc = _umath_linalg.cholesky_up if upper else _umath_linalg.cholesky_lo

a, wrap = _makearray(a)

_assert_stacked_square(a)

t, result_t = _commonType(a)

signature = 'D->D' if isComplexType(t) else 'd->d'

with errstate(call=_raise_linalgerror_nonposdef, invalid='call',

over='ignore', divide='ignore', under='ignore'):

r = gufunc(a, signature=signature)

return wrap(r.astype(result_t, copy=False))

# outer product

def _outer_dispatcher(x1, x2):

return (x1, x2)

@array_function_dispatch(_outer_dispatcher)

def outer(x1, x2, /):

"""

Compute the outer product of two vectors.

This function is Array API compatible. Compared to ``np.outer``

it accepts 1-dimensional inputs only.

Parameters

----------

x1 : (M,) array_like

One-dimensional input array of size ``N``.

Must have a numeric data type.

x2 : (N,) array_like

One-dimensional input array of size ``M``.

Must have a numeric data type.

Returns

-------

out : (M, N) ndarray

``out[i, j] = a[i] * b[j]``

See also

--------

outer

Examples

--------

Make a (*very* coarse) grid for computing a Mandelbrot set:

>>> rl = np.linalg.outer(np.ones((5,)), np.linspace(-2, 2, 5))

>>> rl

array([[-2., -1., 0., 1., 2.],

[-2., -1., 0., 1., 2.],

[-2., -1., 0., 1., 2.],

[-2., -1., 0., 1., 2.],

[-2., -1., 0., 1., 2.]])

>>> im = np.linalg.outer(1j*np.linspace(2, -2, 5), np.ones((5,)))

>>> im

array([[0.+2.j, 0.+2.j, 0.+2.j, 0.+2.j, 0.+2.j],

[0.+1.j, 0.+1.j, 0.+1.j, 0.+1.j, 0.+1.j],

[0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],

[0.-1.j, 0.-1.j, 0.-1.j, 0.-1.j, 0.-1.j],

[0.-2.j, 0.-2.j, 0.-2.j, 0.-2.j, 0.-2.j]])

>>> grid = rl + im

>>> grid

array([[-2.+2.j, -1.+2.j, 0.+2.j, 1.+2.j, 2.+2.j],

[-2.+1.j, -1.+1.j, 0.+1.j, 1.+1.j, 2.+1.j],

[-2.+0.j, -1.+0.j, 0.+0.j, 1.+0.j, 2.+0.j],

[-2.-1.j, -1.-1.j, 0.-1.j, 1.-1.j, 2.-1.j],

[-2.-2.j, -1.-2.j, 0.-2.j, 1.-2.j, 2.-2.j]])

An example using a "vector" of letters:

>>> x = np.array(['a', 'b', 'c'], dtype=np.object_)

>>> np.linalg.outer(x, [1, 2, 3])

array([['a', 'aa', 'aaa'],

['b', 'bb', 'bbb'],

['c', 'cc', 'ccc']], dtype=object)

"""

x1 = asanyarray(x1)

x2 = asanyarray(x2)

if x1.ndim != 1 or x2.ndim != 1:

raise ValueError(

"Input arrays must be one-dimensional, but they are "

f"{x1.ndim=} and {x2.ndim=}."

)

return _core_outer(x1, x2, out=None)

# QR decomposition

def _qr_dispatcher(a, mode=None):

return (a,)

@array_function_dispatch(_qr_dispatcher)

def qr(a, mode='reduced'):

"""

Compute the qr factorization of a matrix.

Factor the matrix `a` as *qr*, where `q` is orthonormal and `r` is

upper-triangular.

Parameters

----------

a : array_like, shape (..., M, N)

An array-like object with the dimensionality of at least 2.

mode : {'reduced', 'complete', 'r', 'raw'}, optional, default: 'reduced'

If K = min(M, N), then

* 'reduced' : returns Q, R with dimensions (..., M, K), (..., K, N)

* 'complete' : returns Q, R with dimensions (..., M, M), (..., M, N)

* 'r' : returns R only with dimensions (..., K, N)

* 'raw' : returns h, tau with dimensions (..., N, M), (..., K,)

The options 'reduced', 'complete, and 'raw' are new in numpy 1.8,

see the notes for more information. The default is 'reduced', and to

maintain backward compatibility with earlier versions of numpy both

it and the old default 'full' can be omitted. Note that array h