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