diff options
Diffstat (limited to 'numpy/linalg/linalg.py')
| -rw-r--r-- | numpy/linalg/linalg.py | 633 |
1 files changed, 290 insertions, 343 deletions
diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py index 34fe9b550..ae0da3685 100644 --- a/numpy/linalg/linalg.py +++ b/numpy/linalg/linalg.py @@ -22,9 +22,9 @@ from numpy.core import array, asarray, zeros, empty, transpose, \ intc, single, double, csingle, cdouble, inexact, complexfloating, \ newaxis, ravel, all, Inf, dot, add, multiply, identity, sqrt, \ maximum, flatnonzero, diagonal, arange, fastCopyAndTranspose, sum, \ - isfinite, size, finfo, absolute, log, exp + isfinite, size, finfo, absolute, log, exp, errstate, geterrobj from numpy.lib import triu -from numpy.linalg import lapack_lite +from numpy.linalg import lapack_lite, _umath_linalg from numpy.matrixlib.defmatrix import matrix_power from numpy.compat import asbytes @@ -67,6 +67,40 @@ class LinAlgError(Exception): """ pass +# Dealing with errors in _umath_linalg + +_linalg_error_extobj = None + +def _determine_error_states(): + global _linalg_error_extobj + errobj = geterrobj() + bufsize = errobj[0] + + with errstate(invalid='call', over='ignore', + divide='ignore', under='ignore'): + invalid_call_errmask = geterrobj()[1] + + _linalg_error_extobj = [bufsize, invalid_call_errmask, None] + +_determine_error_states() + +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 get_linalg_error_extobj(callback): + extobj = list(_linalg_error_extobj) + extobj[2] = callback + return extobj + def _makearray(a): new = asarray(a) wrap = getattr(a, "__array_prepare__", new.__array_wrap__) @@ -124,6 +158,7 @@ def _commonType(*arrays): t = double return t, result_type + # _fastCopyAndTranpose assumes the input is 2D (as all the calls in here are). _fastCT = fastCopyAndTranspose @@ -158,11 +193,22 @@ def _assertRank2(*arrays): raise LinAlgError('%d-dimensional array given. Array must be ' 'two-dimensional' % len(a.shape)) +def _assertRankAtLeast2(*arrays): + for a in arrays: + if len(a.shape) < 2: + raise LinAlgError('%d-dimensional array given. Array must be ' + 'at least two-dimensional' % len(a.shape)) + def _assertSquareness(*arrays): for a in arrays: if max(a.shape) != min(a.shape): raise LinAlgError('Array must be square') +def _assertNdSquareness(*arrays): + for a in arrays: + if max(a.shape[-2:]) != min(a.shape[-2:]): + raise LinAlgError('Last 2 dimensions of the array must be square') + def _assertFinite(*arrays): for a in arrays: if not (isfinite(a).all()): @@ -254,14 +300,14 @@ def solve(a, b): Parameters ---------- - a : (M, M) array_like + a : (..., M, M) array_like Coefficient matrix. - b : {(M,), (M, N)}, array_like + b : {(..., M,), (..., M, K)}, array_like Ordinate or "dependent variable" values. Returns ------- - x : {(M,), (M, N)} ndarray + x : {(..., M,), (..., M, K)} ndarray Solution to the system a x = b. Returned shape is identical to `b`. Raises @@ -271,15 +317,10 @@ def solve(a, b): Notes ----- - `solve` is a wrapper for the LAPACK routines `dgesv`_ and - `zgesv`_, the former being used if `a` is real-valued, the latter if - it is complex-valued. The solution to the system of linear equations - is computed using an LU decomposition [1]_ with partial pivoting and - row interchanges. - - .. _dgesv: http://www.netlib.org/lapack/double/dgesv.f + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. - .. _zgesv: http://www.netlib.org/lapack/complex16/zgesv.f + 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 @@ -308,32 +349,22 @@ def solve(a, b): """ a, _ = _makearray(a) + _assertNonEmpty(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) b, wrap = _makearray(b) - one_eq = len(b.shape) == 1 - if one_eq: - b = b[:, newaxis] - _assertRank2(a, b) - _assertSquareness(a) - n_eq = a.shape[0] - n_rhs = b.shape[1] - if n_eq != b.shape[0]: - raise LinAlgError('Incompatible dimensions') t, result_t = _commonType(a, b) -# lapack_routine = _findLapackRoutine('gesv', t) - if isComplexType(t): - lapack_routine = lapack_lite.zgesv - else: - lapack_routine = lapack_lite.dgesv - a, b = _fastCopyAndTranspose(t, a, b) - a, b = _to_native_byte_order(a, b) - pivots = zeros(n_eq, fortran_int) - results = lapack_routine(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0) - if results['info'] > 0: - raise LinAlgError('Singular matrix') - if one_eq: - return wrap(b.ravel().astype(result_t)) + + if len(b.shape) == len(a.shape) - 1: + gufunc = _umath_linalg.solve1 else: - return wrap(b.transpose().astype(result_t)) + gufunc = _umath_linalg.solve + + signature = 'DD->D' if isComplexType(t) else 'dd->d' + extobj = get_linalg_error_extobj(_raise_linalgerror_singular) + r = gufunc(a, b, signature=signature, extobj=extobj) + + return wrap(r.astype(result_t)) def tensorinv(a, ind=2): @@ -414,24 +445,29 @@ def inv(a): Parameters ---------- - a : (M, M) array_like + a : (..., M, M) array_like Matrix to be inverted. Returns ------- - ainv : (M, M) ndarray or matrix + ainv : (..., M, M) ndarray or matrix (Multiplicative) inverse of the matrix `a`. Raises ------ LinAlgError - If `a` is singular or not square. + If `a` is not square or inversion fails. + + Notes + ----- + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. Examples -------- - >>> from numpy import linalg as LA + >>> from numpy.linalg import inv >>> a = np.array([[1., 2.], [3., 4.]]) - >>> ainv = LA.inv(a) + >>> ainv = inv(a) >>> np.allclose(np.dot(a, ainv), np.eye(2)) True >>> np.allclose(np.dot(ainv, a), np.eye(2)) @@ -439,14 +475,30 @@ def inv(a): If a is a matrix object, then the return value is a matrix as well: - >>> ainv = LA.inv(np.matrix(a)) + >>> 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]], + [[-5. , 2. ], + [ 3. , -1. ]]]) + """ a, wrap = _makearray(a) - return wrap(solve(a, identity(a.shape[0], dtype=a.dtype))) + _assertNonEmpty(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) + t, result_t = _commonType(a) + signature = 'D->D' if isComplexType(t) else 'd->d' + extobj = get_linalg_error_extobj(_raise_linalgerror_singular) + ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj) + return wrap(ainv.astype(result_t)) # Cholesky decomposition @@ -463,15 +515,15 @@ def cholesky(a): Parameters ---------- - a : (M, M) array_like + a : (..., M, M) array_like Hermitian (symmetric if all elements are real), positive-definite input matrix. Returns ------- - L : {(M, M) ndarray, (M, M) matrix} - Lower-triangular Cholesky factor of `a`. Returns a matrix object - if `a` is a matrix object. + L : (..., M, M) array_like + Upper or lower-triangular Cholesky factor of `a`. Returns a + matrix object if `a` is a matrix object. Raises ------ @@ -481,6 +533,9 @@ def cholesky(a): 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} @@ -518,26 +573,14 @@ def cholesky(a): [ 0.+2.j, 1.+0.j]]) """ + extobj = get_linalg_error_extobj(_raise_linalgerror_nonposdef) + gufunc = _umath_linalg.cholesky_lo a, wrap = _makearray(a) - _assertRank2(a) - _assertSquareness(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) t, result_t = _commonType(a) - a = _fastCopyAndTranspose(t, a) - a = _to_native_byte_order(a) - m = a.shape[0] - n = a.shape[1] - if isComplexType(t): - lapack_routine = lapack_lite.zpotrf - else: - lapack_routine = lapack_lite.dpotrf - results = lapack_routine(_L, n, a, m, 0) - if results['info'] > 0: - raise LinAlgError('Matrix is not positive definite - ' - 'Cholesky decomposition cannot be computed') - s = triu(a, k=0).transpose() - if (s.dtype != result_t): - s = s.astype(result_t) - return wrap(s) + signature = 'D->D' if isComplexType(t) else 'd->d' + return wrap(gufunc(a, signature=signature, extobj=extobj).astype(result_t)) # QR decompostion @@ -763,12 +806,12 @@ def eigvals(a): Parameters ---------- - a : (M, M) array_like + a : (..., M, M) array_like A complex- or real-valued matrix whose eigenvalues will be computed. Returns ------- - w : (M,) ndarray + w : (..., M,) ndarray The eigenvalues, each repeated according to its multiplicity. They are not necessarily ordered, nor are they necessarily real for real matrices. @@ -786,9 +829,11 @@ def eigvals(a): Notes ----- - This is a simple interface to the LAPACK routines dgeev and zgeev - that sets those routines' flags to return only the eigenvalues of - general real and complex arrays, respectively. + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + + This is implemented using the _geev LAPACK routines which compute + the eigenvalues and eigenvectors of general square arrays. Examples -------- @@ -817,49 +862,25 @@ def eigvals(a): """ a, wrap = _makearray(a) - _assertRank2(a) - _assertSquareness(a) + _assertNonEmpty(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) _assertFinite(a) t, result_t = _commonType(a) - real_t = _linalgRealType(t) - a = _fastCopyAndTranspose(t, a) - a = _to_native_byte_order(a) - n = a.shape[0] - dummy = zeros((1,), t) - if isComplexType(t): - lapack_routine = lapack_lite.zgeev - w = zeros((n,), t) - rwork = zeros((n,), real_t) - lwork = 1 - work = zeros((lwork,), t) - results = lapack_routine(_N, _N, n, a, n, w, - dummy, 1, dummy, 1, work, -1, rwork, 0) - lwork = int(abs(work[0])) - work = zeros((lwork,), t) - results = lapack_routine(_N, _N, n, a, n, w, - dummy, 1, dummy, 1, work, lwork, rwork, 0) - else: - lapack_routine = lapack_lite.dgeev - wr = zeros((n,), t) - wi = zeros((n,), t) - lwork = 1 - work = zeros((lwork,), t) - results = lapack_routine(_N, _N, n, a, n, wr, wi, - dummy, 1, dummy, 1, work, -1, 0) - lwork = int(work[0]) - work = zeros((lwork,), t) - results = lapack_routine(_N, _N, n, a, n, wr, wi, - dummy, 1, dummy, 1, work, lwork, 0) - if all(wi == 0.): - w = wr + + extobj = get_linalg_error_extobj( + _raise_linalgerror_eigenvalues_nonconvergence) + signature = 'D->D' if isComplexType(t) else 'd->D' + w = _umath_linalg.eigvals(a, signature=signature, extobj=extobj) + + if not isComplexType(t): + if all(w.imag == 0): + w = w.real result_t = _realType(result_t) else: - w = wr+1j*wi result_t = _complexType(result_t) - if results['info'] > 0: - raise LinAlgError('Eigenvalues did not converge') - return w.astype(result_t) + return w.astype(result_t) def eigvalsh(a, UPLO='L'): """ @@ -869,16 +890,16 @@ def eigvalsh(a, UPLO='L'): Parameters ---------- - a : (M, M) array_like + a : (..., M, M) array_like A complex- or real-valued matrix whose eigenvalues are to be computed. UPLO : {'L', 'U'}, optional - Specifies whether the calculation is done with the lower triangular - part of `a` ('L', default) or the upper triangular part ('U'). + Same as `lower`, wth 'L' for lower and 'U' for upper triangular. + Deprecated. Returns ------- - w : (M,) ndarray + w : (..., M,) ndarray The eigenvalues, not necessarily ordered, each repeated according to its multiplicity. @@ -896,9 +917,10 @@ def eigvalsh(a, UPLO='L'): Notes ----- - This is a simple interface to the LAPACK routines dsyevd and zheevd - that sets those routines' flags to return only the eigenvalues of - real symmetric and complex Hermitian arrays, respectively. + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + + The eigenvalues are computed using LAPACK routines _ssyevd, _heevd Examples -------- @@ -908,45 +930,21 @@ def eigvalsh(a, UPLO='L'): array([ 0.17157288+0.j, 5.82842712+0.j]) """ - UPLO = asbytes(UPLO) + + extobj = get_linalg_error_extobj( + _raise_linalgerror_eigenvalues_nonconvergence) + if UPLO == 'L': + gufunc = _umath_linalg.eigvalsh_lo + else: + gufunc = _umath_linalg.eigvalsh_up + a, wrap = _makearray(a) - _assertRank2(a) - _assertSquareness(a) + _assertNonEmpty(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) t, result_t = _commonType(a) - real_t = _linalgRealType(t) - a = _fastCopyAndTranspose(t, a) - a = _to_native_byte_order(a) - n = a.shape[0] - liwork = 5*n+3 - iwork = zeros((liwork,), fortran_int) - if isComplexType(t): - lapack_routine = lapack_lite.zheevd - w = zeros((n,), real_t) - lwork = 1 - work = zeros((lwork,), t) - lrwork = 1 - rwork = zeros((lrwork,), real_t) - results = lapack_routine(_N, UPLO, n, a, n, w, work, -1, - rwork, -1, iwork, liwork, 0) - lwork = int(abs(work[0])) - work = zeros((lwork,), t) - lrwork = int(rwork[0]) - rwork = zeros((lrwork,), real_t) - results = lapack_routine(_N, UPLO, n, a, n, w, work, lwork, - rwork, lrwork, iwork, liwork, 0) - else: - lapack_routine = lapack_lite.dsyevd - w = zeros((n,), t) - lwork = 1 - work = zeros((lwork,), t) - results = lapack_routine(_N, UPLO, n, a, n, w, work, -1, - iwork, liwork, 0) - lwork = int(work[0]) - work = zeros((lwork,), t) - results = lapack_routine(_N, UPLO, n, a, n, w, work, lwork, - iwork, liwork, 0) - if results['info'] > 0: - raise LinAlgError('Eigenvalues did not converge') + signature = 'D->d' if isComplexType(t) else 'd->d' + w = gufunc(a, signature=signature, extobj=extobj) return w.astype(result_t) def _convertarray(a): @@ -964,17 +962,20 @@ def eig(a): Parameters ---------- - a : (M, M) array_like - A square array of real or complex elements. + a : (..., M, M) array + Matrices for which the eigenvalues and right eigenvectors will + be computed Returns ------- - w : (M,) ndarray + w : (..., M) array The eigenvalues, each repeated according to its multiplicity. - The eigenvalues are not necessarily ordered, nor are they - necessarily real for real arrays (though for real arrays - complex-valued eigenvalues should occur in conjugate pairs). - v : (M, M) ndarray + The eigenvalues are not necessarily ordered. The resulting + array will be always be of complex type. When `a` is real + the resulting eigenvalues will be real (0 imaginary part) or + occur in conjugate pairs + + v : (..., M, M) array The normalized (unit "length") eigenvectors, such that the column ``v[:,i]`` is the eigenvector corresponding to the eigenvalue ``w[i]``. @@ -993,9 +994,11 @@ def eig(a): Notes ----- - This is a simple interface to the LAPACK routines dgeev and zgeev - which compute the eigenvalues and eigenvectors of, respectively, - general real- and complex-valued square arrays. + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + + This is implemented using the _geev LAPACK routines which compute + 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 @@ -1066,57 +1069,24 @@ def eig(a): """ a, wrap = _makearray(a) - _assertRank2(a) - _assertSquareness(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) _assertFinite(a) - a, t, result_t = _convertarray(a) # convert to double or cdouble type - a = _to_native_byte_order(a) - real_t = _linalgRealType(t) - n = a.shape[0] - dummy = zeros((1,), t) - if isComplexType(t): - # Complex routines take different arguments - lapack_routine = lapack_lite.zgeev - w = zeros((n,), t) - v = zeros((n, n), t) - lwork = 1 - work = zeros((lwork,), t) - rwork = zeros((2*n,), real_t) - results = lapack_routine(_N, _V, n, a, n, w, - dummy, 1, v, n, work, -1, rwork, 0) - lwork = int(abs(work[0])) - work = zeros((lwork,), t) - results = lapack_routine(_N, _V, n, a, n, w, - dummy, 1, v, n, work, lwork, rwork, 0) + t, result_t = _commonType(a) + + extobj = get_linalg_error_extobj( + _raise_linalgerror_eigenvalues_nonconvergence) + signature = 'D->DD' if isComplexType(t) else 'd->DD' + w, vt = _umath_linalg.eig(a, signature=signature, extobj=extobj) + + if not isComplexType(t) and all(w.imag == 0.0): + w = w.real + vt = vt.real + result_t = _realType(result_t) else: - lapack_routine = lapack_lite.dgeev - wr = zeros((n,), t) - wi = zeros((n,), t) - vr = zeros((n, n), t) - lwork = 1 - work = zeros((lwork,), t) - results = lapack_routine(_N, _V, n, a, n, wr, wi, - dummy, 1, vr, n, work, -1, 0) - lwork = int(work[0]) - work = zeros((lwork,), t) - results = lapack_routine(_N, _V, n, a, n, wr, wi, - dummy, 1, vr, n, work, lwork, 0) - if all(wi == 0.0): - w = wr - v = vr - result_t = _realType(result_t) - else: - w = wr+1j*wi - v = array(vr, w.dtype) - ind = flatnonzero(wi != 0.0) # indices of complex e-vals - for i in range(len(ind)//2): - v[ind[2*i]] = vr[ind[2*i]] + 1j*vr[ind[2*i+1]] - v[ind[2*i+1]] = vr[ind[2*i]] - 1j*vr[ind[2*i+1]] - result_t = _complexType(result_t) + result_t = _complexType(result_t) - if results['info'] > 0: - raise LinAlgError('Eigenvalues did not converge') - vt = v.transpose().astype(result_t) + vt = vt.astype(result_t) return w.astype(result_t), wrap(vt) @@ -1130,17 +1100,18 @@ def eigh(a, UPLO='L'): Parameters ---------- - a : (M, M) array_like - A complex Hermitian or real symmetric matrix. + A : (..., M, M) array + Hermitian/Symmetric matrices whose eigenvalues and + eigenvectors are to be computed. UPLO : {'L', 'U'}, optional Specifies whether the calculation is done with the lower triangular part of `a` ('L', default) or the upper triangular part ('U'). Returns ------- - w : (M,) ndarray + w : (..., M) ndarray The eigenvalues, not necessarily ordered. - v : {(M, M) ndarray, (M, M) matrix} + v : {(..., M, M) ndarray, (..., M, M) matrix} The column ``v[:, i]`` is the normalized eigenvector corresponding to the eigenvalue ``w[i]``. Will return a matrix object if `a` is a matrix object. @@ -1158,9 +1129,11 @@ def eigh(a, UPLO='L'): Notes ----- - This is a simple interface to the LAPACK routines dsyevd and zheevd, - which compute the eigenvalues and eigenvectors of real symmetric and - complex Hermitian arrays, respectively. + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + + The eigenvalues/eigenvectors are computed using LAPACK routines _ssyevd, + _heevd The eigenvalues of real symmetric or complex Hermitian matrices are always real. [1]_ The array `v` of (column) eigenvectors is unitary @@ -1202,46 +1175,24 @@ def eigh(a, UPLO='L'): """ UPLO = asbytes(UPLO) + a, wrap = _makearray(a) - _assertRank2(a) - _assertSquareness(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) t, result_t = _commonType(a) - real_t = _linalgRealType(t) - a = _fastCopyAndTranspose(t, a) - a = _to_native_byte_order(a) - n = a.shape[0] - liwork = 5*n+3 - iwork = zeros((liwork,), fortran_int) - if isComplexType(t): - lapack_routine = lapack_lite.zheevd - w = zeros((n,), real_t) - lwork = 1 - work = zeros((lwork,), t) - lrwork = 1 - rwork = zeros((lrwork,), real_t) - results = lapack_routine(_V, UPLO, n, a, n, w, work, -1, - rwork, -1, iwork, liwork, 0) - lwork = int(abs(work[0])) - work = zeros((lwork,), t) - lrwork = int(rwork[0]) - rwork = zeros((lrwork,), real_t) - results = lapack_routine(_V, UPLO, n, a, n, w, work, lwork, - rwork, lrwork, iwork, liwork, 0) + + extobj = get_linalg_error_extobj( + _raise_linalgerror_eigenvalues_nonconvergence) + if 'L' == UPLO: + gufunc = _umath_linalg.eigh_lo else: - lapack_routine = lapack_lite.dsyevd - w = zeros((n,), t) - lwork = 1 - work = zeros((lwork,), t) - results = lapack_routine(_V, UPLO, n, a, n, w, work, -1, - iwork, liwork, 0) - lwork = int(work[0]) - work = zeros((lwork,), t) - results = lapack_routine(_V, UPLO, n, a, n, w, work, lwork, - iwork, liwork, 0) - if results['info'] > 0: - raise LinAlgError('Eigenvalues did not converge') - at = a.transpose().astype(result_t) - return w.astype(_realType(result_t)), wrap(at) + gufunc = _umath_linalg.eigh_up + + signature = 'D->dD' if isComplexType(t) else 'd->dd' + w, vt = gufunc(a, signature=signature, extobj=extobj) + w = w.astype(_realType(result_t)) + vt = vt.astype(result_t) + return w, wrap(vt) # Singular value decomposition @@ -1255,7 +1206,7 @@ def svd(a, full_matrices=1, compute_uv=1): Parameters ---------- - a : array_like + a : (..., M, N) array_like A real or complex matrix of shape (`M`, `N`) . full_matrices : bool, optional If True (default), `u` and `v` have the shapes (`M`, `M`) and @@ -1267,15 +1218,14 @@ def svd(a, full_matrices=1, compute_uv=1): Returns ------- - u : ndarray - Unitary matrix. The shape of `u` is (`M`, `M`) or (`M`, `K`) - depending on value of ``full_matrices``. - s : ndarray - The singular values, sorted so that ``s[i] >= s[i+1]``. `s` is - a 1-d array of length min(`M`, `N`). - v : ndarray - Unitary matrix of shape (`N`, `N`) or (`K`, `N`), depending on - ``full_matrices``. + u : { (..., M, M), (..., M, K) } array + Unitary matrices. The actual shape depends on the value of + ``full_matrices``. Only returned when ``compute_uv`` is True. + s : (..., K) array + The singular values for every matrix, sorted in descending order. + v : { (..., N, N), (..., K, N) } array + Unitary matrices. The actual shape depends on the value of + ``full_matrices``. Only returned when ``compute_uv`` is True. Raises ------ @@ -1284,6 +1234,11 @@ def svd(a, full_matrices=1, compute_uv=1): Notes ----- + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + + The decomposition is performed using LAPACK routine _gesdd + The SVD is commonly written as ``a = U S V.H``. The `v` returned by this function is ``V.H`` and ``u = U``. @@ -1323,63 +1278,41 @@ def svd(a, full_matrices=1, compute_uv=1): """ a, wrap = _makearray(a) - _assertRank2(a) _assertNonEmpty(a) - m, n = a.shape + _assertRankAtLeast2(a) t, result_t = _commonType(a) - real_t = _linalgRealType(t) - a = _fastCopyAndTranspose(t, a) - a = _to_native_byte_order(a) - s = zeros((min(n, m),), real_t) + + extobj = get_linalg_error_extobj(_raise_linalgerror_svd_nonconvergence) + + m = a.shape[-2] + n = a.shape[-1] if compute_uv: if full_matrices: - nu = m - nvt = n - option = _A + if m < n: + gufunc = _umath_linalg.svd_m_f + else: + gufunc = _umath_linalg.svd_n_f else: - nu = min(n, m) - nvt = min(n, m) - option = _S - u = zeros((nu, m), t) - vt = zeros((n, nvt), t) - else: - option = _N - nu = 1 - nvt = 1 - u = empty((1, 1), t) - vt = empty((1, 1), t) + if m < n: + gufunc = _umath_linalg.svd_m_s + else: + gufunc = _umath_linalg.svd_n_s - iwork = zeros((8*min(m, n),), fortran_int) - if isComplexType(t): - lapack_routine = lapack_lite.zgesdd - lrwork = min(m,n)*max(5*min(m,n)+7, 2*max(m,n)+2*min(m,n)+1) - rwork = zeros((lrwork,), real_t) - lwork = 1 - work = zeros((lwork,), t) - results = lapack_routine(option, m, n, a, m, s, u, m, vt, nvt, - work, -1, rwork, iwork, 0) - lwork = int(abs(work[0])) - work = zeros((lwork,), t) - results = lapack_routine(option, m, n, a, m, s, u, m, vt, nvt, - work, lwork, rwork, iwork, 0) - else: - lapack_routine = lapack_lite.dgesdd - lwork = 1 - work = zeros((lwork,), t) - results = lapack_routine(option, m, n, a, m, s, u, m, vt, nvt, - work, -1, iwork, 0) - lwork = int(work[0]) - work = zeros((lwork,), t) - results = lapack_routine(option, m, n, a, m, s, u, m, vt, nvt, - work, lwork, iwork, 0) - if results['info'] > 0: - raise LinAlgError('SVD did not converge') - s = s.astype(_realType(result_t)) - if compute_uv: - u = u.transpose().astype(result_t) - vt = vt.transpose().astype(result_t) + signature = 'D->DdD' if isComplexType(t) else 'd->ddd' + u, s, vt = gufunc(a, signature=signature, extobj=extobj) + u = u.astype(result_t) + s = s.astype(_realType(result_t)) + vt = vt.astype(result_t) return wrap(u), s, wrap(vt) else: + if m < n: + gufunc = _umath_linalg.svd_m + else: + gufunc = _umath_linalg.svd_n + + signature = 'D->d' if isComplexType(t) else 'd->d' + s = gufunc(a, signature=signature, extobj=extobj) + s = s.astype(_realType(result_t)) return s def cond(x, p=None): @@ -1649,16 +1582,16 @@ def slogdet(a): Parameters ---------- - a : array_like + a : (..., M, M) array_like Input array, has to be a square 2-D array. Returns ------- - sign : float or complex + sign : (...) array_like A number representing the sign of the determinant. For a real matrix, this is 1, 0, or -1. For a complex matrix, this is a complex number with absolute value 1 (i.e., it is on the unit circle), or else 0. - logdet : float + logdet : (...) array_like The natural log of the absolute value of the determinant. If the determinant is zero, then `sign` will be 0 and `logdet` will be @@ -1670,6 +1603,9 @@ def slogdet(a): Notes ----- + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + The determinant is computed via LU factorization using the LAPACK routine z/dgetrf. @@ -1686,6 +1622,17 @@ def slogdet(a): >>> sign * np.exp(logdet) -2.0 + Computing log-determinants for a stack of matrices: + + >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ]) + >>> a.shape + (3, 2, 2) + >>> sign, logdet = np.linalg.slogdet(a) + >>> (sign, logdet) + (array([-1., -1., -1.]), array([ 0.69314718, 1.09861229, 2.07944154])) + >>> sign * np.exp(logdet) + array([-2., -3., -8.]) + This routine succeeds where ordinary `det` does not: >>> np.linalg.det(np.eye(500) * 0.1) @@ -1695,30 +1642,14 @@ def slogdet(a): """ a = asarray(a) - _assertRank2(a) - _assertSquareness(a) + _assertNonEmpty(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) t, result_t = _commonType(a) - a = _fastCopyAndTranspose(t, a) - a = _to_native_byte_order(a) - n = a.shape[0] - if isComplexType(t): - lapack_routine = lapack_lite.zgetrf - else: - lapack_routine = lapack_lite.dgetrf - pivots = zeros((n,), fortran_int) - results = lapack_routine(n, n, a, n, pivots, 0) - info = results['info'] - if (info < 0): - raise TypeError("Illegal input to Fortran routine") - elif (info > 0): - return (t(0.0), _realType(t)(-Inf)) - sign = 1. - 2. * (add.reduce(pivots != arange(1, n + 1)) % 2) - d = diagonal(a) - absd = absolute(d) - sign *= multiply.reduce(d / absd) - log(absd, absd) - logdet = add.reduce(absd, axis=-1) - return sign, logdet + real_t = _realType(result_t) + signature = 'D->Dd' if isComplexType(t) else 'd->dd' + sign, logdet = _umath_linalg.slogdet(a, signature=signature) + return sign.astype(result_t), logdet.astype(real_t) def det(a): """ @@ -1726,12 +1657,12 @@ def det(a): Parameters ---------- - a : (M, M) array_like - Input array. + a : (..., M, M) array_like + Input array to compute determinants for. Returns ------- - det : float + det : (...) array_like Determinant of `a`. See Also @@ -1741,6 +1672,9 @@ def det(a): Notes ----- + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + The determinant is computed via LU factorization using the LAPACK routine z/dgetrf. @@ -1752,9 +1686,22 @@ def det(a): >>> np.linalg.det(a) -2.0 + Computing determinants for a stack of matrices: + + >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ]) + >>> a.shape + (2, 2, 2 + >>> np.linalg.det(a) + array([-2., -3., -8.]) + """ - sign, logdet = slogdet(a) - return sign * exp(logdet) + a = asarray(a) + _assertNonEmpty(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) + t, result_t = _commonType(a) + signature = 'D->D' if isComplexType(t) else 'd->d' + return _umath_linalg.det(a, signature=signature).astype(result_t) # Linear Least Squares |