DOC: linalg: Include information about scipy.linalg. (#14988)

* DOC: Add links to scipy linalg fuctions and compare numpy.linalg vs scipy.linalg.

2 files changed, 60 insertions, 2 deletions

diff --git a/doc/source/reference/routines.linalg.rst b/doc/source/reference/routines.linalg.rst
index d42e77ad8..86e168b26 100644
--- a/doc/source/reference/routines.linalg.rst
+++ b/doc/source/reference/routines.linalg.rst
@@ -18,6 +18,18 @@ or specify the processor architecture.
 .. _OpenBLAS: https://www.openblas.net/
 .. _threadpoolctl: https://github.com/joblib/threadpoolctl
 
+The SciPy library also contains a `~scipy.linalg` submodule, and there is
+overlap in the functionality provided by the SciPy and NumPy submodules.  SciPy
+contains functions not found in `numpy.linalg`, such as functions related to
+LU decomposition and the Schur decomposition, multiple ways of calculating the
+pseudoinverse, and matrix transcendentals such as the matrix logarithm.  Some
+functions that exist in both have augmented functionality in `scipy.linalg`.
+For example, `scipy.linalg.eig` can take a second matrix argument for solving
+generalized eigenvalue problems.  Some functions in NumPy, however, have more
+flexible broadcasting options.  For example, `numpy.linalg.solve` can handle
+"stacked" arrays, while `scipy.linalg.solve` accepts only a single square
+array as its first argument.
+
 .. currentmodule:: numpy
 
 Matrix and vector products
diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py
index 6249a57c2..3d2d17057 100644
--- a/numpy/linalg/linalg.py
+++ b/numpy/linalg/linalg.py
@@ -345,6 +345,10 @@ def solve(a, b):
     LinAlgError
         If `a` is singular or not square.
 
+    See Also
+    --------
+    scipy.linalg.solve : Similar function in SciPy.
+
     Notes
     -----
 
@@ -502,6 +506,10 @@ def inv(a):
     LinAlgError
         If `a` is not square or inversion fails.
 
+    See Also
+    --------
+    scipy.linalg.inv : Similar function in SciPy.
+
     Notes
     -----
 
@@ -700,6 +708,14 @@ def cholesky(a):
        If the decomposition fails, for example, if `a` is not
        positive-definite.
 
+    See Also
+    --------
+    scipy.linalg.cholesky : Similar function in SciPy.
+    scipy.linalg.cholesky_banded : Cholesky decompose a banded Hermitian
+                                   positive-definite matrix.
+    scipy.linalg.cho_factor : Cholesky decomposition of a matrix, to use in
+                              `scipy.linalg.cho_solve`.
+
     Notes
     -----
 
@@ -812,6 +828,11 @@ def qr(a, mode='reduced'):
     LinAlgError
         If factoring fails.
 
+    See Also
+    --------
+    scipy.linalg.qr : Similar function in SciPy.
+    scipy.linalg.rq : Compute RQ decomposition of a matrix.
+
     Notes
     -----
     This is an interface to the LAPACK routines ``dgeqrf``, ``zgeqrf``,
@@ -1004,6 +1025,7 @@ def eigvals(a):
                (conjugate symmetric) arrays.
     eigh : eigenvalues and eigenvectors of real symmetric or complex
            Hermitian (conjugate symmetric) arrays.
+    scipy.linalg.eigvals : Similar function in SciPy.
 
     Notes
     -----
@@ -1105,6 +1127,7 @@ def eigvalsh(a, UPLO='L'):
     eigvals : eigenvalues of general real or complex arrays.
     eig : eigenvalues and right eigenvectors of general real or complex
           arrays.
+    scipy.linalg.eigvalsh : Similar function in SciPy.
 
     Notes
     -----
@@ -1203,12 +1226,12 @@ def eig(a):
     See Also
     --------
     eigvals : eigenvalues of a non-symmetric array.
-
     eigh : eigenvalues and eigenvectors of a real symmetric or complex
            Hermitian (conjugate symmetric) array.
-
     eigvalsh : eigenvalues of a real symmetric or complex Hermitian
                (conjugate symmetric) array.
+    scipy.linalg.eig : Similar function in SciPy (but also solves the
+                       generalized eigenvalue problem).
 
     Notes
     -----
@@ -1355,6 +1378,8 @@ def eigh(a, UPLO='L'):
                (conjugate symmetric) arrays.
     eig : eigenvalues and right eigenvectors for non-symmetric arrays.
     eigvals : eigenvalues of non-symmetric arrays.
+    scipy.linalg.eigh : Similar function in SciPy (but also solves the
+                        generalized eigenvalue problem).
 
     Notes
     -----
@@ -1506,6 +1531,11 @@ def svd(a, full_matrices=True, compute_uv=True, hermitian=False):
     LinAlgError
         If SVD computation does not converge.
 
+    See Also
+    --------
+    scipy.linalg.svd : Similar function in SciPy.
+    scipy.linalg.svdvals : Compute singular values of a matrix.
+
     Notes
     -----
 
@@ -1917,6 +1947,13 @@ def pinv(a, rcond=1e-15, hermitian=False):
     LinAlgError
         If the SVD computation does not converge.
 
+    See Also
+    --------
+    scipy.linalg.pinv : Similar function in SciPy.
+    scipy.linalg.pinv2 : Similar function in SciPy (SVD-based).
+    scipy.linalg.pinvh : Compute the (Moore-Penrose) pseudo-inverse of a
+                         Hermitian matrix.
+
     Notes
     -----
     The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
@@ -2079,6 +2116,7 @@ def det(a):
     --------
     slogdet : Another way to represent the determinant, more suitable
       for large matrices where underflow/overflow may occur.
+    scipy.linalg.det : Similar function in SciPy.
 
     Notes
     -----
@@ -2179,6 +2217,10 @@ def lstsq(a, b, rcond="warn"):
     LinAlgError
         If computation does not converge.
 
+    See Also
+    --------
+    scipy.linalg.lstsq : Similar function in SciPy.
+
     Notes
     -----
     If `b` is a matrix, then all array results are returned as matrices.
@@ -2353,6 +2395,10 @@ def norm(x, ord=None, axis=None, keepdims=False):
     n : float or ndarray
         Norm of the matrix or vector(s).
 
+    See Also
+    --------
+    scipy.linalg.norm : Similar function in SciPy.
+
     Notes
     -----
     For values of ``ord < 1``, the result is, strictly speaking, not a