.. _SVDrank:

Rank from the SVD
=================

.. index:: rank

There are multiple ways to compute the rank of a matrix. An example of
using Gaussian elimination to find the rank of a matrix is shown in
:ref:`rank`. The ``numpy.linalg.matrix_rank`` function uses the SVD. 
The SVD algorithm requires more computation
than some alternatives, but it is the most reliable.

.. index::  rank from the singular values theorem

.. _th-rankSVD:

.. rubric::  rank from the singular values theorem

.. prf:theorem::  Rank from the singular values
   :label: th-rankSVD

   The number of singular values greater than a small tolerance is the
   rank of the matrix. As previously noted, the singular values are
   ordered. The first :math:`r = \text{rank}(\mathbf{A})` singular
   values are positive values, and singular values :math:`r+1` to
   :math:`p = \min\{m, n\}` are zero or near zero.

.. prf:proof::

   Any linearly dependent rows and columns of :math:`\bf{A}` are also
   linearly dependent in :math:`\mathbf{A}^T \mathbf{A}`. The
   eigenvalues are the shift to the diagonal of a matrix needed to make
   it singular. So a singular matrix will have at least one zero-valued
   eigenvalue (:ref:`find-eigvals`). Likewise, if :math:`\bf{A}` is
   singular, :math:`\mathbf{A}^T \mathbf{A}` is singular and will have
   :math:`n-r` zero-valued eigenvalues.
   Recall that the singular values come from the eigenvalues of
   :math:`\mathbf{A}^T \mathbf{A}`.

   :math:`\qed`

We can see why a rank-deficient matrix must have singular values near
zero by reforming :math:`\bf{A}` from its factors using outer
products (:ref:`outerProd`). Dependent rows and columns of
:math:`\bf{A}` do not add new information to the matrix, so the
matrix is complete with the sum of :math:`r` rank-1 outer products.
Only the columns of :math:`\bf{U}` and rows of :math:`\mathbf{V}^T`
that are multiplied by nonzero singular values contribute to the
final values of the matrix.

.. math::

   \mathbf{A} = \sum_{i=1}^{r} \sigma_i\,\bm{u}_i\,\bm{v}_i^T + 
   \sum_{i = r+1}^p (\sigma_i = 0)\,\bm{u}_i\,\bm{v}_i^T

We start with a full rank :math:`3{\times}3` matrix in the following
example and reduce the rank by one. The final singular value is near
zero, so the rank is 2.

::

    In [3]: import numpy as np

    In [4]: from scipy.linalg import svd

    In [5]: A = 3*np.random.rand(3, 3)

    In [6]: A[[2],:] = A[[0],:] - A[[1],:]

    In [7]: np.linalg.matrix_rank(A)
    Out[7]: np.int64(2)

    In [8]: _,S,_ = svd(A); print(S)
    [5.0035 1.2798 0.    ]

    In [9]: print(f"{S[2]:e}")
    2.069106e-16