DOC: Add warning about ill-conditioning to linalg.inv docstring (nump…

…y#24938)

* DOC: Add warning about ill-conditioning to linalg.inv docstring

* DOC: Add wikipedia link and rule of thumb for cond(a)

numpy/linalg/linalg.py

@@ -473,10 +473,10 @@ def _unary_dispatcher(a):
 @array_function_dispatch(_unary_dispatcher)
 def inv(a):
     """
-    Compute the (multiplicative) inverse of a matrix.
+    Compute the inverse of a matrix.
 
     Given a square matrix `a`, return the matrix `ainv` satisfying
-    ``dot(a, ainv) = dot(ainv, a) = eye(a.shape[0])``.
+    ``a @ ainv = ainv @ a = eye(a.shape[0])``.
 
     Parameters
     ----------
@@ -486,7 +486,7 @@ def inv(a):
     Returns
     -------
     ainv : (..., M, M) ndarray or matrix
-        (Multiplicative) inverse of the matrix `a`.
+        Inverse of the matrix `a`.
 
     Raises
     ------
@@ -496,6 +496,8 @@ def inv(a):
     See Also
     --------
     scipy.linalg.inv : Similar function in SciPy.
+    numpy.linalg.cond : Compute the condition number of a matrix.
+    numpy.linalg.svd : Compute the singular value decomposition of a matrix.
 
     Notes
     -----
@@ -505,14 +507,23 @@ def inv(a):
     Broadcasting rules apply, see the `numpy.linalg` documentation for
     details.
 
+    If `a` is detected to be singular, a `LinAlgError` is raised. If `a` is
+    ill-conditioned, a `LinAlgError` may or may not be raised, and results may
+    be innacurate due to floating-point errors.
+
+    References
+    ----------
+    .. [1] Wikipedia, "Condition number",
+           https://en.wikipedia.org/wiki/Condition_number
+    
     Examples
     --------
     >>> from numpy.linalg import inv
     >>> a = np.array([[1., 2.], [3., 4.]])
     >>> ainv = inv(a)
-    >>> np.allclose(np.dot(a, ainv), np.eye(2))
+    >>> np.allclose(a @ ainv, np.eye(2))
     True
-    >>> np.allclose(np.dot(ainv, a), np.eye(2))
+    >>> np.allclose(ainv @ a, np.eye(2))
     True
 
     If a is a matrix object, then the return value is a matrix as well:
@@ -531,6 +542,40 @@ def inv(a):
            [[-1.25,  0.75],
             [ 0.75, -0.25]]])
 
+    If a matrix is close to singular, the computed inverse may not satisfy
+    ``a @ ainv = ainv @ a = eye(a.shape[0])`` even if a `LinAlgError`
+    is not raised:
+
+    >>> a = np.array([[2,4,6],[2,0,2],[6,8,14]])
+    >>> inv(a)  # No errors raised
+    array([[-1.12589991e+15, -5.62949953e+14,  5.62949953e+14],
+       [-1.12589991e+15, -5.62949953e+14,  5.62949953e+14],
+       [ 1.12589991e+15,  5.62949953e+14, -5.62949953e+14]])
+    >>> a @ inv(a)
+    array([[ 0.   , -0.5  ,  0.   ],
+           [-0.5  ,  0.625,  0.25 ],
+           [ 0.   ,  0.   ,  1.   ]])
+
+    To detect ill-conditioned matrices, you can use `numpy.linalg.cond` to
+    compute its *condition number* [1]_. The larger the condition number, the
+    more ill-conditioned the matrix is. As a rule of thumb, if the condition
+    number ``cond(a) = 10**k``, then you may lose up to ``k`` digits of
+    accuracy on top of what would be lost to the numerical method due to loss
+    of precision from arithmetic methods.
+
+    >>> from numpy.linalg import cond
+    >>> cond(a)
+    np.float64(8.659885634118668e+17)
+
+    It is also possible to detect ill-conditioning by inspecting the matrix's
+    singular values directly. The ratio between the largest and the smallest
+    singular value is the condition number:
+
+    >>> from numpy.linalg import svd
+    >>> sigma = svd(a, compute_uv=False)  # Do not compute singular vectors
+    >>> sigma.max()/sigma.min()
+    8.659885634118668e+17
+
     """
     a, wrap = _makearray(a)
     _assert_stacked_2d(a)