la_25_Jordan_canonical_form.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "#### Jordan canonical form"
      ],
      "metadata": {
        "id": "s9JRN0lOHziu"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "`Any` matrix $A\\in \\mathbf{R}^{n \\times n}$ can be put in `Jordan canonical form`(JCF) through a similarity transformation\n",
        "\n",
        "$$T^{-1}AT=J=\\begin{bmatrix}J_1 & & \\\\\n",
        " & \\ddots & \\\\ & & J_q\\end{bmatrix}$$\n",
        "\n",
        "where each `Jordan block`\n",
        "\n",
        "$$J_i=\\begin{bmatrix}\\lambda_i & 1 &  &\\\\\n",
        " & \\lambda_i & \\ddots & \\\\ & & \\ddots &1 \\\\ & & & \\lambda_i\\end{bmatrix}\\in C^{n_i \\times n_i}$$\n",
        "\n",
        "and $n=\\sum_{i=1}^q n_i$\n",
        "\n",
        "`Diagonal matrix` is a special case of Jordan canonical form with $q=n$ and $n_i=1$\n",
        "\n",
        "Jordan canonical form can have `multiple` blocks with `same eigenvalue`"
      ],
      "metadata": {
        "id": "BxzeF_Ogp0bh"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Determinant and eigenvalues"
      ],
      "metadata": {
        "id": "meAO9JfHtl4O"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "For JCF, we can write the `characteristic polynomial` as\n",
        "\n",
        "$$\\det (sI-A)=(s-\\lambda_1)^{n_1}\\cdots(s-\\lambda_q)^{n_q}$$\n",
        "\n",
        "Here, $\\lambda_i's$ are not necessarily distinct, since the same eigenvalue can occupy different Jordan blocks\n",
        "\n",
        "If we set $s=0$ and recall $\\det (cA) = c^n \\det A$, we see that\n",
        "\n",
        "`Determinant of a matrix is the product of all its eigenvalues`"
      ],
      "metadata": {
        "id": "f1A0CGa1qFcW"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Nullspace of $(\\lambda I-A)$"
      ],
      "metadata": {
        "id": "I-UFQJ0Jt-0t"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "If we write\n",
        "\n",
        "$$(\\lambda I-A)=T(\\lambda I -J)T^{-1}$$\n",
        "\n",
        "Then, for $\\lambda I -J$, each Jordan block would look like\n",
        "\n",
        "$$\\begin{bmatrix}\n",
        "\\lambda-\\lambda_i & -1 & & \\\\ & \\ddots & \\ddots & \\\\\n",
        "& & \\lambda-\\lambda_i & -1 \\\\ & & & \\lambda-\\lambda_i\n",
        "\\end{bmatrix}$$\n",
        "\n",
        "We can see if $\\lambda=\\lambda_i$, then for each Jordan block with $\\lambda_i$, the `rank of the matrix will drop by one`\n",
        "\n",
        "Therefore, $\\text{dim} N(\\lambda I-A)$ is the `number of Jordan blocks` with eigenvalue $\\lambda$ (since similarity transformation preserves matrix rank)"
      ],
      "metadata": {
        "id": "ZH2SiR_ouFR4"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Generalized eigenvectors"
      ],
      "metadata": {
        "id": "pAithOIE75Pp"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "If we have $$T^{-1}AT=J=\\text{diag}(J_1, \\cdots, J_q)$$\n",
        "\n",
        "and if we express\n",
        "\n",
        "$$T=\\begin{bmatrix}T_1 & T_2 & \\cdots & T_q\\end{bmatrix}$$\n",
        "\n",
        "where $T_i \\in C^{n \\times n_i}$ are columns of $T$ associated with the `ith Jordan block`\n",
        "\n",
        "then, we have\n",
        "\n",
        "$$AT_i=T_iJ_i$$\n",
        "\n",
        "If we write out\n",
        "\n",
        "$$T_i=\\begin{bmatrix}v_{i1} & v_{i2} \\cdots & v_{in_i}\\end{bmatrix}$$\n",
        "\n",
        "then, we have\n",
        "\n",
        "$$Av_{i1}=\\lambda_iv_{i1}$$\n",
        "\n",
        "That is, the first column of each $T_i$ is an eigenvector associated with eigenvalue $\\lambda_i$\n",
        "\n",
        "And for $j=2, \\cdots, n_i$, we have\n",
        "\n",
        "$$Av_{ij}=v_{i(j-1)}+\\lambda_iv_{ij}$$\n",
        "\n",
        "The columns in $T_i$ are sometimes called `generalized eigenvectors`"
      ],
      "metadata": {
        "id": "dhyqZl_077nx"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Cayley-Hamilton theorem"
      ],
      "metadata": {
        "id": "pJ2OwQud9kI7"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "If we have a polynomial $p(s)=a_0+a_1s+\\cdots+a_ks^k$, we `define` polynomial of matrix $A\\in \\mathbf{R}^{n \\times n}$ as\n",
        "\n",
        "$$p(A)=a_0I+a_1A+\\cdots+a_kA^k$$\n",
        "\n",
        "Cayley-Hamilton theorem says that for any $A\\in \\mathbf{R}^{n \\times n}$ and $p(s)=\\det (sI-A)$, we have\n",
        "\n",
        "$$p(A)=0$$"
      ],
      "metadata": {
        "id": "whLsM1QK-5hy"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "For example, if we have $A=\\begin{bmatrix} 1 & 2 \\\\ 3 & 4\\end{bmatrix}$\n",
        "\n",
        "Then, $p(s)=(s-1)(s-4)-6=s^2-5s-2$\n",
        "\n",
        "We plug in $A$ and get\n",
        "\n",
        "$$p(A)=\\begin{bmatrix} 7 & 10 \\\\ 15 & 22\\end{bmatrix}-5\\begin{bmatrix} 1 & 2 \\\\ 3 & 4\\end{bmatrix}-2I=0$$\n",
        "\n",
        "What this really says is that $A^2$ is a `linear combination` of $A$ and $I$"
      ],
      "metadata": {
        "id": "MYGu2h_N_cMK"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Based on C-H theorem, we further have the following\n",
        "\n",
        "For every positive integer $p$, we have\n",
        "\n",
        "$$\\boxed{A^p\\in \\text{span}\\{I, A, A^2, \\cdots, A^{n-1}\\}}$$"
      ],
      "metadata": {
        "id": "ZhKLyopdCPMw"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "To show this, we divide $s^p$ by $p(s)$ to get\n",
        "\n",
        "$$s^p=q(s)p(s)+r(s)$$\n",
        "\n",
        "where $r(s)=a_0+a_1s+\\cdots+a_{n-1}s^{n-1}$\n",
        "\n",
        "* See example below for when $p\\geq n$\n",
        "* When $p< n$, we can just set $q(s)=0$ and $r(s)=s^p$\n",
        "\n",
        "Then, we immediately have what we want\n",
        "\n",
        "$$\\begin{align*}A^p&=q(A)p(A)+r(A)\\\\\n",
        "&=q(A)\\cdot 0 +r(A)\\\\\n",
        "&=r(A)\\\\\n",
        "&=\\boxed{a_0I+a_1A+\\cdots+a_{n-1}A^{n-1}}\n",
        "\\end{align*}$$"
      ],
      "metadata": {
        "id": "t3yfhHiSC9yk"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "##### Example of polynomial division"
      ],
      "metadata": {
        "id": "QmrYjOr93ORb"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Suppose $n=3$ and $p(s)=s^3+c_2s^2+c_1s+c_0$, and we want to divide $s^5$ by $p(s)$\n",
        "\n",
        "* We first divide by leading term of $p(s)$ to get first term of $q(s)=s^2$\n",
        "* We multiply $s^2$ to $p(s)$ to get\n",
        "$$s^2p(s)=s^5+c_2s^4+c_1s^3+c_0s^2$$\n",
        "* Subtract this from $s^5$ to get `remainder`\n",
        "$$s^5-(s^5+c_2s^4+c_1s^3+c_0s^2)=-c_2s^4-c_1s^3-c_0s^2$$\n",
        "* Divide by leading term $s^3$ to update $q(s)=s^2-c_2s$\n",
        "* Multiply -c_2s to $p(s)$\n",
        "$$-c_2s(s^3+c_2s^2+c_1s+c_0)=-c_2s^4-c_2^2s^3-c_1c_2s^2-c_0c_2s$$\n",
        "* Subtract this from previous `remainder`\n",
        "$$-c_2s^4-c_1s^3-c_0s^2-(-c_2s^4-c_2^2s^3-c_1c_2s^2-c_0c_2s)=(c_2^2-c_1)s^3+(c_1c_2-c_0)s^2+c_0c_2s$$\n",
        "* Divide by leading term $s^3$ to update $q(s)=\\boxed{s^2-c_2s+(c_2^2-c_1)}$\n",
        "* Multiply $(c_2^2-c_1)$ to $p(s)$\n",
        "$$(c_2^2-c_1)(s^3+c_2s^2+c_1s+c_0)$$\n",
        "* Subtract this from previous `remainder`\n",
        "$$(c_2^2-c_1)s^3+(c_1c_2-c_0)s^2+c_0c_2s-(c_2^2-c_1)(s^3+c_2s^2+c_1s+c_0)$$\n",
        "* This will be a polynomial of degree `less than 3`, and it will be $r(s)$\n",
        "\n",
        "We can use SymPy to verify that $q(s)p(s)+r(s)=s^5$"
      ],
      "metadata": {
        "id": "QpIQUwxt3UHv"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from sympy import symbols, expand, collect"
      ],
      "metadata": {
        "id": "yppnBwvtBgao"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "s, c0, c1, c2 = symbols('s c0 c1 c2')\n",
        "\n",
        "p_s = s**3 + c2*s**2 + c1*s + c0\n",
        "q_s = s**2 - c2*s + (c2**2 - c1)\n",
        "\n",
        "r_s = (c2**2 - c1)*s**3 + (c1*c2 - c0)*s**2 + c0*c2*s - (c2**2 - c1)*p_s\n",
        "\n",
        "q_p_expanded = expand(q_s * p_s)\n",
        "final_expression = expand(q_p_expanded + r_s)\n",
        "final_collected = collect(final_expression, s)\n",
        "\n",
        "print(final_collected)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "x0vB-R22BIPq",
        "outputId": "579c198c-5a54-4b0e-cc97-7cac6b276a72"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "s**5\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "##### Relation to inverse of matrices"
      ],
      "metadata": {
        "id": "JjdRX32PFCgF"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "$$p(A)=a_nA^n+a_{n-1}A^{n-1}+\\cdots+a_0I=0$$\n",
        "\n",
        "If we separate $I$ and assume $A$ is `invertible`, we get\n",
        "\n",
        "$$\\begin{align*}I&=A\\left(-\\frac{a_1}{a_0}I-\\frac{a_2}{a_0}A\\cdots -\\frac{a_n}{a_0}A^{n-1}\\right)\\\\\n",
        "&=AA^{-1}\n",
        "\\end{align*}$$\n",
        "\n",
        "Therefore, the `inverse` of $A$ is a linear combination of $I, A, \\cdots, A^{n-1}$, and the `coefficients are from the characteristic polynomial`\n",
        "\n",
        "Of course, it requires that $a_0\\neq 0$\n",
        "\n",
        "But since $a_0=p(0)=\\det (0I-A)=\\det (-A)$, therefore, when $A$ is invertible, we know that\n",
        "\n",
        "$$a_0=\\det (-A)=(-1)^n \\det A \\neq 0$$"
      ],
      "metadata": {
        "id": "jhGNMqqqFJ5w"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "##### Intuition behind C-H theorem"
      ],
      "metadata": {
        "id": "Th6VR6gcIuQd"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "First, the case where $A$ has distinct eigenvalues, then, it is diagonalizable $A=T\\Lambda T^{-1}$, we want to show first that $p(A)=Tp(\\Lambda)T^{-1}$\n",
        "\n",
        "$$\\begin{align*}\n",
        "p(A)&=a_0I+a_1A+\\cdots+a_nA^n \\\\\n",
        "&= a_0I+a_1T\\Lambda T^{-1}+a_2T\\Lambda^2 T^{-1}+\\cdots+a_nT\\Lambda^n T^{-1} \\\\\n",
        "&=T(a_0I+a_1\\Lambda +a_2\\Lambda^2 +\\cdots+a_n \\Lambda^n )T^{-1} \\\\\n",
        "&=Tp(\\Lambda)T^{-1}\n",
        "\\end{align*}$$"
      ],
      "metadata": {
        "id": "n5qQ3K3pKY0k"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We write it using eigenvalues to get what we want\n",
        "\n",
        "$$\\begin{align*}p(A)&=Tp(\\Lambda)T^{-1} \\\\\n",
        "& p(s)=(s-\\lambda_1)\\cdots(s-\\lambda_n) \\\\\n",
        "&=T(\\Lambda-\\lambda_1 I)\\cdots(\\Lambda-\\lambda_n I)T^{-1}\\\\\n",
        "&=T\\text{diag}(0, \\lambda_2-\\lambda_1, \\cdots, \\lambda_n-\\lambda_1)\\cdots \\text{diag}(\\lambda_1-\\lambda_n, \\cdots, \\lambda_{n-1}-\\lambda_n, 0) T^{-1}\\\\\n",
        "&=0\n",
        "\\end{align*}$$"
      ],
      "metadata": {
        "id": "F29j-t5ly6pV"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "For more general case, we have\n",
        "\n",
        "$$T^{-1}AT=J$$\n",
        "\n",
        "We just need to show for any $J_i$, $p(J_i)=0$\n",
        "\n",
        "Plug into $p(s)=(s-\\lambda_i)^{n_1}\\cdots(s-\\lambda_q)^{n_q}$, we have\n",
        "\n",
        "$$\\begin{align*}p(J_i)&=(J_i-\\lambda_1 I)^{n_1}\\cdots(J_i-\\lambda_i I)^{n_i}\\cdots(J_i-\\lambda_q I)^{n_q}\\\\\n",
        "&=(J_i-\\lambda_1 I)^{n_1}\\cdots\\begin{bmatrix}0 & 1 & 0 & \\cdots \\\\0 & 0 & 1 & \\cdots \\\\ 0 & 0 & 0 & \\ddots\\end{bmatrix}^{n_i}\\cdots(J_i-\\lambda_q I)^{n_q}\n",
        "\\end{align*}$$\n",
        "\n",
        "and for a $n_i \\times n_i$ matrix $\\begin{bmatrix}0 & 1 & 0 & \\cdots \\\\0 & 0 & 1 & \\cdots \\\\ 0 & 0 & 0 & \\ddots\\end{bmatrix}^{n_i}$, we can see that every multiplication `shifts` the 1's towards top right and after power of $n_i$, it becomes 0. Therefore $p(J_i)=0$"
      ],
      "metadata": {
        "id": "Yc71sVwGMaXS"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "HmogTCItMNRI"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}