Merge branch 'master' of https://github.com/numpy/numpy into patch-2

1 files changed, 85 insertions, 59 deletions

diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py
index aba656b5e..a7f1c074d 100644
--- a/numpy/linalg/linalg.py
+++ b/numpy/linalg/linalg.py
@@ -250,15 +250,15 @@ def solve(a, b):
 
     Parameters
     ----------
-    a : array_like, shape (M, M)
+    a : (M, M) array_like
         Coefficient matrix.
-    b : array_like, shape (M,) or (M, N)
+    b : {(M,), (M, N)}, array_like
         Ordinate or "dependent variable" values.
 
     Returns
     -------
-    x : ndarray, shape (M,) or (M, N) depending on b
-        Solution to the system a x = b
+    x : {(M,), (M, N)} ndarray
+        Solution to the system a x = b.  Returned shape is identical to `b`.
 
     Raises
     ------
@@ -410,12 +410,12 @@ def inv(a):
 
     Parameters
     ----------
-    a : array_like, shape (M, M)
+    a : (M, M) array_like
         Matrix to be inverted.
 
     Returns
     -------
-    ainv : ndarray or matrix, shape (M, M)
+    ainv : (M, M) ndarray or matrix
         (Multiplicative) inverse of the matrix `a`.
 
     Raises
@@ -459,14 +459,15 @@ def cholesky(a):
 
     Parameters
     ----------
-    a : array_like, shape (M, M)
+    a : (M, M) array_like
         Hermitian (symmetric if all elements are real), positive-definite
         input matrix.
 
     Returns
     -------
-    L : ndarray, or matrix object if `a` is, shape (M, M)
-        Lower-triangular Cholesky factor of a.
+    L : {(M, M) ndarray, (M, M) matrix}
+        Lower-triangular Cholesky factor of `a`.  Returns a matrix object
+        if `a` is a matrix object.
 
     Raises
     ------
@@ -709,12 +710,12 @@ def eigvals(a):
 
     Parameters
     ----------
-    a : array_like, shape (M, M)
+    a : (M, M) array_like
         A complex- or real-valued matrix whose eigenvalues will be computed.
 
     Returns
     -------
-    w : ndarray, shape (M,)
+    w : (M,) ndarray
         The eigenvalues, each repeated according to its multiplicity.
         They are not necessarily ordered, nor are they necessarily
         real for real matrices.
@@ -815,7 +816,7 @@ def eigvalsh(a, UPLO='L'):
 
     Parameters
     ----------
-    a : array_like, shape (M, M)
+    a : (M, M) array_like
         A complex- or real-valued matrix whose eigenvalues are to be
         computed.
     UPLO : {'L', 'U'}, optional
@@ -824,7 +825,7 @@ def eigvalsh(a, UPLO='L'):
 
     Returns
     -------
-    w : ndarray, shape (M,)
+    w : (M,) ndarray
         The eigenvalues, not necessarily ordered, each repeated according to
         its multiplicity.
 
@@ -910,18 +911,17 @@ def eig(a):
 
     Parameters
     ----------
-    a : array_like, shape (M, M)
+    a : (M, M) array_like
         A square array of real or complex elements.
 
     Returns
     -------
-    w : ndarray, shape (M,)
+    w : (M,) ndarray
         The eigenvalues, each repeated according to its multiplicity.
         The eigenvalues are not necessarily ordered, nor are they
         necessarily real for real arrays (though for real arrays
         complex-valued eigenvalues should occur in conjugate pairs).
-
-    v : ndarray, shape (M, M)
+    v : (M, M) ndarray
         The normalized (unit "length") eigenvectors, such that the
         column ``v[:,i]`` is the eigenvector corresponding to the
         eigenvalue ``w[i]``.
@@ -1077,7 +1077,7 @@ def eigh(a, UPLO='L'):
 
     Parameters
     ----------
-    a : array_like, shape (M, M)
+    a : (M, M) array_like
         A complex Hermitian or real symmetric matrix.
     UPLO : {'L', 'U'}, optional
         Specifies whether the calculation is done with the lower triangular
@@ -1085,11 +1085,12 @@ def eigh(a, UPLO='L'):
 
     Returns
     -------
-    w : ndarray, shape (M,)
+    w : (M,) ndarray
         The eigenvalues, not necessarily ordered.
-    v : ndarray, or matrix object if `a` is, shape (M, M)
+    v : {(M, M) ndarray, (M, M) matrix}
         The column ``v[:, i]`` is the normalized eigenvector corresponding
-        to the eigenvalue ``w[i]``.
+        to the eigenvalue ``w[i]``.  Will return a matrix object if `a` is
+        a matrix object.
 
     Raises
     ------
@@ -1338,7 +1339,7 @@ def cond(x, p=None):
 
     Parameters
     ----------
-    x : array_like, shape (M, N)
+    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:
@@ -1419,38 +1420,64 @@ 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`.
+    Rank of the array is the number of SVD singular values of the array that are
+    greater than `tol`.
 
     Parameters
     ----------
-    M : array_like
+    M : {(M,), (M, N)} array_like
         array of <=2 dimensions
-    tol : {None, float}
+    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() * eps``.
+       set to ``S.max() * max(M.shape) * eps``.
 
     Notes
     -----
-    Golub and van Loan [1]_ define "numerical rank deficiency" as using
-    tol=eps*S[0] (where S[0] is the maximum singular value and thus the
-    2-norm of the matrix). This is one definition of rank deficiency,
-    and the one we use here.  When floating point roundoff is the main
-    concern, then "numerical rank deficiency" is a reasonable choice. In
-    some cases you may prefer other definitions. 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.
+    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] G. H. Golub and C. F. Van Loan, *Matrix Computations*.
-       Baltimore: Johns Hopkins University Press, 1996.
+    .. [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
     --------
@@ -1464,7 +1491,6 @@ def matrix_rank(M, tol=None):
     1
     >>> matrix_rank(np.zeros((4,)))
     0
-
     """
     M = asarray(M)
     if M.ndim > 2:
@@ -1473,7 +1499,7 @@ def matrix_rank(M, tol=None):
         return int(not all(M==0))
     S = svd(M, compute_uv=False)
     if tol is None:
-        tol = S.max() * finfo(S.dtype).eps
+        tol = S.max() * max(M.shape) * finfo(S.dtype).eps
     return sum(S > tol)
 
 
@@ -1489,7 +1515,7 @@ def pinv(a, rcond=1e-15 ):
 
     Parameters
     ----------
-    a : array_like, shape (M, N)
+    a : (M, N) array_like
       Matrix to be pseudo-inverted.
     rcond : float
       Cutoff for small singular values.
@@ -1499,7 +1525,7 @@ def pinv(a, rcond=1e-15 ):
 
     Returns
     -------
-    B : ndarray, shape (N, M)
+    B : (N, M) ndarray
       The pseudo-inverse of `a`. If `a` is a `matrix` instance, then so
       is `B`.
 
@@ -1647,14 +1673,19 @@ def det(a):
 
     Parameters
     ----------
-    a : array_like, shape (M, M)
+    a : (M, M) array_like
         Input array.
 
     Returns
     -------
-    det : ndarray
+    det : float
         Determinant of `a`.
 
+    See Also
+    --------
+    slogdet : Another way to representing the determinant, more suitable
+      for large matrices where underflow/overflow may occur.
+
     Notes
     -----
     The determinant is computed via LU factorization using the LAPACK
@@ -1668,11 +1699,6 @@ def det(a):
     >>> np.linalg.det(a)
     -2.0
 
-    See Also
-    --------
-    slogdet : Another way to representing the determinant, more suitable
-      for large matrices where underflow/overflow may occur.
-
     """
     sign, logdet = slogdet(a)
     return sign * exp(logdet)
@@ -1693,9 +1719,9 @@ def lstsq(a, b, rcond=-1):
 
     Parameters
     ----------
-    a : array_like, shape (M, N)
+    a : (M, N) array_like
         "Coefficient" matrix.
-    b : array_like, shape (M,) or (M, K)
+    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`.
@@ -1706,18 +1732,18 @@ def lstsq(a, b, rcond=-1):
 
     Returns
     -------
-    x : ndarray, shape (N,) or (N, K)
+    x : {(M,), (M, K)} ndarray
         Least-squares solution.  The shape of `x` depends on the shape of
         `b`.
-    residues : ndarray, shape (), (1,), or (K,)
-        Sums of residues; squared Euclidean 2-norm for each column in
+    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, 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 : ndarray, shape (min(M,N),)
+    s : (min(M, N),) ndarray
         Singular values of `a`.
 
     Raises
@@ -1849,7 +1875,7 @@ def norm(x, ord=None):
 
     Parameters
     ----------
-    x : array_like, shape (M,) or (M, N)
+    x : {(M,), (M, N)} array_like
         Input array.
     ord : {non-zero int, inf, -inf, 'fro'}, optional
         Order of the norm (see table under ``Notes``). inf means numpy's