Mercurial > hg > vamp-build-and-test
diff DEPENDENCIES/mingw32/Python27/Lib/site-packages/numpy/linalg/linalg.py @ 87:2a2c65a20a8b
Add Python libs and headers
| author | Chris Cannam |
|---|---|
| date | Wed, 25 Feb 2015 14:05:22 +0000 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DEPENDENCIES/mingw32/Python27/Lib/site-packages/numpy/linalg/linalg.py Wed Feb 25 14:05:22 2015 +0000 @@ -0,0 +1,2136 @@ +"""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. +""" +from __future__ import division, absolute_import, print_function + + +__all__ = ['matrix_power', 'solve', 'tensorsolve', 'tensorinv', 'inv', + 'cholesky', 'eigvals', 'eigvalsh', 'pinv', 'slogdet', 'det', + 'svd', 'eig', 'eigh', 'lstsq', 'norm', 'qr', 'cond', 'matrix_rank', + 'LinAlgError'] + +import warnings + +from numpy.core import ( + array, asarray, zeros, empty, empty_like, transpose, intc, single, double, + csingle, cdouble, inexact, complexfloating, newaxis, ravel, all, Inf, dot, + add, multiply, sqrt, maximum, fastCopyAndTranspose, sum, isfinite, size, + finfo, errstate, geterrobj, longdouble, rollaxis, amin, amax, product, abs, + broadcast + ) +from numpy.lib import triu, asfarray +from numpy.linalg import lapack_lite, _umath_linalg +from numpy.matrixlib.defmatrix import matrix_power +from numpy.compat import asbytes + +# For Python2/3 compatibility +_N = asbytes('N') +_V = asbytes('V') +_A = asbytes('A') +_S = asbytes('S') +_L = asbytes('L') + +fortran_int = intc + +# Error object +class LinAlgError(Exception): + """ + Generic Python-exception-derived object raised by linalg functions. + + General purpose exception class, derived from Python's exception.Exception + 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 + + """ + 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__) + 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 _linalgRealType(t): + """Cast the type t to either double or cdouble.""" + return double + +_complex_types_map = {single : csingle, + double : cdouble, + csingle : csingle, + cdouble : cdouble} + +def _commonType(*arrays): + # in lite version, use higher precision (always double or cdouble) + result_type = single + is_complex = False + for a in arrays: + if issubclass(a.dtype.type, inexact): + if isComplexType(a.dtype.type): + is_complex = True + rt = _realType(a.dtype.type, default=None) + if rt is None: + # unsupported inexact scalar + raise TypeError("array type %s is unsupported in linalg" % + (a.dtype.name,)) + else: + rt = double + if rt is double: + result_type = double + if is_complex: + t = cdouble + result_type = _complex_types_map[result_type] + else: + t = double + return t, result_type + + +# _fastCopyAndTranpose assumes the input is 2D (as all the calls in here are). + +_fastCT = fastCopyAndTranspose + +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 _fastCopyAndTranspose(type, *arrays): + cast_arrays = () + for a in arrays: + if a.dtype.type is type: + cast_arrays = cast_arrays + (_fastCT(a),) + else: + cast_arrays = cast_arrays + (_fastCT(a.astype(type)),) + if len(cast_arrays) == 1: + return cast_arrays[0] + else: + return cast_arrays + +def _assertRank2(*arrays): + for a in arrays: + if len(a.shape) != 2: + 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()): + raise LinAlgError("Array must not contain infs or NaNs") + +def _assertNoEmpty2d(*arrays): + for a in arrays: + if a.size == 0 and product(a.shape[-2:]) == 0: + raise LinAlgError("Arrays cannot be empty") + + +# Linear equations + +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=len(b.shape))``. + + 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 + -------- + tensordot, tensorinv, einsum + + Examples + -------- + >>> a = np.eye(2*3*4) + >>> a.shape = (2*3, 4, 2, 3, 4) + >>> b = np.random.randn(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(0, 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 + + a = a.reshape(-1, prod) + b = b.ravel() + res = wrap(solve(a, b)) + res.shape = oldshape + return res + +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 identical to `b`. + + Raises + ------ + LinAlgError + If `a` is singular or not square. + + 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. + + 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 ``3 * x0 + x1 = 9`` and ``x0 + 2 * x1 = 8``: + + >>> a = np.array([[3,1], [1,2]]) + >>> b = np.array([9,8]) + >>> x = np.linalg.solve(a, b) + >>> x + array([ 2., 3.]) + + Check that the solution is correct: + + >>> np.allclose(np.dot(a, x), b) + True + + """ + a, _ = _makearray(a) + _assertRankAtLeast2(a) + _assertNdSquareness(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 == a.ndim - 1: + if a.shape[-1] == 0 and b.shape[-1] == 0: + # Legal, but the ufunc cannot handle the 0-sized inner dims + # let the ufunc handle all wrong cases. + a = a.reshape(a.shape[:-1]) + bc = broadcast(a, b) + return wrap(empty(bc.shape, dtype=result_t)) + + gufunc = _umath_linalg.solve1 + else: + if b.size == 0: + if (a.shape[-1] == 0 and b.shape[-2] == 0) or b.shape[-1] == 0: + a = a[:,:1].reshape(a.shape[:-1] + (1,)) + bc = broadcast(a, b) + return wrap(empty(bc.shape, dtype=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): + """ + 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 + -------- + tensordot, tensorsolve + + Examples + -------- + >>> a = np.eye(4*6) + >>> a.shape = (4, 6, 8, 3) + >>> ainv = np.linalg.tensorinv(a, ind=2) + >>> ainv.shape + (8, 3, 4, 6) + >>> b = np.random.randn(4, 6) + >>> np.allclose(np.tensordot(ainv, b), np.linalg.tensorsolve(a, b)) + True + + >>> a = np.eye(4*6) + >>> a.shape = (24, 8, 3) + >>> ainv = np.linalg.tensorinv(a, ind=1) + >>> ainv.shape + (8, 3, 24) + >>> b = np.random.randn(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 inv(a): + """ + Compute the (multiplicative) inverse of a matrix. + + Given a square matrix `a`, return the matrix `ainv` satisfying + ``dot(a, ainv) = dot(ainv, a) = eye(a.shape[0])``. + + Parameters + ---------- + a : (..., M, M) array_like + Matrix to be inverted. + + Returns + ------- + ainv : (..., M, M) ndarray or matrix + (Multiplicative) inverse of the matrix `a`. + + Raises + ------ + LinAlgError + If `a` is not square or inversion fails. + + Notes + ----- + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + + Examples + -------- + >>> from numpy.linalg import inv + >>> a = np.array([[1., 2.], [3., 4.]]) + >>> ainv = inv(a) + >>> np.allclose(np.dot(a, ainv), np.eye(2)) + True + >>> np.allclose(np.dot(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]], + [[-5. , 2. ], + [ 3. , -1. ]]]) + + """ + a, wrap = _makearray(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) + t, result_t = _commonType(a) + + if a.shape[-1] == 0: + # The inner array is 0x0, the ufunc cannot handle this case + return wrap(empty_like(a, dtype=result_t)) + + 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 + +def cholesky(a): + """ + Cholesky decomposition. + + 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. + + Parameters + ---------- + a : (..., M, M) array_like + Hermitian (symmetric if all elements are real), positive-definite + input matrix. + + Returns + ------- + L : (..., M, M) array_like + Upper or lower-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. + + 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 + -------- + >>> 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 + >>> LA.cholesky(np.matrix(A)) + matrix([[ 1.+0.j, 0.+0.j], + [ 0.+2.j, 1.+0.j]]) + + """ + extobj = get_linalg_error_extobj(_raise_linalgerror_nonposdef) + gufunc = _umath_linalg.cholesky_lo + a, wrap = _makearray(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) + t, result_t = _commonType(a) + signature = 'D->D' if isComplexType(t) else 'd->d' + return wrap(gufunc(a, signature=signature, extobj=extobj).astype(result_t)) + +# QR decompostion + +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) + Matrix to be factored. + mode : {'reduced', 'complete', 'r', 'raw', 'full', 'economic'}, optional + If K = min(M, N), then + + 'reduced' : returns q, r with dimensions (M, K), (K, N) (default) + '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,) + 'full' : alias of 'reduced', deprecated + 'economic' : returns h from 'raw', deprecated. + + 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 + returned in 'raw' mode is transposed for calling Fortran. The + 'economic' mode is deprecated. The modes 'full' and 'economic' may + be passed using only the first letter for backwards compatibility, + but all others must be spelled out. See the Notes for more + explanation. + + + Returns + ------- + q : ndarray of float or complex, optional + A matrix with orthonormal columns. When mode = 'complete' the + result is an orthogonal/unitary matrix depending on whether or not + a is real/complex. The determinant may be either +/- 1 in that + case. + r : ndarray of float or complex, optional + The upper-triangular matrix. + (h, tau) : ndarrays of np.double or np.cdouble, optional + The array h contains the Householder reflectors that generate q + along with r. The tau array contains scaling factors for the + reflectors. In the deprecated 'economic' mode only h is returned. + + Raises + ------ + LinAlgError + If factoring fails. + + Notes + ----- + This is an interface to the LAPACK routines dgeqrf, zgeqrf, + dorgqr, and zungqr. + + For more information on the qr factorization, see for example: + http://en.wikipedia.org/wiki/QR_factorization + + Subclasses of `ndarray` are preserved except for the 'raw' mode. So if + `a` is of type `matrix`, all the return values will be matrices too. + + New 'reduced', 'complete', and 'raw' options for mode were added in + Numpy 1.8 and the old option 'full' was made an alias of 'reduced'. In + addition the options 'full' and 'economic' were deprecated. Because + 'full' was the previous default and 'reduced' is the new default, + backward compatibility can be maintained by letting `mode` default. + The 'raw' option was added so that LAPACK routines that can multiply + arrays by q using the Householder reflectors can be used. Note that in + this case the returned arrays are of type np.double or np.cdouble and + the h array is transposed to be FORTRAN compatible. No routines using + the 'raw' return are currently exposed by numpy, but some are available + in lapack_lite and just await the necessary work. + + Examples + -------- + >>> a = np.random.randn(9, 6) + >>> q, r = np.linalg.qr(a) + >>> np.allclose(a, np.dot(q, r)) # a does equal qr + True + >>> r2 = np.linalg.qr(a, mode='r') + >>> r3 = np.linalg.qr(a, mode='economic') + >>> np.allclose(r, r2) # mode='r' returns the same r as mode='full' + True + >>> # But only triu parts are guaranteed equal when mode='economic' + >>> np.allclose(r, np.triu(r3[:6,:6], k=0)) + True + + Example illustrating a common use of `qr`: solving of least squares + problems + + What are the least-squares-best `m` and `y0` in ``y = y0 + mx`` for + the following data: {(0,1), (1,0), (1,2), (2,1)}. (Graph the points + and you'll see that it should be y0 = 0, m = 1.) The answer is provided + by solving the over-determined matrix equation ``Ax = b``, where:: + + A = array([[0, 1], [1, 1], [1, 1], [2, 1]]) + x = array([[y0], [m]]) + b = array([[1], [0], [2], [1]]) + + If A = qr such that q is orthonormal (which is always possible via + Gram-Schmidt), then ``x = inv(r) * (q.T) * b``. (In numpy practice, + however, we simply use `lstsq`.) + + >>> A = np.array([[0, 1], [1, 1], [1, 1], [2, 1]]) + >>> A + array([[0, 1], + [1, 1], + [1, 1], + [2, 1]]) + >>> b = np.array([1, 0, 2, 1]) + >>> q, r = LA.qr(A) + >>> p = np.dot(q.T, b) + >>> np.dot(LA.inv(r), p) + array([ 1.1e-16, 1.0e+00]) + + """ + if mode not in ('reduced', 'complete', 'r', 'raw'): + if mode in ('f', 'full'): + msg = "".join(( + "The 'full' option is deprecated in favor of 'reduced'.\n", + "For backward compatibility let mode default.")) + warnings.warn(msg, DeprecationWarning) + mode = 'reduced' + elif mode in ('e', 'economic'): + msg = "The 'economic' option is deprecated.", + warnings.warn(msg, DeprecationWarning) + mode = 'economic' + else: + raise ValueError("Unrecognized mode '%s'" % mode) + + a, wrap = _makearray(a) + _assertRank2(a) + _assertNoEmpty2d(a) + m, n = a.shape + t, result_t = _commonType(a) + a = _fastCopyAndTranspose(t, a) + a = _to_native_byte_order(a) + mn = min(m, n) + tau = zeros((mn,), t) + if isComplexType(t): + lapack_routine = lapack_lite.zgeqrf + routine_name = 'zgeqrf' + else: + lapack_routine = lapack_lite.dgeqrf + routine_name = 'dgeqrf' + + # calculate optimal size of work data 'work' + lwork = 1 + work = zeros((lwork,), t) + results = lapack_routine(m, n, a, m, tau, work, -1, 0) + if results['info'] != 0: + raise LinAlgError('%s returns %d' % (routine_name, results['info'])) + + # do qr decomposition + lwork = int(abs(work[0])) + work = zeros((lwork,), t) + results = lapack_routine(m, n, a, m, tau, work, lwork, 0) + if results['info'] != 0: + raise LinAlgError('%s returns %d' % (routine_name, results['info'])) + + # handle modes that don't return q + if mode == 'r': + r = _fastCopyAndTranspose(result_t, a[:, :mn]) + return wrap(triu(r)) + + if mode == 'raw': + return a, tau + + if mode == 'economic': + if t != result_t : + a = a.astype(result_t) + return wrap(a.T) + + # generate q from a + if mode == 'complete' and m > n: + mc = m + q = empty((m, m), t) + else: + mc = mn + q = empty((n, m), t) + q[:n] = a + + if isComplexType(t): + lapack_routine = lapack_lite.zungqr + routine_name = 'zungqr' + else: + lapack_routine = lapack_lite.dorgqr + routine_name = 'dorgqr' + + # determine optimal lwork + lwork = 1 + work = zeros((lwork,), t) + results = lapack_routine(m, mc, mn, q, m, tau, work, -1, 0) + if results['info'] != 0: + raise LinAlgError('%s returns %d' % (routine_name, results['info'])) + + # compute q + lwork = int(abs(work[0])) + work = zeros((lwork,), t) + results = lapack_routine(m, mc, mn, q, m, tau, work, lwork, 0) + if results['info'] != 0: + raise LinAlgError('%s returns %d' % (routine_name, results['info'])) + + q = _fastCopyAndTranspose(result_t, q[:mc]) + r = _fastCopyAndTranspose(result_t, a[:, :mc]) + + return wrap(q), wrap(triu(r)) + + +# Eigenvalues + + +def eigvals(a): + """ + Compute the eigenvalues of a general matrix. + + Main difference between `eigvals` and `eig`: the eigenvectors aren't + returned. + + Parameters + ---------- + a : (..., M, M) array_like + A complex- or real-valued matrix whose eigenvalues will be computed. + + Returns + ------- + w : (..., M,) ndarray + The eigenvalues, each repeated according to its multiplicity. + They are not necessarily ordered, nor are they necessarily + real for real matrices. + + Raises + ------ + LinAlgError + If the eigenvalue computation does not converge. + + See Also + -------- + eig : eigenvalues and right eigenvectors of general arrays + eigvalsh : eigenvalues of symmetric or Hermitian arrays. + eigh : eigenvalues and eigenvectors of symmetric/Hermitian arrays. + + Notes + ----- + 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 + -------- + Illustration, using the fact that the eigenvalues of a diagonal matrix + are its diagonal elements, that multiplying a matrix on the left + by an orthogonal matrix, `Q`, and on the right by `Q.T` (the transpose + of `Q`), preserves the eigenvalues of the "middle" matrix. In other words, + if `Q` is orthogonal, then ``Q * A * Q.T`` has the same eigenvalues as + ``A``: + + >>> from numpy import linalg as LA + >>> x = np.random.random() + >>> Q = np.array([[np.cos(x), -np.sin(x)], [np.sin(x), np.cos(x)]]) + >>> LA.norm(Q[0, :]), LA.norm(Q[1, :]), np.dot(Q[0, :],Q[1, :]) + (1.0, 1.0, 0.0) + + Now multiply a diagonal matrix by Q on one side and by Q.T on the other: + + >>> D = np.diag((-1,1)) + >>> LA.eigvals(D) + array([-1., 1.]) + >>> A = np.dot(Q, D) + >>> A = np.dot(A, Q.T) + >>> LA.eigvals(A) + array([ 1., -1.]) + + """ + a, wrap = _makearray(a) + _assertNoEmpty2d(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) + _assertFinite(a) + t, result_t = _commonType(a) + + 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: + result_t = _complexType(result_t) + + return w.astype(result_t) + +def eigvalsh(a, UPLO='L'): + """ + Compute the eigenvalues of a Hermitian or real symmetric matrix. + + Main difference from eigh: the eigenvectors are not computed. + + Parameters + ---------- + a : (..., M, M) array_like + A complex- or real-valued matrix whose eigenvalues are to be + computed. + UPLO : {'L', 'U'}, optional + Same as `lower`, with 'L' for lower and 'U' for upper triangular. + Deprecated. + + Returns + ------- + w : (..., M,) ndarray + The eigenvalues, not necessarily ordered, each repeated according to + its multiplicity. + + Raises + ------ + LinAlgError + If the eigenvalue computation does not converge. + + See Also + -------- + eigh : eigenvalues and eigenvectors of symmetric/Hermitian arrays. + eigvals : eigenvalues of general real or complex arrays. + eig : eigenvalues and right eigenvectors of general real or complex + arrays. + + Notes + ----- + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + + The eigenvalues are computed using LAPACK routines _ssyevd, _heevd + + Examples + -------- + >>> from numpy import linalg as LA + >>> a = np.array([[1, -2j], [2j, 5]]) + >>> LA.eigvalsh(a) + array([ 0.17157288+0.j, 5.82842712+0.j]) + + """ + UPLO = UPLO.upper() + if UPLO not in ('L', 'U'): + raise ValueError("UPLO argument must be 'L' or 'U'") + + 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) + _assertNoEmpty2d(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) + t, result_t = _commonType(a) + signature = 'D->d' if isComplexType(t) else 'd->d' + w = gufunc(a, signature=signature, extobj=extobj) + return w.astype(_realType(result_t)) + +def _convertarray(a): + t, result_t = _commonType(a) + a = _fastCT(a.astype(t)) + return a, t, result_t + + +# Eigenvectors + + +def eig(a): + """ + Compute the eigenvalues and right eigenvectors of a square array. + + Parameters + ---------- + a : (..., M, M) array + Matrices for which the eigenvalues and right eigenvectors will + be computed + + Returns + ------- + w : (..., M) array + The eigenvalues, each repeated according to its multiplicity. + 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]``. + + Raises + ------ + LinAlgError + If the eigenvalue computation does not converge. + + See Also + -------- + eigvalsh : eigenvalues of a symmetric or Hermitian (conjugate symmetric) + array. + + eigvals : eigenvalues of a non-symmetric array. + + Notes + ----- + 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 + `v` satisfy the equations ``dot(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`. + + 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* + eigenvector of `a`, and, in general, the left and right eigenvectors + of a matrix are not necessarily the (perhaps conjugate) transposes + of each other. + + References + ---------- + G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, FL, + Academic Press, Inc., 1980, Various pp. + + Examples + -------- + >>> from numpy import linalg as LA + + (Almost) trivial example with real e-values and e-vectors. + + >>> w, v = LA.eig(np.diag((1, 2, 3))) + >>> w; v + array([ 1., 2., 3.]) + array([[ 1., 0., 0.], + [ 0., 1., 0.], + [ 0., 0., 1.]]) + + Real matrix possessing complex e-values and e-vectors; note that the + e-values are complex conjugates of each other. + + >>> w, v = LA.eig(np.array([[1, -1], [1, 1]])) + >>> w; v + array([ 1. + 1.j, 1. - 1.j]) + array([[ 0.70710678+0.j , 0.70710678+0.j ], + [ 0.00000000-0.70710678j, 0.00000000+0.70710678j]]) + + Complex-valued matrix with real e-values (but complex-valued e-vectors); + note that a.conj().T = a, i.e., a is Hermitian. + + >>> a = np.array([[1, 1j], [-1j, 1]]) + >>> w, v = LA.eig(a) + >>> w; v + array([ 2.00000000e+00+0.j, 5.98651912e-36+0.j]) # i.e., {2, 0} + array([[ 0.00000000+0.70710678j, 0.70710678+0.j ], + [ 0.70710678+0.j , 0.00000000+0.70710678j]]) + + Be careful about round-off error! + + >>> a = np.array([[1 + 1e-9, 0], [0, 1 - 1e-9]]) + >>> # Theor. e-values are 1 +/- 1e-9 + >>> w, v = LA.eig(a) + >>> w; v + array([ 1., 1.]) + array([[ 1., 0.], + [ 0., 1.]]) + + """ + a, wrap = _makearray(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) + _assertFinite(a) + 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: + result_t = _complexType(result_t) + + vt = vt.astype(result_t) + return w.astype(result_t), wrap(vt) + + +def eigh(a, UPLO='L'): + """ + Return the eigenvalues and eigenvectors of a Hermitian or symmetric matrix. + + Returns two objects, a 1-D array containing the eigenvalues of `a`, and + a 2-D square array or matrix (depending on the input type) of the + corresponding eigenvectors (in columns). + + Parameters + ---------- + 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 + The eigenvalues, not necessarily ordered. + 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. + + Raises + ------ + LinAlgError + If the eigenvalue computation does not converge. + + See Also + -------- + eigvalsh : eigenvalues of symmetric or Hermitian arrays. + eig : eigenvalues and right eigenvectors for non-symmetric arrays. + eigvals : eigenvalues of non-symmetric arrays. + + Notes + ----- + 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 + and `a`, `w`, and `v` satisfy the equations + ``dot(a, v[:, i]) = w[i] * v[:, i]``. + + References + ---------- + .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, + FL, Academic Press, Inc., 1980, pg. 222. + + Examples + -------- + >>> from numpy import linalg as LA + >>> a = np.array([[1, -2j], [2j, 5]]) + >>> a + array([[ 1.+0.j, 0.-2.j], + [ 0.+2.j, 5.+0.j]]) + >>> w, v = LA.eigh(a) + >>> w; v + array([ 0.17157288, 5.82842712]) + array([[-0.92387953+0.j , -0.38268343+0.j ], + [ 0.00000000+0.38268343j, 0.00000000-0.92387953j]]) + + >>> np.dot(a, v[:, 0]) - w[0] * v[:, 0] # verify 1st e-val/vec pair + array([2.77555756e-17 + 0.j, 0. + 1.38777878e-16j]) + >>> np.dot(a, v[:, 1]) - w[1] * v[:, 1] # verify 2nd e-val/vec pair + array([ 0.+0.j, 0.+0.j]) + + >>> A = np.matrix(a) # what happens if input is a matrix object + >>> A + matrix([[ 1.+0.j, 0.-2.j], + [ 0.+2.j, 5.+0.j]]) + >>> w, v = LA.eigh(A) + >>> w; v + array([ 0.17157288, 5.82842712]) + matrix([[-0.92387953+0.j , -0.38268343+0.j ], + [ 0.00000000+0.38268343j, 0.00000000-0.92387953j]]) + + """ + UPLO = UPLO.upper() + if UPLO not in ('L', 'U'): + raise ValueError("UPLO argument must be 'L' or 'U'") + + a, wrap = _makearray(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) + t, result_t = _commonType(a) + + extobj = get_linalg_error_extobj( + _raise_linalgerror_eigenvalues_nonconvergence) + if UPLO == 'L': + gufunc = _umath_linalg.eigh_lo + else: + 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 + +def svd(a, full_matrices=1, compute_uv=1): + """ + Singular Value Decomposition. + + Factors the matrix `a` as ``u * np.diag(s) * v``, where `u` and `v` + are unitary and `s` is a 1-d array of `a`'s singular values. + + Parameters + ---------- + 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 + (`N`, `N`), respectively. Otherwise, the shapes are (`M`, `K`) + and (`K`, `N`), respectively, where `K` = min(`M`, `N`). + compute_uv : bool, optional + Whether or not to compute `u` and `v` in addition to `s`. True + by default. + + Returns + ------- + 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 + ------ + LinAlgError + If SVD computation does not converge. + + 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``. + + If ``U`` is a unitary matrix, it means that it + satisfies ``U.H = inv(U)``. + + The rows of `v` are the eigenvectors of ``a.H a``. The columns + of `u` are the eigenvectors of ``a a.H``. For row ``i`` in + `v` and column ``i`` in `u`, the corresponding eigenvalue is + ``s[i]**2``. + + If `a` is a `matrix` object (as opposed to an `ndarray`), then so + are all the return values. + + Examples + -------- + >>> a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6) + + Reconstruction based on full SVD: + + >>> U, s, V = np.linalg.svd(a, full_matrices=True) + >>> U.shape, V.shape, s.shape + ((9, 9), (6, 6), (6,)) + >>> S = np.zeros((9, 6), dtype=complex) + >>> S[:6, :6] = np.diag(s) + >>> np.allclose(a, np.dot(U, np.dot(S, V))) + True + + Reconstruction based on reduced SVD: + + >>> U, s, V = np.linalg.svd(a, full_matrices=False) + >>> U.shape, V.shape, s.shape + ((9, 6), (6, 6), (6,)) + >>> S = np.diag(s) + >>> np.allclose(a, np.dot(U, np.dot(S, V))) + True + + """ + a, wrap = _makearray(a) + _assertNoEmpty2d(a) + _assertRankAtLeast2(a) + t, result_t = _commonType(a) + + extobj = get_linalg_error_extobj(_raise_linalgerror_svd_nonconvergence) + + m = a.shape[-2] + n = a.shape[-1] + if compute_uv: + if full_matrices: + if m < n: + gufunc = _umath_linalg.svd_m_f + else: + gufunc = _umath_linalg.svd_n_f + else: + if m < n: + gufunc = _umath_linalg.svd_m_s + else: + gufunc = _umath_linalg.svd_n_s + + 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): + """ + Compute the condition number of a matrix. + + This function is capable of returning the condition number using + one of seven different norms, depending on the value of `p` (see + Parameters below). + + Parameters + ---------- + x : (M, N) array_like + The matrix whose condition number is sought. + p : {None, 1, -1, 2, -2, inf, -inf, 'fro'}, optional + Order of the norm: + + ===== ============================ + p norm for matrices + ===== ============================ + None 2-norm, computed directly using the ``SVD`` + 'fro' Frobenius norm + inf max(sum(abs(x), axis=1)) + -inf min(sum(abs(x), axis=1)) + 1 max(sum(abs(x), axis=0)) + -1 min(sum(abs(x), axis=0)) + 2 2-norm (largest sing. value) + -2 smallest singular value + ===== ============================ + + inf means the numpy.inf object, and the Frobenius norm is + the root-of-sum-of-squares norm. + + Returns + ------- + c : {float, inf} + The condition number of the matrix. May be infinite. + + See Also + -------- + numpy.linalg.norm + + Notes + ----- + The condition number of `x` is defined as the norm of `x` times the + norm of the inverse of `x` [1]_; the norm can be the usual L2-norm + (root-of-sum-of-squares) or one of a number of other matrix norms. + + References + ---------- + .. [1] G. Strang, *Linear Algebra and Its Applications*, Orlando, FL, + Academic Press, Inc., 1980, pg. 285. + + Examples + -------- + >>> from numpy import linalg as LA + >>> a = np.array([[1, 0, -1], [0, 1, 0], [1, 0, 1]]) + >>> a + array([[ 1, 0, -1], + [ 0, 1, 0], + [ 1, 0, 1]]) + >>> LA.cond(a) + 1.4142135623730951 + >>> LA.cond(a, 'fro') + 3.1622776601683795 + >>> LA.cond(a, np.inf) + 2.0 + >>> LA.cond(a, -np.inf) + 1.0 + >>> LA.cond(a, 1) + 2.0 + >>> LA.cond(a, -1) + 1.0 + >>> LA.cond(a, 2) + 1.4142135623730951 + >>> LA.cond(a, -2) + 0.70710678118654746 + >>> min(LA.svd(a, compute_uv=0))*min(LA.svd(LA.inv(a), compute_uv=0)) + 0.70710678118654746 + + """ + x = asarray(x) # in case we have a matrix + if p is None: + s = svd(x, compute_uv=False) + return s[0]/s[-1] + else: + return norm(x, p)*norm(inv(x), p) + + +def matrix_rank(M, tol=None): + """ + Return matrix rank of array using SVD method + + Rank of the array is the number of SVD singular values of the array that are + greater than `tol`. + + Parameters + ---------- + M : {(M,), (M, N)} array_like + array of <=2 dimensions + tol : {None, float}, optional + threshold below which SVD values are considered zero. If `tol` is + None, and ``S`` is an array with singular values for `M`, and + ``eps`` is the epsilon value for datatype of ``S``, then `tol` is + set to ``S.max() * max(M.shape) * eps``. + + Notes + ----- + The default threshold to detect rank deficiency is a test on the magnitude + of the singular values of `M`. By default, we identify singular values less + than ``S.max() * max(M.shape) * eps`` as indicating rank deficiency (with + the symbols defined above). This is the algorithm MATLAB uses [1]. It also + appears in *Numerical recipes* in the discussion of SVD solutions for linear + least squares [2]. + + This default threshold is designed to detect rank deficiency accounting for + the numerical errors of the SVD computation. Imagine that there is a column + in `M` that is an exact (in floating point) linear combination of other + columns in `M`. Computing the SVD on `M` will not produce a singular value + exactly equal to 0 in general: any difference of the smallest SVD value from + 0 will be caused by numerical imprecision in the calculation of the SVD. + Our threshold for small SVD values takes this numerical imprecision into + account, and the default threshold will detect such numerical rank + deficiency. The threshold may declare a matrix `M` rank deficient even if + the linear combination of some columns of `M` is not exactly equal to + another column of `M` but only numerically very close to another column of + `M`. + + We chose our default threshold because it is in wide use. Other thresholds + are possible. For example, elsewhere in the 2007 edition of *Numerical + recipes* there is an alternative threshold of ``S.max() * + np.finfo(M.dtype).eps / 2. * np.sqrt(m + n + 1.)``. The authors describe + this threshold as being based on "expected roundoff error" (p 71). + + The thresholds above deal with floating point roundoff error in the + calculation of the SVD. However, you may have more information about the + sources of error in `M` that would make you consider other tolerance values + to detect *effective* rank deficiency. The most useful measure of the + tolerance depends on the operations you intend to use on your matrix. For + example, if your data come from uncertain measurements with uncertainties + greater than floating point epsilon, choosing a tolerance near that + uncertainty may be preferable. The tolerance may be absolute if the + uncertainties are absolute rather than relative. + + References + ---------- + .. [1] MATLAB reference documention, "Rank" + http://www.mathworks.com/help/techdoc/ref/rank.html + .. [2] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, + "Numerical Recipes (3rd edition)", Cambridge University Press, 2007, + page 795. + + Examples + -------- + >>> from numpy.linalg import matrix_rank + >>> matrix_rank(np.eye(4)) # Full rank matrix + 4 + >>> I=np.eye(4); I[-1,-1] = 0. # rank deficient matrix + >>> matrix_rank(I) + 3 + >>> matrix_rank(np.ones((4,))) # 1 dimension - rank 1 unless all 0 + 1 + >>> matrix_rank(np.zeros((4,))) + 0 + """ + M = asarray(M) + if M.ndim > 2: + raise TypeError('array should have 2 or fewer dimensions') + if M.ndim < 2: + return int(not all(M==0)) + S = svd(M, compute_uv=False) + if tol is None: + tol = S.max() * max(M.shape) * finfo(S.dtype).eps + return sum(S > tol) + + +# Generalized inverse + +def pinv(a, rcond=1e-15 ): + """ + Compute the (Moore-Penrose) pseudo-inverse of a matrix. + + Calculate the generalized inverse of a matrix using its + singular-value decomposition (SVD) and including all + *large* singular values. + + Parameters + ---------- + a : (M, N) array_like + Matrix to be pseudo-inverted. + rcond : float + Cutoff for small singular values. + Singular values smaller (in modulus) than + `rcond` * largest_singular_value (again, in modulus) + are set to zero. + + Returns + ------- + B : (N, M) ndarray + The pseudo-inverse of `a`. If `a` is a `matrix` instance, then so + is `B`. + + Raises + ------ + LinAlgError + If the SVD computation does not converge. + + Notes + ----- + The pseudo-inverse of a matrix A, denoted :math:`A^+`, is + defined as: "the matrix that 'solves' [the least-squares problem] + :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then + :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`. + + It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular + value decomposition of A, then + :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are + orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting + of A's so-called singular values, (followed, typically, by + zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix + consisting of the reciprocals of A's singular values + (again, followed by zeros). [1]_ + + References + ---------- + .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, + FL, Academic Press, Inc., 1980, pp. 139-142. + + Examples + -------- + The following example checks that ``a * a+ * a == a`` and + ``a+ * a * a+ == a+``: + + >>> a = np.random.randn(9, 6) + >>> B = np.linalg.pinv(a) + >>> np.allclose(a, np.dot(a, np.dot(B, a))) + True + >>> np.allclose(B, np.dot(B, np.dot(a, B))) + True + + """ + a, wrap = _makearray(a) + _assertNoEmpty2d(a) + a = a.conjugate() + u, s, vt = svd(a, 0) + m = u.shape[0] + n = vt.shape[1] + cutoff = rcond*maximum.reduce(s) + for i in range(min(n, m)): + if s[i] > cutoff: + s[i] = 1./s[i] + else: + s[i] = 0.; + res = dot(transpose(vt), multiply(s[:, newaxis], transpose(u))) + return wrap(res) + +# Determinant + +def slogdet(a): + """ + Compute the sign and (natural) logarithm of the determinant of an array. + + If an array has a very small or very large determinant, than a call to + `det` may overflow or underflow. This routine is more robust against such + issues, because it computes the logarithm of the determinant rather than + the determinant itself. + + Parameters + ---------- + a : (..., M, M) array_like + Input array, has to be a square 2-D array. + + Returns + ------- + 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 : (...) 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 + -Inf. In all cases, the determinant is equal to ``sign * np.exp(logdet)``. + + See Also + -------- + det + + Notes + ----- + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + + The determinant is computed via LU factorization using the LAPACK + routine z/dgetrf. + + .. versionadded:: 1.6.0. + + Examples + -------- + The determinant of a 2-D array ``[[a, b], [c, d]]`` is ``ad - bc``: + + >>> a = np.array([[1, 2], [3, 4]]) + >>> (sign, logdet) = np.linalg.slogdet(a) + >>> (sign, logdet) + (-1, 0.69314718055994529) + >>> 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) + 0.0 + >>> np.linalg.slogdet(np.eye(500) * 0.1) + (1, -1151.2925464970228) + + """ + a = asarray(a) + _assertNoEmpty2d(a) + _assertRankAtLeast2(a) + _assertNdSquareness(a) + t, result_t = _commonType(a) + 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): + """ + Compute the determinant of an array. + + Parameters + ---------- + a : (..., M, M) array_like + Input array to compute determinants for. + + Returns + ------- + det : (...) array_like + Determinant of `a`. + + See Also + -------- + slogdet : Another way to representing the determinant, more suitable + for large matrices where underflow/overflow may occur. + + Notes + ----- + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + + The determinant is computed via LU factorization using the LAPACK + routine z/dgetrf. + + Examples + -------- + The determinant of a 2-D array [[a, b], [c, d]] is ad - bc: + + >>> a = np.array([[1, 2], [3, 4]]) + >>> 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.]) + + """ + a = asarray(a) + _assertNoEmpty2d(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 + +def lstsq(a, b, rcond=-1): + """ + Return the least-squares solution to a linear matrix equation. + + Solves the equation `a x = b` by computing a vector `x` that + minimizes the Euclidean 2-norm `|| b - a x ||^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). If `a` + is square and of full rank, then `x` (but for round-off error) is + the "exact" solution of the equation. + + Parameters + ---------- + a : (M, N) array_like + "Coefficient" matrix. + b : {(M,), (M, K)} array_like + Ordinate or "dependent variable" values. If `b` is two-dimensional, + the least-squares solution is calculated for each of the `K` columns + of `b`. + rcond : float, optional + Cut-off ratio for small singular values of `a`. + Singular values are set to zero if they are smaller than `rcond` + times the largest singular value of `a`. + + Returns + ------- + x : {(N,), (N, K)} ndarray + Least-squares solution. If `b` is two-dimensional, + the solutions are in the `K` columns of `x`. + residuals : {(), (1,), (K,)} ndarray + Sums of residuals; squared Euclidean 2-norm for each column in + ``b - a*x``. + If the rank of `a` is < N or M <= N, this is an empty array. + If `b` is 1-dimensional, this is a (1,) shape array. + Otherwise the shape is (K,). + rank : int + Rank of matrix `a`. + s : (min(M, N),) ndarray + Singular values of `a`. + + Raises + ------ + LinAlgError + If computation does not converge. + + Notes + ----- + If `b` is a matrix, then all array results are returned as matrices. + + Examples + -------- + Fit a line, ``y = mx + c``, through some noisy data-points: + + >>> x = np.array([0, 1, 2, 3]) + >>> y = np.array([-1, 0.2, 0.9, 2.1]) + + By examining the coefficients, we see that the line should have a + gradient of roughly 1 and cut the y-axis at, more or less, -1. + + We can rewrite the line equation as ``y = Ap``, where ``A = [[x 1]]`` + and ``p = [[m], [c]]``. Now use `lstsq` to solve for `p`: + + >>> A = np.vstack([x, np.ones(len(x))]).T + >>> A + array([[ 0., 1.], + [ 1., 1.], + [ 2., 1.], + [ 3., 1.]]) + + >>> m, c = np.linalg.lstsq(A, y)[0] + >>> print m, c + 1.0 -0.95 + + Plot the data along with the fitted line: + + >>> import matplotlib.pyplot as plt + >>> plt.plot(x, y, 'o', label='Original data', markersize=10) + >>> plt.plot(x, m*x + c, 'r', label='Fitted line') + >>> plt.legend() + >>> plt.show() + + """ + import math + a, _ = _makearray(a) + b, wrap = _makearray(b) + is_1d = len(b.shape) == 1 + if is_1d: + b = b[:, newaxis] + _assertRank2(a, b) + m = a.shape[0] + n = a.shape[1] + n_rhs = b.shape[1] + ldb = max(n, m) + if m != b.shape[0]: + raise LinAlgError('Incompatible dimensions') + t, result_t = _commonType(a, b) + result_real_t = _realType(result_t) + real_t = _linalgRealType(t) + bstar = zeros((ldb, n_rhs), t) + bstar[:b.shape[0], :n_rhs] = b.copy() + a, bstar = _fastCopyAndTranspose(t, a, bstar) + a, bstar = _to_native_byte_order(a, bstar) + s = zeros((min(m, n),), real_t) + nlvl = max( 0, int( math.log( float(min(m, n))/2. ) ) + 1 ) + iwork = zeros((3*min(m, n)*nlvl+11*min(m, n),), fortran_int) + if isComplexType(t): + lapack_routine = lapack_lite.zgelsd + lwork = 1 + rwork = zeros((lwork,), real_t) + work = zeros((lwork,), t) + results = lapack_routine(m, n, n_rhs, a, m, bstar, ldb, s, rcond, + 0, work, -1, rwork, iwork, 0) + lwork = int(abs(work[0])) + rwork = zeros((lwork,), real_t) + a_real = zeros((m, n), real_t) + bstar_real = zeros((ldb, n_rhs,), real_t) + results = lapack_lite.dgelsd(m, n, n_rhs, a_real, m, + bstar_real, ldb, s, rcond, + 0, rwork, -1, iwork, 0) + lrwork = int(rwork[0]) + work = zeros((lwork,), t) + rwork = zeros((lrwork,), real_t) + results = lapack_routine(m, n, n_rhs, a, m, bstar, ldb, s, rcond, + 0, work, lwork, rwork, iwork, 0) + else: + lapack_routine = lapack_lite.dgelsd + lwork = 1 + work = zeros((lwork,), t) + results = lapack_routine(m, n, n_rhs, a, m, bstar, ldb, s, rcond, + 0, work, -1, iwork, 0) + lwork = int(work[0]) + work = zeros((lwork,), t) + results = lapack_routine(m, n, n_rhs, a, m, bstar, ldb, s, rcond, + 0, work, lwork, iwork, 0) + if results['info'] > 0: + raise LinAlgError('SVD did not converge in Linear Least Squares') + resids = array([], result_real_t) + if is_1d: + x = array(ravel(bstar)[:n], dtype=result_t, copy=True) + if results['rank'] == n and m > n: + if isComplexType(t): + resids = array([sum(abs(ravel(bstar)[n:])**2)], + dtype=result_real_t) + else: + resids = array([sum((ravel(bstar)[n:])**2)], + dtype=result_real_t) + else: + x = array(transpose(bstar)[:n,:], dtype=result_t, copy=True) + if results['rank'] == n and m > n: + if isComplexType(t): + resids = sum(abs(transpose(bstar)[n:,:])**2, axis=0).astype( + result_real_t) + else: + resids = sum((transpose(bstar)[n:,:])**2, axis=0).astype( + result_real_t) + + st = s[:min(n, m)].copy().astype(result_real_t) + return wrap(x), wrap(resids), results['rank'], st + + +def _multi_svd_norm(x, row_axis, col_axis, op): + """Compute the extreme singular values of the 2-D matrices in `x`. + + This is a private utility function used by numpy.linalg.norm(). + + Parameters + ---------- + x : ndarray + row_axis, col_axis : int + The axes of `x` that hold the 2-D matrices. + op : callable + This should be either numpy.amin or numpy.amax. + + Returns + ------- + result : float or ndarray + If `x` is 2-D, the return values is a float. + Otherwise, it is an array with ``x.ndim - 2`` dimensions. + The return values are either the minimum or maximum of the + singular values of the matrices, depending on whether `op` + is `numpy.amin` or `numpy.amax`. + + """ + if row_axis > col_axis: + row_axis -= 1 + y = rollaxis(rollaxis(x, col_axis, x.ndim), row_axis, -1) + result = op(svd(y, compute_uv=0), axis=-1) + return result + + +def norm(x, ord=None, axis=None): + """ + Matrix or vector norm. + + This function is able to return one of seven different matrix norms, + or one of an infinite number of vector norms (described below), depending + on the value of the ``ord`` parameter. + + Parameters + ---------- + x : array_like + Input array. If `axis` is None, `x` must be 1-D or 2-D. + ord : {non-zero int, inf, -inf, 'fro'}, optional + Order of the norm (see table under ``Notes``). inf means numpy's + `inf` object. + axis : {int, 2-tuple of ints, None}, 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. + + Returns + ------- + n : float or ndarray + Norm of the matrix or vector(s). + + Notes + ----- + For values of ``ord <= 0``, the result is, strictly speaking, not a + mathematical 'norm', but it may still be useful for various numerical + purposes. + + The following norms can be calculated: + + ===== ============================ ========================== + ord norm for matrices norm for vectors + ===== ============================ ========================== + None Frobenius norm 2-norm + 'fro' Frobenius norm -- + inf max(sum(abs(x), axis=1)) max(abs(x)) + -inf min(sum(abs(x), axis=1)) min(abs(x)) + 0 -- sum(x != 0) + 1 max(sum(abs(x), axis=0)) as below + -1 min(sum(abs(x), axis=0)) as below + 2 2-norm (largest sing. value) as below + -2 smallest singular value as below + other -- sum(abs(x)**ord)**(1./ord) + ===== ============================ ========================== + + The Frobenius norm is given by [1]_: + + :math:`||A||_F = [\\sum_{i,j} abs(a_{i,j})^2]^{1/2}` + + References + ---------- + .. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*, + Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15 + + Examples + -------- + >>> from numpy import linalg as LA + >>> a = np.arange(9) - 4 + >>> a + array([-4, -3, -2, -1, 0, 1, 2, 3, 4]) + >>> b = a.reshape((3, 3)) + >>> b + array([[-4, -3, -2], + [-1, 0, 1], + [ 2, 3, 4]]) + + >>> LA.norm(a) + 7.745966692414834 + >>> LA.norm(b) + 7.745966692414834 + >>> LA.norm(b, 'fro') + 7.745966692414834 + >>> LA.norm(a, np.inf) + 4 + >>> LA.norm(b, np.inf) + 9 + >>> LA.norm(a, -np.inf) + 0 + >>> LA.norm(b, -np.inf) + 2 + + >>> LA.norm(a, 1) + 20 + >>> LA.norm(b, 1) + 7 + >>> LA.norm(a, -1) + -4.6566128774142013e-010 + >>> LA.norm(b, -1) + 6 + >>> LA.norm(a, 2) + 7.745966692414834 + >>> LA.norm(b, 2) + 7.3484692283495345 + + >>> LA.norm(a, -2) + nan + >>> LA.norm(b, -2) + 1.8570331885190563e-016 + >>> LA.norm(a, 3) + 5.8480354764257312 + >>> LA.norm(a, -3) + nan + + Using the `axis` argument to compute vector norms: + + >>> c = np.array([[ 1, 2, 3], + ... [-1, 1, 4]]) + >>> LA.norm(c, axis=0) + array([ 1.41421356, 2.23606798, 5. ]) + >>> LA.norm(c, axis=1) + array([ 3.74165739, 4.24264069]) + >>> LA.norm(c, ord=1, axis=1) + array([6, 6]) + + Using the `axis` argument to compute matrix norms: + + >>> m = np.arange(8).reshape(2,2,2) + >>> LA.norm(m, axis=(1,2)) + array([ 3.74165739, 11.22497216]) + >>> LA.norm(m[0, :, :]), LA.norm(m[1, :, :]) + (3.7416573867739413, 11.224972160321824) + + """ + x = asarray(x) + + # Check the default case first and handle it immediately. + if ord is None and axis is None: + x = x.ravel(order='K') + if isComplexType(x.dtype.type): + sqnorm = dot(x.real, x.real) + dot(x.imag, x.imag) + else: + sqnorm = dot(x, x) + return sqrt(sqnorm) + + # Normalize the `axis` argument to a tuple. + nd = x.ndim + if axis is None: + axis = tuple(range(nd)) + elif not isinstance(axis, tuple): + axis = (axis,) + + if len(axis) == 1: + if ord == Inf: + return abs(x).max(axis=axis) + elif ord == -Inf: + return abs(x).min(axis=axis) + elif ord == 0: + # Zero norm + return (x != 0).sum(axis=axis) + elif ord == 1: + # special case for speedup + return add.reduce(abs(x), axis=axis) + elif ord is None or ord == 2: + # special case for speedup + s = (x.conj() * x).real + return sqrt(add.reduce(s, axis=axis)) + else: + try: + ord + 1 + except TypeError: + raise ValueError("Invalid norm order for vectors.") + if x.dtype.type is longdouble: + # Convert to a float type, so integer arrays give + # float results. Don't apply asfarray to longdouble arrays, + # because it will downcast to float64. + absx = abs(x) + else: + absx = x if isComplexType(x.dtype.type) else asfarray(x) + if absx.dtype is x.dtype: + absx = abs(absx) + else: + # if the type changed, we can safely overwrite absx + abs(absx, out=absx) + absx **= ord + return add.reduce(absx, axis=axis) ** (1.0 / ord) + elif len(axis) == 2: + row_axis, col_axis = axis + if not (-nd <= row_axis < nd and -nd <= col_axis < nd): + raise ValueError('Invalid axis %r for an array with shape %r' % + (axis, x.shape)) + if row_axis % nd == col_axis % nd: + raise ValueError('Duplicate axes given.') + if ord == 2: + return _multi_svd_norm(x, row_axis, col_axis, amax) + elif ord == -2: + return _multi_svd_norm(x, row_axis, col_axis, amin) + elif ord == 1: + if col_axis > row_axis: + col_axis -= 1 + return add.reduce(abs(x), axis=row_axis).max(axis=col_axis) + elif ord == Inf: + if row_axis > col_axis: + row_axis -= 1 + return add.reduce(abs(x), axis=col_axis).max(axis=row_axis) + elif ord == -1: + if col_axis > row_axis: + col_axis -= 1 + return add.reduce(abs(x), axis=row_axis).min(axis=col_axis) + elif ord == -Inf: + if row_axis > col_axis: + row_axis -= 1 + return add.reduce(abs(x), axis=col_axis).min(axis=row_axis) + elif ord in [None, 'fro', 'f']: + return sqrt(add.reduce((x.conj() * x).real, axis=axis)) + else: + raise ValueError("Invalid norm order for matrices.") + else: + raise ValueError("Improper number of dimensions to norm.")