summaryrefslogtreecommitdiff
path: root/numpy/linalg/linalg.py
diff options
context:
space:
mode:
Diffstat (limited to 'numpy/linalg/linalg.py')
-rw-r--r--numpy/linalg/linalg.py222
1 files changed, 138 insertions, 84 deletions
diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py
index 816a200eb..eac6267c6 100644
--- a/numpy/linalg/linalg.py
+++ b/numpy/linalg/linalg.py
@@ -8,8 +8,6 @@ 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.
"""
-from __future__ import division, absolute_import, print_function
-
__all__ = ['matrix_power', 'solve', 'tensorsolve', 'tensorinv', 'inv',
'cholesky', 'eigvals', 'eigvalsh', 'pinv', 'slogdet', 'det',
@@ -26,7 +24,7 @@ from numpy.core import (
add, multiply, sqrt, fastCopyAndTranspose, sum, isfinite,
finfo, errstate, geterrobj, moveaxis, amin, amax, product, abs,
atleast_2d, intp, asanyarray, object_, matmul,
- swapaxes, divide, count_nonzero, isnan, sign
+ swapaxes, divide, count_nonzero, isnan, sign, argsort, sort
)
from numpy.core.multiarray import normalize_axis_index
from numpy.core.overrides import set_module
@@ -39,13 +37,6 @@ array_function_dispatch = functools.partial(
overrides.array_function_dispatch, module='numpy.linalg')
-# For Python2/3 compatibility
-_N = b'N'
-_V = b'V'
-_A = b'A'
-_S = b'S'
-_L = b'L'
-
fortran_int = intc
@@ -194,37 +185,33 @@ def _fastCopyAndTranspose(type, *arrays):
else:
return cast_arrays
-def _assertRank2(*arrays):
+def _assert_2d(*arrays):
for a in arrays:
if a.ndim != 2:
raise LinAlgError('%d-dimensional array given. Array must be '
'two-dimensional' % a.ndim)
-def _assertRankAtLeast2(*arrays):
+def _assert_stacked_2d(*arrays):
for a in arrays:
if a.ndim < 2:
raise LinAlgError('%d-dimensional array given. Array must be '
'at least two-dimensional' % a.ndim)
-def _assertNdSquareness(*arrays):
+def _assert_stacked_square(*arrays):
for a in arrays:
m, n = a.shape[-2:]
if m != n:
raise LinAlgError('Last 2 dimensions of the array must be square')
-def _assertFinite(*arrays):
+def _assert_finite(*arrays):
for a in arrays:
- if not (isfinite(a).all()):
+ if not isfinite(a).all():
raise LinAlgError("Array must not contain infs or NaNs")
-def _isEmpty2d(arr):
+def _is_empty_2d(arr):
# check size first for efficiency
return arr.size == 0 and product(arr.shape[-2:]) == 0
-def _assertNoEmpty2d(*arrays):
- for a in arrays:
- if _isEmpty2d(a):
- raise LinAlgError("Arrays cannot be empty")
def transpose(a):
"""
@@ -349,6 +336,10 @@ def solve(a, b):
LinAlgError
If `a` is singular or not square.
+ See Also
+ --------
+ scipy.linalg.solve : Similar function in SciPy.
+
Notes
-----
@@ -386,8 +377,8 @@ def solve(a, b):
"""
a, _ = _makearray(a)
- _assertRankAtLeast2(a)
- _assertNdSquareness(a)
+ _assert_stacked_2d(a)
+ _assert_stacked_square(a)
b, wrap = _makearray(b)
t, result_t = _commonType(a, b)
@@ -506,6 +497,10 @@ def inv(a):
LinAlgError
If `a` is not square or inversion fails.
+ See Also
+ --------
+ scipy.linalg.inv : Similar function in SciPy.
+
Notes
-----
@@ -542,8 +537,8 @@ def inv(a):
"""
a, wrap = _makearray(a)
- _assertRankAtLeast2(a)
- _assertNdSquareness(a)
+ _assert_stacked_2d(a)
+ _assert_stacked_square(a)
t, result_t = _commonType(a)
signature = 'D->D' if isComplexType(t) else 'd->d'
@@ -622,8 +617,8 @@ def matrix_power(a, n):
"""
a = asanyarray(a)
- _assertRankAtLeast2(a)
- _assertNdSquareness(a)
+ _assert_stacked_2d(a)
+ _assert_stacked_square(a)
try:
n = operator.index(n)
@@ -683,8 +678,10 @@ def cholesky(a):
Return the Cholesky decomposition, `L * L.H`, of the square matrix `a`,
where `L` is lower-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. Only `L` is
- actually returned.
+ 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-triangular and diagonal elements of `a`
+ are used. Only `L` is actually returned.
Parameters
----------
@@ -704,6 +701,14 @@ def cholesky(a):
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
-----
@@ -752,15 +757,15 @@ def cholesky(a):
extobj = get_linalg_error_extobj(_raise_linalgerror_nonposdef)
gufunc = _umath_linalg.cholesky_lo
a, wrap = _makearray(a)
- _assertRankAtLeast2(a)
- _assertNdSquareness(a)
+ _assert_stacked_2d(a)
+ _assert_stacked_square(a)
t, result_t = _commonType(a)
signature = 'D->D' if isComplexType(t) else 'd->d'
r = gufunc(a, signature=signature, extobj=extobj)
return wrap(r.astype(result_t, copy=False))
-# QR decompostion
+# QR decomposition
def _qr_dispatcher(a, mode=None):
return (a,)
@@ -816,6 +821,11 @@ def qr(a, mode='reduced'):
LinAlgError
If factoring fails.
+ See Also
+ --------
+ scipy.linalg.qr : Similar function in SciPy.
+ scipy.linalg.rq : Compute RQ decomposition of a matrix.
+
Notes
-----
This is an interface to the LAPACK routines ``dgeqrf``, ``zgeqrf``,
@@ -895,7 +905,7 @@ def qr(a, mode='reduced'):
raise ValueError("Unrecognized mode '%s'" % mode)
a, wrap = _makearray(a)
- _assertRank2(a)
+ _assert_2d(a)
m, n = a.shape
t, result_t = _commonType(a)
a = _fastCopyAndTranspose(t, a)
@@ -1008,6 +1018,7 @@ def eigvals(a):
(conjugate symmetric) arrays.
eigh : eigenvalues and eigenvectors of real symmetric or complex
Hermitian (conjugate symmetric) arrays.
+ scipy.linalg.eigvals : Similar function in SciPy.
Notes
-----
@@ -1047,9 +1058,9 @@ def eigvals(a):
"""
a, wrap = _makearray(a)
- _assertRankAtLeast2(a)
- _assertNdSquareness(a)
- _assertFinite(a)
+ _assert_stacked_2d(a)
+ _assert_stacked_square(a)
+ _assert_finite(a)
t, result_t = _commonType(a)
extobj = get_linalg_error_extobj(
@@ -1109,6 +1120,7 @@ def eigvalsh(a, UPLO='L'):
eigvals : eigenvalues of general real or complex arrays.
eig : eigenvalues and right eigenvectors of general real or complex
arrays.
+ scipy.linalg.eigvalsh : Similar function in SciPy.
Notes
-----
@@ -1157,8 +1169,8 @@ def eigvalsh(a, UPLO='L'):
gufunc = _umath_linalg.eigvalsh_up
a, wrap = _makearray(a)
- _assertRankAtLeast2(a)
- _assertNdSquareness(a)
+ _assert_stacked_2d(a)
+ _assert_stacked_square(a)
t, result_t = _commonType(a)
signature = 'D->d' if isComplexType(t) else 'd->d'
w = gufunc(a, signature=signature, extobj=extobj)
@@ -1207,12 +1219,14 @@ def eig(a):
See Also
--------
eigvals : eigenvalues of a non-symmetric array.
-
eigh : eigenvalues and eigenvectors of a real symmetric or complex
Hermitian (conjugate symmetric) array.
-
eigvalsh : eigenvalues of a real symmetric or complex Hermitian
(conjugate symmetric) array.
+ scipy.linalg.eig : Similar function in SciPy that also solves the
+ generalized eigenvalue problem.
+ scipy.linalg.schur : Best choice for unitary and other non-Hermitian
+ normal matrices.
Notes
-----
@@ -1226,21 +1240,26 @@ def eig(a):
the eigenvalues and eigenvectors of general square arrays.
The number `w` is an eigenvalue of `a` if there exists a vector
- `v` such that ``dot(a,v) = w * v``. Thus, the arrays `a`, `w`, and
- `v` satisfy the equations ``dot(a[:,:], v[:,i]) = w[i] * v[:,i]``
+ `v` such that ``a @ v = w * v``. Thus, the arrays `a`, `w`, and
+ `v` satisfy the equations ``a @ v[:,i] = w[i] * v[:,i]``
for :math:`i \\in \\{0,...,M-1\\}`.
The array `v` of eigenvectors may not be of maximum rank, that is, some
of the columns may be linearly dependent, although round-off error may
obscure that fact. If the eigenvalues are all different, then theoretically
- the eigenvectors are linearly independent. Likewise, the (complex-valued)
- matrix of eigenvectors `v` is unitary if the matrix `a` is normal, i.e.,
- if ``dot(a, a.H) = dot(a.H, a)``, where `a.H` denotes the conjugate
- transpose of `a`.
+ the eigenvectors are linearly independent and `a` can be diagonalized by
+ a similarity transformation using `v`, i.e, ``inv(v) @ a @ v`` is diagonal.
+
+ For non-Hermitian normal matrices the SciPy function `scipy.linalg.schur`
+ is preferred because the matrix `v` is guaranteed to be unitary, which is
+ not the case when using `eig`. The Schur factorization produces an
+ upper triangular matrix rather than a diagonal matrix, but for normal
+ matrices only the diagonal of the upper triangular matrix is needed, the
+ rest is roundoff error.
Finally, it is emphasized that `v` consists of the *right* (as in
right-hand side) eigenvectors of `a`. A vector `y` satisfying
- ``dot(y.T, a) = z * y.T`` for some number `z` is called a *left*
+ ``y.T @ a = z * y.T`` for some number `z` is called a *left*
eigenvector of `a`, and, in general, the left and right eigenvectors
of a matrix are not necessarily the (perhaps conjugate) transposes
of each other.
@@ -1294,9 +1313,9 @@ def eig(a):
"""
a, wrap = _makearray(a)
- _assertRankAtLeast2(a)
- _assertNdSquareness(a)
- _assertFinite(a)
+ _assert_stacked_2d(a)
+ _assert_stacked_square(a)
+ _assert_finite(a)
t, result_t = _commonType(a)
extobj = get_linalg_error_extobj(
@@ -1359,6 +1378,8 @@ def eigh(a, UPLO='L'):
(conjugate symmetric) arrays.
eig : eigenvalues and right eigenvectors for non-symmetric arrays.
eigvals : eigenvalues of non-symmetric arrays.
+ scipy.linalg.eigh : Similar function in SciPy (but also solves the
+ generalized eigenvalue problem).
Notes
-----
@@ -1435,8 +1456,8 @@ def eigh(a, UPLO='L'):
raise ValueError("UPLO argument must be 'L' or 'U'")
a, wrap = _makearray(a)
- _assertRankAtLeast2(a)
- _assertNdSquareness(a)
+ _assert_stacked_2d(a)
+ _assert_stacked_square(a)
t, result_t = _commonType(a)
extobj = get_linalg_error_extobj(
@@ -1510,6 +1531,11 @@ def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
LinAlgError
If SVD computation does not converge.
+ See Also
+ --------
+ scipy.linalg.svd : Similar function in SciPy.
+ scipy.linalg.svdvals : Compute singular values of a matrix.
+
Notes
-----
@@ -1589,26 +1615,31 @@ def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
True
"""
+ import numpy as _nx
a, wrap = _makearray(a)
if hermitian:
- # note: lapack returns eigenvalues in reverse order to our contract.
- # reversing is cheap by design in numpy, so we do so to be consistent
+ # note: lapack svd returns eigenvalues with s ** 2 sorted descending,
+ # but eig returns s sorted ascending, so we re-order the eigenvalues
+ # and related arrays to have the correct order
if compute_uv:
s, u = eigh(a)
- s = s[..., ::-1]
- u = u[..., ::-1]
- # singular values are unsigned, move the sign into v
- vt = transpose(u * sign(s)[..., None, :]).conjugate()
+ sgn = sign(s)
s = abs(s)
+ sidx = argsort(s)[..., ::-1]
+ sgn = _nx.take_along_axis(sgn, sidx, axis=-1)
+ s = _nx.take_along_axis(s, sidx, axis=-1)
+ u = _nx.take_along_axis(u, sidx[..., None, :], axis=-1)
+ # singular values are unsigned, move the sign into v
+ vt = transpose(u * sgn[..., None, :]).conjugate()
return wrap(u), s, wrap(vt)
else:
s = eigvalsh(a)
s = s[..., ::-1]
s = abs(s)
- return s
+ return sort(s)[..., ::-1]
- _assertRankAtLeast2(a)
+ _assert_stacked_2d(a)
t, result_t = _commonType(a)
extobj = get_linalg_error_extobj(_raise_linalgerror_svd_nonconvergence)
@@ -1729,7 +1760,8 @@ def cond(x, p=None):
"""
x = asarray(x) # in case we have a matrix
- _assertNoEmpty2d(x)
+ if _is_empty_2d(x):
+ raise LinAlgError("cond is not defined on empty arrays")
if p is None or p == 2 or p == -2:
s = svd(x, compute_uv=False)
with errstate(all='ignore'):
@@ -1740,8 +1772,8 @@ def cond(x, p=None):
else:
# Call inv(x) ignoring errors. The result array will
# contain nans in the entries where inversion failed.
- _assertRankAtLeast2(x)
- _assertNdSquareness(x)
+ _assert_stacked_2d(x)
+ _assert_stacked_square(x)
t, result_t = _commonType(x)
signature = 'D->D' if isComplexType(t) else 'd->d'
with errstate(all='ignore'):
@@ -1920,6 +1952,13 @@ def pinv(a, rcond=1e-15, hermitian=False):
LinAlgError
If the SVD computation does not converge.
+ See Also
+ --------
+ scipy.linalg.pinv : Similar function in SciPy.
+ scipy.linalg.pinv2 : Similar function in SciPy (SVD-based).
+ scipy.linalg.pinvh : Compute the (Moore-Penrose) pseudo-inverse of a
+ Hermitian matrix.
+
Notes
-----
The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
@@ -1956,7 +1995,7 @@ def pinv(a, rcond=1e-15, hermitian=False):
"""
a, wrap = _makearray(a)
rcond = asarray(rcond)
- if _isEmpty2d(a):
+ if _is_empty_2d(a):
m, n = a.shape[-2:]
res = empty(a.shape[:-2] + (n, m), dtype=a.dtype)
return wrap(res)
@@ -2052,8 +2091,8 @@ def slogdet(a):
"""
a = asarray(a)
- _assertRankAtLeast2(a)
- _assertNdSquareness(a)
+ _assert_stacked_2d(a)
+ _assert_stacked_square(a)
t, result_t = _commonType(a)
real_t = _realType(result_t)
signature = 'D->Dd' if isComplexType(t) else 'd->dd'
@@ -2082,6 +2121,7 @@ def det(a):
--------
slogdet : Another way to represent the determinant, more suitable
for large matrices where underflow/overflow may occur.
+ scipy.linalg.det : Similar function in SciPy.
Notes
-----
@@ -2112,8 +2152,8 @@ def det(a):
"""
a = asarray(a)
- _assertRankAtLeast2(a)
- _assertNdSquareness(a)
+ _assert_stacked_2d(a)
+ _assert_stacked_square(a)
t, result_t = _commonType(a)
signature = 'D->D' if isComplexType(t) else 'd->d'
r = _umath_linalg.det(a, signature=signature)
@@ -2132,13 +2172,13 @@ def lstsq(a, b, rcond="warn"):
r"""
Return the least-squares solution to a linear matrix equation.
- Solves the equation :math:`a x = b` by computing a vector `x` that
- minimizes the squared Euclidean 2-norm :math:`\| b - a x \|^2_2`.
- The equation may be under-, well-, or over-determined (i.e., the
- number of linearly independent rows of `a` can be less than, equal
- to, or greater than its number of linearly independent columns).
+ Computes the vector x that approximatively solves the equation
+ ``a @ x = b``. The equation may be under-, well-, or over-determined
+ (i.e., the number of linearly independent rows of `a` can be less than,
+ equal to, or greater than its number of linearly independent columns).
If `a` is square and of full rank, then `x` (but for round-off error)
- is the "exact" solution of the equation.
+ is the "exact" solution of the equation. Else, `x` minimizes the
+ Euclidean 2-norm :math:`|| b - a x ||`.
Parameters
----------
@@ -2182,6 +2222,10 @@ def lstsq(a, b, rcond="warn"):
LinAlgError
If computation does not converge.
+ See Also
+ --------
+ scipy.linalg.lstsq : Similar function in SciPy.
+
Notes
-----
If `b` is a matrix, then all array results are returned as matrices.
@@ -2224,7 +2268,7 @@ def lstsq(a, b, rcond="warn"):
is_1d = b.ndim == 1
if is_1d:
b = b[:, newaxis]
- _assertRank2(a, b)
+ _assert_2d(a, b)
m, n = a.shape[-2:]
m2, n_rhs = b.shape[-2:]
if m != m2:
@@ -2328,16 +2372,19 @@ def norm(x, ord=None, axis=None, keepdims=False):
Parameters
----------
x : array_like
- Input array. If `axis` is None, `x` must be 1-D or 2-D.
+ Input array. If `axis` is None, `x` must be 1-D or 2-D, unless `ord`
+ is None. If both `axis` and `ord` are None, the 2-norm of
+ ``x.ravel`` will be returned.
ord : {non-zero int, inf, -inf, 'fro', 'nuc'}, optional
Order of the norm (see table under ``Notes``). inf means numpy's
- `inf` object.
- axis : {int, 2-tuple of ints, None}, optional
+ `inf` object. The default is None.
+ axis : {None, int, 2-tuple of ints}, optional.
If `axis` is an integer, it specifies the axis of `x` along which to
compute the vector norms. If `axis` is a 2-tuple, it specifies the
axes that hold 2-D matrices, and the matrix norms of these matrices
are computed. If `axis` is None then either a vector norm (when `x`
- is 1-D) or a matrix norm (when `x` is 2-D) is returned.
+ is 1-D) or a matrix norm (when `x` is 2-D) is returned. The default
+ is None.
.. versionadded:: 1.8.0
@@ -2353,9 +2400,13 @@ def norm(x, ord=None, axis=None, keepdims=False):
n : float or ndarray
Norm of the matrix or vector(s).
+ See Also
+ --------
+ scipy.linalg.norm : Similar function in SciPy.
+
Notes
-----
- For values of ``ord <= 0``, the result is, strictly speaking, not a
+ For values of ``ord < 1``, the result is, strictly speaking, not a
mathematical 'norm', but it may still be useful for various numerical
purposes.
@@ -2383,6 +2434,9 @@ def norm(x, ord=None, axis=None, keepdims=False):
The nuclear norm is the sum of the singular values.
+ Both the Frobenius and nuclear norm orders are only defined for
+ matrices and raise a ValueError when ``x.ndim != 2``.
+
References
----------
.. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*,
@@ -2505,11 +2559,11 @@ def norm(x, ord=None, axis=None, keepdims=False):
# special case for speedup
s = (x.conj() * x).real
return sqrt(add.reduce(s, axis=axis, keepdims=keepdims))
+ # None of the str-type keywords for ord ('fro', 'nuc')
+ # are valid for vectors
+ elif isinstance(ord, str):
+ raise ValueError(f"Invalid norm order '{ord}' for vectors")
else:
- try:
- ord + 1
- except TypeError:
- raise ValueError("Invalid norm order for vectors.")
absx = abs(x)
absx **= ord
ret = add.reduce(absx, axis=axis, keepdims=keepdims)
@@ -2657,7 +2711,7 @@ def multi_dot(arrays):
arrays[0] = atleast_2d(arrays[0])
if arrays[-1].ndim == 1:
arrays[-1] = atleast_2d(arrays[-1]).T
- _assertRank2(*arrays)
+ _assert_2d(*arrays)
# _multi_dot_three is much faster than _multi_dot_matrix_chain_order
if n == 3: