"""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 "", line 1, in 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.")