.. _diagonalization:

Diagonalization
================

Here, we will establish a factoring of a matrix based on its eigenvalues
and eigenvectors. The factoring is called *diagonalization* or the
*eigendecomposition*.

.. index:: eigen-decomposition
.. index:: diagonalization


The eigenpairs of the :math:`n{\times}n` matrix :math:`\bf{A}` give us
:math:`n` equations.

.. math::

   \begin{array}{rl}
   \mathbf{A}\,\bm{x}_1 &= \lambda_1\,\bm{x}_1 \\
   \mathbf{A}\,\bm{x}_2 &= \lambda_2\,\bm{x}_2 \\
     &\vdots \\
   \mathbf{A}\,\bm{x}_n &= \lambda_n\,\bm{x}_n \\
   \end{array}

We wish to combine the :math:`n` equations into one matrix equation.
Consider the product of :math:`\bf{A}` and the :math:`\bf{X}` matrix
containing the eigenvectors.

.. math::

   \begin{array}{ll}
     \mathbf{A\,X} &= \mathbf{A}
         \begin{bmatrix}
            \vertbar{}  & \vertbar{} & {} & \vertbar{} \\ 
            \bm{x}_1  & \bm{x}_2 & \cdots{} & \bm{x}_n \\ 
            \vertbar{}  & \vertbar{} & {} & \vertbar{}
         \end{bmatrix} \\
       &{} \\
       &= \begin{bmatrix}
             \vertbar{}  & \vertbar{} & {} & \vertbar{} \\ 
             {\lambda_1\,\bm{x}_1}  & {\lambda_2\,\bm{x}_2} & \cdots{} & {\lambda_n\,\bm{x}_n} \\ 
             \vertbar{}  & \vertbar{} & {} & \vertbar{}
          \end{bmatrix} \\[18pt]
       &= \begin{bmatrix}
             \vertbar{}  & \vertbar{} & {} & \vertbar{} \\ 
             \bm{x}_1  & \bm{x}_2 & \cdots{} & \bm{x}_n \\ 
             \vertbar{}  & \vertbar{} & {} & \vertbar{}
          \end{bmatrix}
          \begin{bmatrix}
             \lambda_1  & 0 & \cdots{} & 0 \\ 
             0  & \lambda_2 & \cdots{} & 0 \\ 
             \vdots{}  & \vdots{} & \ddots{} & \vdots{} \\ 
             0  & 0 & \cdots{} & \lambda_n
          \end{bmatrix} \\
       &= \mathbf{X\,\Lambda}
   \end{array}

The matrix :math:`\bf{\Lambda}` is a diagonal eigenvalue matrix. If the
matrix has linearly independent eigenvectors, which will be the case
when all eigenvalues are different (no repeating :math:`\lambda`\ s),
then :math:`\bf{X}` is invertible.

.. math::

   \begin{array}{c}
         \mathbf{A\,X} = \mathbf{X\,\Lambda} \\
         \mathbf{X}^{-1}\mathbf{A\,X} = \mathbf{\Lambda} \\
         {}\\
   \boxed{\mathbf{A} = \mathbf{X}\,\mathbf{\Lambda}\,\mathbf{X}^{-1}}
         \end{array}

.. index:: eig

The ``scipy.linalg`` module has a function called ``eig`` that calculates a
matrix’s eigenvalues and eigenvectors. 
The eigenvectors are the columns of matrix ``X``. 
The eigenvalues are returned as a row vector. 
We often want the eigenvalues as a diagonal matrix, so we use the
``diag`` function from NumPy.

::

    In [26]: print(A)
    [[ 2  2]
     [ 2 -1]]

    In [27]: L, X = eig(A); print(f'L = \n{L}\nX =\n{X}')
    L = 
    [ 3.+0.j -2.+0.j]
    X =
    [[ 0.8944 -0.4472]
     [ 0.4472  0.8944]]

    In [28]: L = np.diag(L); print(L)
    [[ 3.+0.j  0.+0.j]
     [ 0.+0.j -2.+0.j]]

When Does Diagonalization Not Work?
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Notice that the inverse of the eigenvector matrix must exist for the
eigendecomposition to be valid. So each eigenvector *must be
independent* of the other eigenvectors.

One example of a matrix with repeating eigenvectors is a
:math:`2{\times}2` matrix with the same values on the forward diagonal
and a zero on the backward diagonal. The determinant of the
characteristic matrix has repeating roots, so it will also have
repeating eigenvectors. Thus, the eigenvector matrix is singular and may
not be inverted.

.. math::

    \begin{bmatrix}
      {a-\lambda}  & n \\ 
      0  & {a-\lambda}
    \end{bmatrix}
    = (a - \lambda) (a - \lambda)

::

    In [36]: A = np.array([[5, 7], [0, 5]])

    In [37]: L, X = eig(A); print(X)
    [[ 1. -1.]
     [ 0.  0.]]

    In [39]: np.linalg.matrix_rank(X)
    Out[39]: np.int64(1)
   
.. _symdiag:

Diagonalization of a Symmetric Matrix
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Symmetric matrices have a simplified diagonalization because the matrix
of eigenvectors of a symmetric matrix, :math:`\bf{S}`, is orthogonal
(:ref:`th-spect-decomp`). Recall also from :ref:`matTranspose` that from
the spectral theorem, orthogonal matrices have the property
:math:`\mathbf{Q}^{T} = \mathbf{Q}^{-1}`. Thus, the diagonalization of a
symmetric matrix is shown below and is also called the 
*spectral decomposition* of a symmetric matrix.

.. math:: \boxed{\mathbf{S} = \mathbf{Q\,\Lambda\,Q}^T}.

.. note::

   Recall that the columns of orthonormal matrices must be unit vectors
   (length of 1), which the ``eig`` function returns.

.. index:: spectral decomposition