Chapter 13 - Diagonalization.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sympy as sy\n",
    "sy.init_printing()  \n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('ggplot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Similarity </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $A = PBP^{-1}$, we say $A$ is _similar_ to $B$, decomposing $A$ into $PBP^{-1}$ is also called a _similarity transformation_.\n",
    "\n",
    "If $n\\times n$ matrices $A$ and $B$ are similar, they have the _same eigenvalues_.\n",
    "\n",
    "The _diagnoalization_, which we will explain below,  is a process of finding similar matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Diagonalizable Matrix</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $A$ be an $n\\times n$ matrix. If there exists an $n\\times n$ invertible matrix $P$ and a diagonal matrix $D$, such that\n",
    "\n",
    "$$\n",
    "A=PDP^{-1}\n",
    "$$\n",
    "\n",
    "then matrix $A$ is called a _diagonalizable matrix_."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And further, the columns of $P$ are linearly independent eigenvectors of $A$, and its corresponding eigenvalues are on the diagonal of $D$. In other words, $A$ is diagonalizable if and only if the dimension of eigenspace basis is $n$</font>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's show why this equation holds."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let\n",
    "\n",
    "$$\n",
    "P = \\left[\\begin{array}{llll}\n",
    "{v}_{1} & {v}_{2} & \\cdots & {v}_{n}\n",
    "\\end{array}\\right]\\\\\n",
    "$$\n",
    "\n",
    "$$\n",
    "D = \\left[\\begin{array}{cccc}\n",
    "\\lambda_{1} & 0 & \\cdots & 0 \\\\\n",
    "0 & \\lambda_{2} & \\cdots & 0 \\\\\n",
    "\\vdots & \\vdots & & \\vdots \\\\\n",
    "0 & 0 & \\cdots & \\lambda_{n}\n",
    "\\end{array}\\right]\n",
    "$$\n",
    "\n",
    "where $v_i, i \\in (1, 2, ...n)$ is an eigenvector of $A$, $\\lambda_i, i \\in (1, 2, ...n)$ is an eigenvalue of $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "AP = A\\left[\\begin{array}{llll}\n",
    "{v}_{1} & {v}_{2} & \\cdots & {v}_{n}\n",
    "\\end{array}\\right]=\\left[\\begin{array}{llll}\n",
    "A {v}_{1} & A {v}_{2} & \\cdots & A {v}_{n}\n",
    "\\end{array}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$P D=P\\left[\\begin{array}{cccc}\n",
    "\\lambda_{1} & 0 & \\cdots & 0 \\\\\n",
    "0 & \\lambda_{2} & \\cdots & 0 \\\\\n",
    "\\vdots & \\vdots & & \\vdots \\\\\n",
    "0 & 0 & \\cdots & \\lambda_{n}\n",
    "\\end{array}\\right]=\\left[\\begin{array}{lllll}\n",
    "\\lambda_{1} {v}_{1} & \\lambda_{2} {v}_{2} & \\cdots & \\lambda_{n} {v}_{n}\n",
    "\\end{array}\\right]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know that $A{v}_i = \\lambda_i{v}_i$, i.e.\n",
    "\n",
    "$$\n",
    "AP = PD\n",
    "$$\n",
    "\n",
    "Since $P$ has all independent eigenvectors, then\n",
    "\n",
    "$$\n",
    "A = PDP^{-1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Strictly speaking, if $A$ is symmetric, i.e. $A=A^T$, what we have just shown is called **Spectral decomposition**, the similar matrix $D$ holds all the eigenvalues on its diagonal. And $P$ is orthogonal matrix, which means any of of its two columns are perpendicular. Therefore it could be rewritten as \n",
    "$$\n",
    "A = PDP^{T}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Spectral Decomposition Visualization</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What these plots tells are actually the functions of each decomposed matrix, $P$ and $P^T$ are for rotation because they orthogonal, and $D$ are for scaling because it's diagonal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4MAAAD7CAYAAADQDkHLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA11UlEQVR4nO3de3xU9Z3/8fdkJhcSSEhmJBgIAgEUrIgxSIuAFQaKrttSl7qggogsaqisbUVBBbeLSKrwwGKlasEUae1WuwpbWy1NrZefl4cEgheqQggKyiUmQ7iEBJjM+f1BSaG5Qmbme87M6/l49FEyczLfdybzTvLxnDnHZVmWJQAAAABAXEkwHQAAAAAAEH0MgwAAAAAQhxgGAQAAACAOMQwCAAAAQBxiGAQAAACAOMQwCAAAAABxiGEQAAAAAOIQwyAAAAAAxCGGQQBwmF/84hcaMmSIhgwZIpfL1fjv2267zXQ0IKbRPcAc+hcZLsuyLNMhAABnbt++fbr44ou1d+9e01GAuEL3AHPoX3ixZxAAHOrDDz/URRddZDoGEHfoHmAO/QsvhkEAsKHvfe97crlccrlcSkhIUI8ePXTXXXepoaGhcRt+IQLtc2qf3G63srKyNGrUKK1atUqhUKjFbeke0Lq2+hLu7kn0L9wYBgHAhkpLSzVr1izt2bNHn3/+uebOnaulS5fqqaeeatyGX4hA+5SWlur73/++9uzZo88++0wvv/yyRowYoVmzZunWW29tsi3dA9qnrb6Eu3sS/Qs3hkEAsJnq6mp99tln+sY3vqHu3bsrNzdXd9xxh3r27KkNGzY0bscvRKBtJ/s0YsSIxj4NGzZMDz30kP77v/9bK1eu1Mcff3zatnQPaFtbfYlE9yT6F24MgwBgMyd/8V1yySWNt9XV1emrr75St27dJEmhUEiffPKJBg0aZCQj4BQn+1RQUNDkvquuukqSVFZWdtq2dA9oW1t9CXf3JPoXCQyDAGAzpaWl6tSpk84//3xJJ86cNmPGDCUkJOjmm2+WJG3fvl3nnnuuUlNTTUYFbK+0tFSZmZnKy8trcl9iYuJp/0/3gPZrqy/h7p5E/yKBYRAAbGbDhg2qr69XRkaGUlNT1bNnT+3evVuvv/564y/K/v37a+vWrYaTAva3YcMGXXrppc3e9+mnn0pS414Guge0X1t9CXf3JPoXCVxnEABspkePHrrqqqs0d+5cJSUlKScnRx6Px3QswJF69OihKVOmqKioqMl9EydO1Mcff6wtW7Y0bkv3gPZpqy90zxl4lgHARvbs2aPdu3fL7/erX79+puMAjnayT829Z+npp5/WCy+8oJdeeum0beke0La2+kL3nINhEABs5OSb6IcNG2Y4CeB8J/uUm5urvXv36siRI/r444/1zDPP6MUXX9Tjjz+uq6+++rRt6R7Qtrb6Qvecg2EQAGyktLRUPp9Pffr0MR0FcLzS0lJJ0te//nUlJCSoS5cu6tOnj0aNGqUPP/zwtPci0T2g/drqC91zDt4zCAAAAABxiLOJAgAAAEAcYhgEAAAAgDhki/cMVlVV6fHHH1dNTY1cLpf8fn/jm0pPsixLxcXFKisrU3JysgoLC9W3b19DiQFnomsAAAA4yRbDoNvt1pQpU9S3b1/V1dVp7ty5Gjx4sHr27Nm4TVlZmfbu3avly5dr27ZtWrlypR566CGDqQHnoWsAAAA4yRaHiWZmZjbueejUqZN69OihQCBw2jalpaUaNWqUXC6XBgwYoNraWu3fv99EXMCx6BoAAABOssWewVNVVlZqx44dTS46GQgE5PP5Gj/2er0KBALKzMxs8hglJSUqKSmRJBUVFUU2MOBQHe0aPQMAAHA2Ww2D9fX1Wrp0qaZNm6bU1NTT7mvuChgul6vZx/H7/fL7/Y0f7969O7xBI8zn86mqqsp0jDPmxNxOzCxJOTk5Hfr8cHSNnpnhxNxOzCx1vGeREI2e2eH7ZTqD6fXtkCFa69uxZ1J8dM30+nbIYHr9aGZoqWu2OExUkoLBoJYuXaqRI0dq2LBhTe73er2nPVHV1dXN7hUE0Dq6BgAAAMkmw6BlWXriiSfUo0cPXXPNNc1uU1BQoDfeeEOWZWnr1q1KTU3lD1TgDNE1AAAAnGSLw0Q//fRTvfHGG+rVq5fmzJkjSZo8eXLj3olx48bpkksu0aZNmzR79mwlJSWpsLDQZGTAkegaAAAATrLFMHjBBRfoueeea3Ubl8ulGTNmRCkREJvoGgAAAE6yxWGiAAAAAIDoYhgEAAAAgDjEMAgAAAAAcYhhEAAAAADikC1OIAMAQDyqra3VE088oV27dsnlcun222/XgAEDTMcCYgo9A1rGMAgAgCHFxcUaMmSIfvSjHykYDOro0aOmIwExh54BLeMwUQAADDhy5Ig+/vhjjR49WpLk8XiUlpZmOBUQW+gZ0Dr2DAIAYEBlZaXS09O1YsUKff755+rbt6+mTZumlJSUxm1KSkpUUlIiSSoqKpLP54t4Lo/HE5V17JzB9Pp2yGB6/XBpT8+k+Oya6fXtkMH0+nbIwDAIAIABDQ0N2rFjh6ZPn67+/furuLhYa9eu1aRJkxq38fv98vv9jR9XVVVFPJfP54vKOnbOYHp9O2SI1vo5OTkRffz29EyKz66ZXt8OGUyvH80MLXWNw0QBADDA6/XK6/Wqf//+kqSvf/3r2rFjh+FUQGyhZ0DrGAYBADCga9eu8nq92r17tyTpww8/VM+ePQ2nAmILPQNax2GiAAAYMn36dC1fvlzBYFDdunVTYWGh6UhAzKFnQMsYBgEAMKR3794qKioyHQOIafQMaBmHiQIAAABAHLLNnsEVK1Zo06ZNysjI0NKlS5vcv2XLFj388MPq1q2bJGnYsGGaOHFitGMCjkbPAAAAcJJthsFvfvObGj9+vB5//PEWtxk4cKDmzp0bxVRAbKFnAAAAOMk2h4kOGjRInTt3Nh0DiGn0DAAAACfZZs9ge2zdulVz5sxRZmampkyZotzc3Ga3KykpUUlJiSSpqKhIPp8vmjE7zOPxOC6z5MzcTswcafTM3pyY24mZAQCIB44ZBvv06aMVK1YoJSVFmzZt0iOPPKLly5c3u63f75ff72/8uKqqKloxw8Ln8zkus+TM3E7MLEk5OTkReVx6Zn9OzO3EzFLkegYAgF3Y5jDRtqSmpiolJUWSlJ+fr4aGBh08eNBwKiC20DMAAID44ZhhsKamRpZlSZLKy8sVCoXUpUsXw6mA2ELPAAAA4odtDhN99NFH9be//U2HDh3Sbbfdpuuuu07BYFCSNG7cOL377rtav3693G63kpKSdOedd8rlchlODTgLPQMAAMBJthkG77zzzlbvHz9+vMaPHx+dMECMomcAAAA4yTGHiQIAAAAAwodhEAAAAADiEMMgAAAAAMQhhkEAAAAAiEMMgwAAAAAQhxgGAQAAACAOMQwCAAAAQBxiGAQAwJBQKKS7775bRUVFpqMAMY2uAc1jGAQAwJA//vGP6tGjh+kYQMyja0DzGAYBADCgurpamzZt0pgxY0xHAWIaXQNa5jEdAACAePTLX/5SN954o+rq6lrcpqSkRCUlJZKkoqIi+Xy+iOfyeDxRWcfOGUyvb4cMptcPJ7pmz/XtkMH0+nbIwDAIAECUbdy4URkZGerbt6+2bNnS4nZ+v19+v7/x46qqqohn8/l8UVnHzhlMr2+HDNFaPycnJ6KPT9fsu74dMpheP5oZWuoawyAAAFH26aefqrS0VGVlZTp27Jjq6uq0fPlyzZ4923Q0IKbQNaB1thkGV6xYoU2bNikjI0NLly5tcr9lWSouLlZZWZmSk5NVWFiovn37GkgKOBc9A+zh+uuv1/XXXy9J2rJli37/+9/zxykQAXQNaJ1tTiDzzW9+U/fee2+L95eVlWnv3r1avny5Zs6cqZUrV0YxHRAb6BkAAABOss2ewUGDBqmysrLF+0tLSzVq1Ci5XC4NGDBAtbW12r9/vzIzM6OYEnA2egbYz4UXXqgLL7zQdAwg5tE1oCnbDINtCQQCp51px+v1KhAINPtHqokzQoWT6bMKnS0n5nZi5kiiZ/b2+eeS1+u83E58rgEAiAeOGQYty2pym8vlanZbE2eECic7nNnobDgxtxMzS5E7+xo9s7c1a9I0YkSqvvY1Z+V24nMtRf4shwAAmOaYYdDr9Z72x0R1dTWHrgFhRs/sraQkRQcPJuhrXzOdBAAAxALbnECmLQUFBXrjjTdkWZa2bt2q1NRU/kgNk2DQdALYBT2zr5oal957L0l/+INjfmwDAACbs82ewUcffVR/+9vfdOjQId1222267rrrFPz7lDJu3Dhdcskl2rRpk2bPnq2kpCQVFhYaThw7li/vrP/8z8Nyu00nQaTRM+f6619T1NDgUnm5VF7uVr9+DaYjAQAAh7PNMHjnnXe2er/L5dKMGTOiEybObNmSqKKiLrrvvkOmoyDC6Jlz/fnPyY3/LilJUb9+tQbTAACAWMDxRlB2dkgrVnTRCy90Mh0FQDOOH5f27HErO7tBF14Y0gcfJJqOBAAAYgDDIJSdfeJwszlzuur99/kjE7Abl0v63/+tVk5Og0aMsPToozWmIwEAgBjAMAh1735iGKyvd2n69CxVVvKyAOzE45ESTqllUpK5LAAAIHbwVz+UnR1Sp04hSdL06bU6eLD568oBAAAAiB0Mg1BublC/+U1AAwce14cfJnKWQgAAACAOMAxCeXkNGjr0mKZOrdXLL6do3z5eFgAAAECs469+NLr22jqlpFh69tlU01EAAAAARBjDIBp17mxp4sQ6/epXafr7dcgBAAAAxCiGQZxm6tRa7d3r1p//nGI6CgAAAIAIYhjEac4/P6hvfOOoVq9OMx0FAAAAQAQxDKKJqVNr9eabySovd5uOAgAAACBCPKYDwH7Gj69Xt24NWrMmTT/+8UHTcQAgJlVVVenxxx9XTU2NXC6X/H6/rr76atOxgJhCz4DWMQyiiaQk6frrj+jpp9N0zz2HlJpqmY4EADHH7XZrypQp6tu3r+rq6jR37lwNHjxYPXv2NB0NiBn0DGgdh4miWTfcUKvaWpfWretkOgoAxKTMzEz17dtXktSpUyf16NFDgUDAcCogttAzoHW22TO4efNmFRcXKxQKacyYMZowYcJp92/ZskUPP/ywunXrJkkaNmyYJk6caCBpfMjJCelb36rXL3+ZqkmTjsjlMp0I4ULXAPuprKzUjh071K9fv9NuLykpUUlJiSSpqKhIPp8v4lk8Hk9U1rFzBtPr2yGD6fUjoaWeSfHZNdPr2yGD6fXtkMEWw2AoFNKqVat0//33y+v1at68eSooKGiyC3/gwIGaO3euoZTxZ+rUWk2a5FNZWaLy84+bjoMwoGuA/dTX12vp0qWaNm2aUlNTT7vP7/fL7/c3flxVVRXxPD6fLyrr2DmD6fXtkCFa6+fk5ER8Dan1nknx2TXT69shg+n1o5mhpa7Z4jDR8vJyde/eXdnZ2fJ4PBo+fLg2bNhgOlbcGzHimPLyjnOZiRhC1wB7CQaDWrp0qUaOHKlhw4aZjgPEJHoGtMwWewYDgYC8Xm/jx16vV9u2bWuy3datWzVnzhxlZmZqypQpys3NbfbxTOzqDyfTu4tPdfvtLt13Xyf99KcetRXJTrnby4mZOyKcXaNn0ZeY6FFCgstxuZ34XEeDZVl64okn1KNHD11zzTWm4wAxiZ4BrbPFMGhZTc9W6fqnN6n16dNHK1asUEpKijZt2qRHHnlEy5cvb/bxTOzqDyc77LI+6eqrXbr//mz9/Od1uv322la3tVPu9nJiZunsD6sJZ9foWfQdP+5TKOR2XG4nPtdS5A9f+/TTT/XGG2+oV69emjNnjiRp8uTJys/Pj+i6QDyhZ0DrbDEMer1eVVdXN35cXV2tzMzM07Y59fju/Px8rVq1SgcPHlR6enrUcsajjAxL115bp2eeSdOtt9YqwRYHFuNs0TXAPi644AI999xzpmMAMY2eAa2zxZ/2eXl52rNnjyorKxUMBvX222+roKDgtG1qamoa92qUl5crFAqpS5cuJuLGnalTj2jnTo9eey3ZdBR0EF0DAADASbbYM+h2uzV9+nQtWrRIoVBIV155pXJzc7V+/XpJ0rhx4/Tuu+9q/fr1crvdSkpK0p133tnk8DZExkUXHVd+/jGtXp2m0aOPmo6DDqBrAAAAOMkWw6B04nC0fz5+e9y4cY3/Hj9+vMaPHx/tWPi7m26q1Z13dtWuXW7l5jaYjtM2y5Jr/3659+2Te98+JezbJ3dlpeq++101/NNlFOINXQMAAIBkk8NEYX/XXFOnrl1D+tWvml6bx27cO3eq2/DhOveii9TN75f3hhuUsWCBjp9/ftwPggAAAMBJDINol5QUafLkI3r22VQdtfGRoq4jR5T85puykpIabwv27q2q3/9eR0/Z+wUAAADEO4ZBtNuNNx7R/v0J+sMfOpmO0oSnvFzpCxYoOz9fGfPmKXj++WrIzFT9qFH66qWXFBwwwHREAAAAwFZs855B2N955zXoyiuPavXqNF17bZ3pONLx40pZv15pq1cr+a231NC9uw7PnKkjkycrdO656vzYYzp8++2Sh5c5AAAA8M/4Kxln5KabanXTTV599JFHX/ta0EiGhL17lfrss0r79a/l3rtXRy+/XIGnnlL9uHFSYmLjdofvuMNIPgAAAMAJGAZxRq688qhyc4N65pk0PfzwgegtbFlKevttpa1erZRXXpGVmqoj112nI1OmKNi/f/RyAAAAADGCYRBnxO2Wpkw5omXLOuu++w4qI8OK6HqugweV+rvfKfWZZ5S4bZuODxqkA4sXq27CBFlpaRFdGwAAAIhlnEAGZ2zSpCNqaHDpd7+L3GUmPFu2KOPuu5V96aVKX7hQxy+6SF+tXauv1q/XkRtuYBAEAAAAOog9gzhjXm9I11xTp9WrUzV9eq1crjA98NGj6vTHPyp19Wolb9igYM+eOvyf/6kjkyYp5POFaREAAAAAEsMgztJNN9XqO985R2+9laQRI4516LHcX3yh1DVrlPqb3yghENDRK69U9S9/qaOjR584LhUAAABA2DEM4qxceulxXXjhca1enXZ2w2AopOTXXz9xWYiSElkZGToyebJqb7xRDb17hz0vAAAAgNMxDOKsuFzStGm1mjs3Q3v2JOjcc0Pt+7xAQKnPPae0Z56R5/PPdeySS1SzbJnqrrlG6mS/i9kDAAAAsYoTyOCsTZhQp7Q0S88+28bJXCxLiWVl6nrnnepeUKAujzyio8OH66uXX1bVSy+p7nvfYxAEAAAAoow9gzhrqamWvve9I/r1r1M1e/ahJve76uqUsm6d0lavVtIHHyjYp48OzpunI9/7nqyuXaMfGAAAAEAj2wyDmzdvVnFxsUKhkMaMGaMJEyacdr9lWSouLlZZWZmSk5NVWFiovn37mgmLRlOn1mrVqs76059SNG3aidvc27crbc0apT73nFyHDqn+W99S1bx5OjZihJTAzmjT6BoiwXXwoKxOnaTERNNRHKWtPiK20BMz6BnQsnb/Zb569Wp99tlnEQkRCoW0atUq3XvvvVq2bJneeustffHFF6dtU1ZWpr1792r58uWaOXOmVq5cGZEsODP9+jVoxIijWv3LVLnWrVPW5MnKHjVKndauVe3NN2vfu+9q/8qVOjZqFINgO7z44otNXvvhRNcQKa7aWmVfdpm6zpqlTmvXylVTYzpSh0Xy957Uvj4itriOHDnRk8JCdXrxxZjoSUfRM8Csdu8ZbGho0KJFi5Senq6RI0dq5MiR8nq9YQlRXl6u7t27Kzs7W5I0fPhwbdiwQT179mzcprS0VKNGjZLL5dKAAQNUW1ur/fv3KzMzMywZcPZuuqlWj9x9TPXXTVfS17+mwIoVqr/qKikpyXQ0xwmFQnriiSfUuXNnFRQUqKCgQF3DeEgtXXO23ufWqudzT6vbn35mOkqzEmpqlLp2rVLXrpXlduvYZZepftw4afZs09HOSiR/70nt6yOiLxBIUNWM+zTo499F5PETamqUum6dUtet+0dPxo7VkRtukNW5c0TWtDN6Fp+OHpW2/3m3vn7eVumiAabjxLV2D4PTp0/XtGnTVFZWpjfffFMvvPCC+vfvr1GjRmnYsGFKSUk56xCBQOC04nu9Xm3btq3JNr5TLjzu9XoVCASa/QO1pKREJSUlkqSioqLTPs8JPB6PozLfeKM0ZdTnSqx/S8ELLlBnSU75dWa35/rWW29VKBTS+++/r7feektFRUXKy8vTiBEjNHTo0A71TApv1+hZ9P3PM4eU+NMDCoVuMB2leY89Jh07Jis1VdaYMXL/y7+o01VXydO1q3zBoOl0ZyySv/ek9vXRRM/s0A2TGX7zmwQtev1BffWD5Mgs8LOfNduTTt27n7aZ6e9DtNa3Q8+k+OyayfXLylwae2uONnUeqQv/38+kgQON5DD9PbBDhjN6z2BCQoIuvfRSXXrppdq1a5eWL1+uFStWaOXKlbr88st13XXXKSsr64xDWJbV5DaXy3XG25zk9/vl9/sbP66qqjrjTCb5fD7HZVbnNPl6n+e43HZ9rnNzczVp0iRdccUVWrNmjZ566ikVFxcrPz9f48eP16BBg87qccPZNXpmhu/++22Z27Ntm9J27VL92LE6Ony4dMofcL5g0JaZ25KTkxOx33tS+7pmomd26IbJDLW1aVJyuvYVFob9sT3l5UrbufMfPTn1TNr/9PWa/j5Ea3079EyKz66ZXH/TphRJWerXJ6iEf/1XVf3+9wqdc07Uc5j+HkQzQ05OTrO3n9EweOTIEb377rt688039fnnn2vYsGG65ZZb5PP59NJLL+mhhx7SkiVLzjic1+tVdXV148fV1dVN9kJ4vd7TnqjmtgFiQX19vTZv3qzS0lLt3r1bF198sSZOnKjMzEz99a9/1ZNPPqmf/vSnZ/XYdA2REuzfXwcWLzYdI+wi9XtPal8fEVuC/frpQFGR6Ri2Q8/iT0WFR+ee26Dk//uNXMOHK+vmm1X9/PMnTrCEqGr3MLh06VK9//77GjhwoMaOHauhQ4cq8ZSzYU2dOlXTTp5O8gzl5eVpz549qqysVFZWlt5++23N/qf3mBQUFOiVV17R5Zdfrm3btik1NZUyI+YUFxfrk08+UV5enoYPH67BgwfL4/lHTSdMmKB58+ad9ePTNaD9Ivl7T2pfH4FYR8/iU0WFR337BqWcHFWvWSPfhAnqescd2v/kk5LbbTpeXGn3MNi/f3/dcsstLZ7MIiEhQb/4xS/OKoTb7db06dO1aNEihUIhXXnllcrNzdX69eslSePGjdMll1yiTZs2afbs2UpKSlJhBA7fAEw777zz9G//9m9KT09v9v6EhAQtXLjwrB+frgHtF8nfe1LLfQTiCT2LT9u3ezR48HFJiQoOHKj9v/iFsqZMUfrChTr4X/9lOl5cafcw+O1vf7vNbZKTz/7N1vn5+crPzz/ttnHjxjX+2+VyacaMGWf9+IATjB49us1tkjp4lla6BrRPpH/vSc33EYgn9Cz+WNaJPYMTJtRJOrEX+OioUTpQVKSud92lhl69VDt9utmQccQ2F50HAAAAENuqqhJ06FDCicNET3Fk8mS5d+5U+gMPKNizp46e8h+qETlcBRwAAABAVFRUnNgX9c/DoCQduvtu1X3nO8osLFTi++9HO1pcYhgEAAAAEBUVFR55PJZycxua3ulyqWbpUh0fPFhZN90k9xdfRD9gnGEYBAAAABAVFRVu9erVoFNOGnu65GQFVq1SKD1dWVOmyHXgQFTzxRuGQQAAAABRUVHhUV5e00NET2VlZiqwZo0SqquV9R//IR07FqV08YdhEAAAAEBUNF5jsA0N552nQHGxkjZuVNe77z5xGlKEHcMgAAAAgIhraJA++6x9w6AkHb/0Uu1fvlydfvc7dX700ciGi1MMgwAAAAAi7ssv3Tp2zNXuYVCS6v/lX3Rw/nylL1miTs8/H8F08YnrDAIAAACIuNYuK9Ga2pkz5dm5U13nzFFDTo6OXX55JOLFJfYMAgAAAIi47ds9Sk0NKTs7dGaf6HLpwI9/rKNXXKGsGTPk2bo1MgHjEMMgAAAAgIg7efIYl+ssPtnj0f4VKxQ87zxlTZmihMrKsOeLRwyDAAAAACKuosKtvn2budh8O1lpaQqsXi2FQsqaNk2uI0fCmC4+MQwCAAAAiLj2XlaiNaHsbAXWrJGnokJdv//9E6coxVljGAQAAAAQUXV1J84m2tFhUJKCF1ygwFNPKeUvf1H6j38chnTxy/jZRA8fPqxly5bpq6++0jnnnKMf/OAH6ty5c5PtZs2apZSUFCUkJMjtdquoqMhAWsC56BoAADDl8889siyX8vI6PgxK0rFRo1Tz8MPK/OEP1dCrl2pnzAjL48Yb48Pg2rVrddFFF2nChAlau3at1q5dqxtvvLHZbR944AGlp6dHOSEQG+gaAAAw5eRlJfr0Cc8wKEl1//7v8uzcqfT/+i815Oaq/lvfCttjxwvjh4lu2LBBV1xxhSTpiiuu0IYNGwwnAmITXQMAAKZUVHjk8zUoI8MK6+Meuusu1X33u+paWKjEzZvD+tjxwPiewQMHDigzM1OSlJmZqYMHD7a47aJFiyRJY8eOld/vb3G7kpISlZSUSJKKiork8/nCmDjyPB6P4zJLzsztxMxnK9xdo2dmODG3EzNH2po1a7Rx40Z5PB5lZ2ersLBQaWlppmMBMYWe2Us4Th7TLJdLNUuWyLtnj7JuuklVv/+9Gnr1Cv86MSoqw+DChQtVU1PT5PZJkyad0WNkZWXpwIEDevDBB5WTk6NBgwY1u63f7z/tD9iqqqozzmySz+dzXGbJmbmdmFmScnJymr09ml2jZ2Y4MbcTM0st9ywcBg8erOuvv15ut1u/+tWv9OKLL7Z42DaAs0PP7KWiwh229ws2kZyswMqV8n3nO8qaOlVVa9fKvXevghdcEJn1YkhUhsH58+e3eF9GRob279+vzMxM7d+/v8X3KWVlZTVuP3ToUJWXl7c4DALxiq4BznDxxRc3/nvAgAF69913DaYBYhM9s5ft2z0aO/ZoxB7f6tpVgTVr5PvXf1XWjBly79un6mefVUNubsTWjAXGDxMtKCjQ66+/rgkTJuj111/X0KFDm2xTX18vy7LUqVMn1dfX64MPPtDEiRMNpAWci64B9vTqq69q+PDhzd5n4nBsOxzWazJDWtqJ0ynE83Ngh/XDrbWeSfHZtWiuHwhIgYBbQ4Z0ks+XErkMmZmyFixQ8uzZkiTfr3+thiVLWtzc9PfADhmMD4MTJkzQsmXL9Oqrr8rn8+mHP/yhJCkQCOjJJ5/UvHnzdODAAS35+zeyoaFBI0aM0JAhQwymBpyHrgHR1dph2yf/Y8wLL7wgt9utkSNHNvsYJg7HtsNhvSYz1NamSUqP6+cgmut39HDscPRMis+uRXP9TZsSJZ0jn2+/qqr+cahouDMk7N2rjD/8oXHAcT39tKpvv11WRkaz25v+HkQzQ0tdMz4MdunSRQsWLGhye1ZWlubNmydJys7O1iOPPBLtaEBMoWtAdLV22LYkvfbaa9q4caMWLFggl8sVpVRAbKFnzlBR4ZHLZem88yL0nsG/C3Xvrv2rVqnu5ZeVcf/9cu/bp9Rnn1Xt7bdHdF0nM35pCQAA4s3mzZu1bt063XPPPUpOTjYdB4hJ9Mw+Kio8ys1tUEpK29t2mMul+quvVuVrr6l26lSlPf20dOxYFBZ2JuN7BgEAiDerVq1SMBjUwoULJUn9+/fXzJkzDacCYgs9s4+IXVaiFVZ6ug4sXqyka69V4ief6PjgwVFd3ykYBgEAiLLHHnvMdAQg5tEz+6io8GjYsMidSbQ1x5o5YR7+gcNEAQAAAESEZZ24xmC09wyifRgGAQAAAETE3r0JqqtLUN++DaajoBkMgwAAAAAioqLixLvS2DNoTwyDAAAAACKiosKj5GRLOTnsGbQjhkEAAAAAEbF9u0e9ewfldptOguYwDAIAAACICBOXlUD7MQwCAAAAiIiKCo/y8hgG7YphEAAAAEDYHT8u7dzJZSXsjGEQAAAAQNjt3OlWQ4OLy0rYGMMgAAAAgLDjshL25zEd4J133tHzzz+vL7/8Ug899JDy8vKa3W7z5s0qLi5WKBTSmDFjNGHChOgGBRyOrgEAgGiqqPAoIyOkrKyQ6ShogfE9g7m5ubrrrrs0cODAFrcJhUJatWqV7r33Xi1btkxvvfWWvvjiiyimBJyPrgEAgGg6eSZRl8t0ErTE+J7Bnj17trlNeXm5unfvruzsbEnS8OHDtWHDhnZ9LoAT6BoAAIimigqP+vThEFE7Mz4MtkcgEJDX62382Ov1atu2bS1uX1JSopKSEklSUVGRfD5fxDOGk8fjcVxmyZm5nZg5ks6ka/TMDCfmdmJmAEDHVVR4NHz4UdMx0IqoDIMLFy5UTU1Nk9snTZqkoUOHtvn5lmU1uc3Vyv5mv98vv9/f+HFVVVX7gtqEz+dzXGbJmbmdmFmScnJymr09ml2jZ2Y4MbcTM0st9wwA0LbaWpf27uWyEnYXlWFw/vz5Hfp8r9er6urqxo+rq6uVmZnZ0VhAzKFrAADADnbscEsSF5y3OeMnkGmPvLw87dmzR5WVlQoGg3r77bdVUFBgOhYQc+gaAAAIh+3bT+xz6tOHawzamfFh8L333tNtt92mrVu3qqioSIsWLZJ04r1LixcvliS53W5Nnz5dixYt0g9+8AN94xvfUG5ursnYgOPQNcB+/u///k/XXXedDh48aDoKELPomRkVFR51796gtLSmb0GBfRg/gcxll12myy67rMntWVlZmjdvXuPH+fn5ys/Pj2Y0IKbQNcBeqqqq9OGHH3JyHSCC6Jk5Jy8rAXszvmcQAIB4tHr1at1www2tnhANQMfQM3N27GAYdALjewYBAIg3paWlysrKUu/evVvdzsQlXOxwKRCTGdLSTvx38nh+Duywfji0t2dSfHYtkutblrRjR6ImTUpodY1Yfg6ckoFhEACACGjtUi8vvvii7r///jYfw8QlXOxwKRCTGWpr0ySlx/VzEM31O3oJl3D0TIrPrkVy/erqBNXUdFd29gFVVbV8ncFYfg7slqGlrjEMAgAQAS1d6mXnzp2qrKzUnDlzJJ24hMs999yjxYsXq2vXrlFMCDgfPbOniooTIwaHidofwyAAAFHUq1cvrVy5svHjWbNmafHixUpPTzeYCogt9Mysigq33G5LvXpxWQm74wQyAAAAAMKmosKjXr0alJhoOgnawp5BAAAMevzxx01HAGIePYuuigqP8vI4RNQJ2DMIAAAAIGy2b+eyEk7BMAgAAAAgLBoapM8+Yxh0CoZBAAAAAGGxe7dbR4+6GAYdgmEQAAAAQFhwWQlnYRgEAAAAEBYVFW516hRS9+4h01HQDgyDAAAAAMKiosKjvn0b5HKZToL2YBgEAAAAEBYnhkEOEXUK49cZfOedd/T888/ryy+/1EMPPaS8vLxmt5s1a5ZSUlKUkJAgt9utoqKiKCcFnI2uAQCASKuo8Oi7360zHQPtZHwYzM3N1V133aWnnnqqzW0feOABpaenRyEVEHvoGgAAiKSjR6Vdu9zsGXQQ48Ngz549TUcA4gJdAwAAkVJb69KXX7plWS7l5TEMOoXxYfBMLFq0SJI0duxY+f3+FrcrKSlRSUmJJKmoqEg+ny8q+cLF4/E4LrPkzNxOzBwN7ekaPTPDibmdmBkAcGZWr07Tm28mSZLefz9RKSmWBg5kKLS7qAyDCxcuVE1NTZPbJ02apKFDh7b7MbKysnTgwAE9+OCDysnJ0aBBg5rd1u/3n/YHbFVV1VnlNsXn8zkus+TM3E7MLEk5OTnN3h7NrtEzM5yY24mZpZZ7BgBo6rzzglq06MRbTJYv76J33tlnOBHaIyrD4Pz58zv8GFlZWZKkjIwMDR06VOXl5S0Og0C8omsAAMCEPn3+sRdw5szDSk42GAbt5ohLS9TX16uurq7x3x988IF69eplOBUQe+gaAAA4G336NEiSMjJCuvHGI4bToL2Mv2fwvffe09NPP62DBw+qqKhIvXv31n333adAIKAnn3xS8+bN04EDB7RkyRJJUkNDg0aMGKEhQ4aYDQ44DF0DAACR0qmTpe7dG/Tv/35EnTtbpuOgnYwPg5dddpkuu+yyJrdnZWVp3rx5kqTs7Gw98sgj0Y4GxBS6BgAAImnQoOO65ZZa0zFwBhxxmCgAAECkFRenasuWRAWD0s9+1lm7drlNRwIc5Z57DsrrDZmOgTNgfM8gAADx6OWXX9Yrr7wit9ut/Px83XjjjaYjxb3kZOm551IlSatXp2rWrMOGE6Gj6Fl0fe1rXErCaRgGAQCIso8++kilpaVasmSJEhMTdeDAAdORIGnMmPrGf/v9R+VyGQyDDqNnQNs4TBQAgChbv369vvOd7ygxMVHSiUu5wLzs7JCGDDkmSRo7tr6NrWF39AxoG3sGAQCIsj179uiTTz7R//zP/ygxMVFTpkxRv379mmxXUlKikpISSVJRUZF8Pl/Es3k8nqisY9cM3/52grZutfTtb3dRSkoXIxkk898H0+uHQ3t7JsVn10yvb4cMpte3QwaGQQAAImDhwoWqqalpcvukSZMUCoV0+PBhLVq0SNu3b9eyZcv0s5/9TK5/Oi7R7/fL7/c3flxVVRXp2PL5fFFZx64ZLr/co7Iynw4frtJhg28ZNP19iNb6OTk5Hfr8cPRMis+umV7fDhlMrx/NDC11jWEQAIAImD9/fov3rV+/XsOGDZPL5VK/fv2UkJCgQ4cOKT09PYoJ0ZwLLwzq1lsbTMdAO9EzoGN4zyAAAFE2dOhQffTRR5Kk3bt3KxgMqksXc4ck4h9cLmnsWC6YHQvoGdA29gwCABBlo0eP1ooVK/SjH/1IHo9Hs2bNavbQNQBnj54BbWMYBAAgyjwej2bPnm06BhDT6BnQNg4TBQAAAIA4xDAIAAAAAHGIYRAAAAAA4hDDIAAAAADEIeMnkFmzZo02btwoj8ej7OxsFRYWKi0trcl2mzdvVnFxsUKhkMaMGaMJEyZEPyzgYHQNAAAApzK+Z3Dw4MFaunSplixZonPPPVcvvvhik21CoZBWrVqle++9V8uWLdNbb72lL774wkBawLnoGgAAAE5lfBi8+OKL5Xa7JUkDBgxQIBBosk15ebm6d++u7OxseTweDR8+XBs2bIh2VMDR6BoAAABOZfww0VO9+uqrGj58eJPbA4GAvF5v48der1fbtm1r8XFKSkpUUlIiSSoqKlJOTk74w0aYEzNLzsztxMwdFY6u0TNznJjbiZntKFrPox2+X6YzmF7fDhlMr29SvHTN9Pp2yGB6fdMZorJncOHChfrRj37U5H+n7nF44YUX5Ha7NXLkyCafb1lWk9tcLleL6/n9fhUVFamoqEhz584NzxcRRU7MLDkztxMzSy3njmbX6JkZTsztxMySc3N3lB2+btMZTK9vhwym148Hpp9j0+vbIYPp9e2QISp7BufPn9/q/a+99po2btyoBQsWNPuHp9frVXV1dePH1dXVyszMDHtOwOnoGgAAANrL+HsGN2/erHXr1umee+5RcnJys9vk5eVpz549qqysVDAY1Ntvv62CgoIoJwWcja4BAADgVMbfM7hq1SoFg0EtXLhQktS/f3/NnDlTgUBATz75pObNmye3263p06dr0aJFCoVCuvLKK5Wbm9uux/f7/ZGMHxFOzCw5M7cTM0tnlzuSXXPi8+jEzJIzczsxs+Tc3B1lh6/bdAbT69shg+n144Hp59j0+nbIYHp9O2RwWc29SQgAAAAAENOMHyYKAAAAAIg+hkEAAAAAiEPG3zMYaWvWrNHGjRvl8XiUnZ2twsJCpaWlmY7VpnfeeUfPP/+8vvzySz300EPKy8szHalFmzdvVnFxsUKhkMaMGaMJEyaYjtSmFStWaNOmTcrIyNDSpUtNx2mXqqoqPf7446qpqZHL5ZLf79fVV19tOlYjJ3aNnkUWPXOG5557Tn/5y1+Unp4uSZo8ebLy8/ObbBfJ12B7f37MmjVLKSkpSkhIkNvtVlFRUYfWbetrsixLxcXFKisrU3JysgoLC9W3b98OrXlSe15rW7Zs0cMPP6xu3bpJkoYNG6aJEyeGZf2T2npOI/kcxBvTXYvHnkn26Jqte2bFuM2bN1vBYNCyLMtas2aNtWbNGsOJ2mfXrl3Wl19+aT3wwANWeXm56TgtamhosL7//e9be/futY4fP27ddddd1q5du0zHatOWLVus7du3Wz/84Q9NR2m3QCBgbd++3bIsyzpy5Ig1e/ZsWz3XTuwaPYsseuYMv/3tb61169a1uk2kX4Pt/flRWFhoHThwICxrtudr2rhxo7Vo0SIrFApZn376qTVv3rywrG1Z7XutffTRR9bixYvDtmZz2npOI/kcxBvTXYvHnlmWPbpm557F/GGiF198sdxutyRpwIABCgQChhO1T8+ePZWTk2M6RpvKy8vVvXt3ZWdny+PxaPjw4add4NyuBg0apM6dO5uOcUYyMzMb/ytRp06d1KNHD1u9np3YNXoWWfQsdkT6NWji50d7vqbS0lKNGjVKLpdLAwYMUG1trfbv3x+W9Z3yWovkc4CmItm1eOyZ5IyumexZzB8meqpXX31Vw4cPNx0jpgQCAXm93saPvV6vtm3bZjBRfKisrNSOHTvUr18/01GaRdfCi56ZYfeehdOf/vQnvfHGG+rbt6+mTp3aZIiP5muwrZ8fixYtkiSNHTu2Q6dkb8/XFAgE5PP5TtsmEAgoMzPzrNdtTmuvta1bt2rOnDnKzMzUlClT2n1prTPR2nMarecgXtila/HYM8ls1+zas5gYBhcuXKiampomt0+aNElDhw6VJL3wwgtyu90aOXJklNO1rD257c5q5sokLpfLQJL4UV9fr6VLl2ratGlKTU2N6tpO7Bo9w9kw2bNIaK0H48aNa3xvzG9/+1s988wzKiwsPG27cLwGw/HzY+HChcrKytKBAwf04IMPKicnR4MGDTqjHCe152uKRvdae6316dNHK1asUEpKijZt2qRHHnlEy5cvD+v6bT2n/Pw5M6a7Rs9aZrJrdu5ZTAyD8+fPb/X+1157TRs3btSCBQts9QOsrdxO4PV6VV1d3fhxdXU1/7UwgoLBoJYuXaqRI0dq2LBhUV/fiV2jZzhTpnsWCe3twZgxY/STn/ykye3heA2G4+dHVlaWJCkjI0NDhw5VeXn5Wf+R2p6vyev1qqqqqtVtOqKt19qpf7Dm5+dr1apVOnjwYOMJSMKhrec00s9BrDHdNXrWPNNds3PPYv49g5s3b9a6det0zz33KDk52XScmJOXl6c9e/aosrJSwWBQb7/9tgoKCkzHikmWZemJJ55Qjx49dM0115iO0wRdixx6Fj1271kknPq+lPfee6/ZQ6Mi/Rpsz8+P+vp61dXVNf77gw8+UK9evc56zfZ8TQUFBXrjjTdkWZa2bt2q1NTUsP2B1p7XWk1NTeMeg/LycoVCIXXp0iUs60vte04j+RzEG9Ndi8eeSea7Zveeuazm9kvGkDvuuEPBYLDxmOz+/ftr5syZhlO17b333tPTTz+tgwcPKi0tTb1799Z9991nOlazNm3apNWrVysUCunKK6/UtddeazpSmx599FH97W9/06FDh5SRkaHrrrtOo0ePNh2rVZ988okWLFigXr16Nf7XvJZOS22CE7tGzyKLnjnDY489ps8++0wul0vnnHOOZs6cqczMTAUCAT355JOaN2+epMi+Blv6+XFqhn379mnJkiWSpIaGBo0YMaLDGZr7mtavXy9JGjdunCzL0qpVq/T+++8rKSlJhYWFYbsETUuvtZN7B8aNG6dXXnlF69evl9vtVlJSkqZOnarzzz8/LOtLavE5jdZzEG9Mdy0eeyaZ75rdexbzwyAAAAAAoKmYP0wUAAAAANAUwyAAAAAAxCGGQQAAAACIQwyDAAAAABCHGAYBAAAAIA4xDAIAAABAHGIYBAAAAIA4xDAIAAAAAHGIYRBnbe/evbr55ptVUVEhSQoEArrlllu0ZcsWw8mA2ELXgMijZ0Dk0TP7YRjEWevevbtuuOEGPfbYYzp69Kh+/vOf64orrtCFF15oOhoQU+gaEHn0DIg8emY/LsuyLNMh4Gw/+clPVFlZKZfLpcWLFysxMdF0JCAm0TUg8ugZEHn0zD7YM4gOGzNmjHbt2qXx48dTZiCC6BoQefQMiDx6Zh8Mg+iQ+vp6rV69WqNHj9bzzz+vw4cPm44ExCS6BkQePQMij57ZC8MgOqS4uFh9+vTRbbfdpvz8fD311FOmIwExia4BkUfPgMijZ/bCMIiztmHDBm3evFkzZ86UJN10003asWOH3nzzTcPJgNhC14DIo2dA5NEz++EEMgAAAAAQh9gzCAAAAABxiGEQAAAAAOIQwyAAAAAAxCGGQQAAAACIQwyDAAAAABCHGAYBAAAAIA4xDAIAAABAHGIYBAAAAIA49P8B91HjoipHnDgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x288 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Define matrix A\n",
    "A = np.array([[1, 3], [3, -5]])\n",
    "\n",
    "# Calculate eigenvalues and eigenvectors\n",
    "eigenvalues, eigenvectors = np.linalg.eig(A)\n",
    "\n",
    "# Plot eigenvectors\n",
    "fig, ax = plt.subplots(figsize=(15, 4), nrows=1, ncols=4)\n",
    "for i in range(2):\n",
    "    ax[0].quiver(0, 0, eigenvectors[:,i][0], eigenvectors[:,i][1], angles='xy', \n",
    "                 scale_units='xy', scale=1, color=['r','b'][i])\n",
    "ax[0].set_xlim(-2, 2)\n",
    "ax[0].set_ylim(-2, 2)\n",
    "\n",
    "for i in range(4):\n",
    "    ax[i].set_xlabel('x')\n",
    "    ax[i].set_ylabel('y')\n",
    "    ax[i].set_aspect('equal')\n",
    "\n",
    "ax[1].quiver(0, 0, 1, 0, angles='xy', \n",
    "             scale_units='xy', scale=1, color=['r'])\n",
    "ax[1].quiver(0, 0, 0, 1, angles='xy', \n",
    "             scale_units='xy', scale=1, color=['b'])\n",
    "ax[1].set_xlim(-2, 2)\n",
    "ax[1].set_ylim(-2, 2)\n",
    "ax[1].set_title('$P^T$')\n",
    "\n",
    "\n",
    "ax[2].quiver(0, 0, eigenvalues[0], 0, angles='xy', \n",
    "             scale_units='xy', scale=1, color=['r'])\n",
    "ax[2].quiver(0, 0, 0, eigenvalues[1], angles='xy', \n",
    "             scale_units='xy', scale=1, color=['b'])\n",
    "ax[2].set_xlim(-7, 7)\n",
    "ax[2].set_ylim(-7, 7)\n",
    "ax[2].set_title('$DP^T$')\n",
    "\n",
    "temp = np.array([[eigenvalues[0], 0],\n",
    "                 [0, eigenvalues[1]]])\n",
    "temp1 = temp@eigenvectors\n",
    "\n",
    "ax[3].quiver(0, 0, temp1[:,0][0], temp1[:,0][1], angles='xy', \n",
    "             scale_units='xy', scale=1, color=['r'])\n",
    "ax[3].quiver(0, 0, temp1[:,1][0], temp1[:,1][1], angles='xy', \n",
    "             scale_units='xy', scale=1, color=['b'])\n",
    "ax[3].set_xlim(-7, 7)\n",
    "ax[3].set_ylim(-7, 7)\n",
    "ax[3].set_title('$PDP^T$')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, in reality we don't have many chances of working on symmetric matrices, so spectral decomposition is mostly for theoretical demonstration, not so much application in reality."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Diagonalizing a Matrix</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a matrix\n",
    "\n",
    "$$A=\\left[\\begin{array}{rrr}\n",
    "1 & 3 & 3 \\\\\n",
    "-3 & -5 & -3 \\\\\n",
    "3 & 3 & 1\n",
    "\\end{array}\\right]$$\n",
    "\n",
    "We seek to diagonalize the matrix $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following these steps:\n",
    "\n",
    "1. Compute the eigenvalues of $A$\n",
    "2. Compute the eigenvectors of $A$\n",
    "3. Construct $P$.\n",
    "4. Construct $D$ from the corresponding columns of $P$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[ \\left( -2, \\  2, \\  \\left[ \\left[\\begin{matrix}-1\\\\1\\\\0\\end{matrix}\\right], \\  \\left[\\begin{matrix}-1\\\\0\\\\1\\end{matrix}\\right]\\right]\\right), \\  \\left( 1, \\  1, \\  \\left[ \\left[\\begin{matrix}1\\\\-1\\\\1\\end{matrix}\\right]\\right]\\right)\\right]$"
      ],
      "text/plain": [
       "⎡⎛       ⎡⎡-1⎤  ⎡-1⎤⎤⎞  ⎛      ⎡⎡1 ⎤⎤⎞⎤\n",
       "⎢⎜       ⎢⎢  ⎥  ⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟⎥\n",
       "⎢⎜-2, 2, ⎢⎢1 ⎥, ⎢0 ⎥⎥⎟, ⎜1, 1, ⎢⎢-1⎥⎥⎟⎥\n",
       "⎢⎜       ⎢⎢  ⎥  ⎢  ⎥⎥⎟  ⎜      ⎢⎢  ⎥⎥⎟⎥\n",
       "⎣⎝       ⎣⎣0 ⎦  ⎣1 ⎦⎦⎠  ⎝      ⎣⎣1 ⎦⎦⎠⎦"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sy.Matrix([[1,3,3], [-3, -5, -3], [3,3,1]])\n",
    "eig = sy.matrices.matrices.MatrixEigen.eigenvects(A)\n",
    "eig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reminder the return value takes the form ```[(eigenval, multiplicity, eigenspace), ...]```."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct $P$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-1 & -1 & 1\\\\1 & 0 & -1\\\\0 & 1 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡-1  -1  1 ⎤\n",
       "⎢          ⎥\n",
       "⎢1   0   -1⎥\n",
       "⎢          ⎥\n",
       "⎣0   1   1 ⎦"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P = sy.zeros(3, 3)\n",
    "P[:, 0] = eig[0][2][0]\n",
    "P[:, 1] = eig[0][2][1]\n",
    "P[:, 2] = eig[1][2][0]\n",
    "P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct $D$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-2 & 0 & 0\\\\0 & -2 & 0\\\\0 & 0 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡-2  0   0⎤\n",
       "⎢         ⎥\n",
       "⎢0   -2  0⎥\n",
       "⎢         ⎥\n",
       "⎣0   0   1⎦"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D = sy.diag(eig[0][0], eig[0][0], eig[1][0])\n",
    "D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can verify if $PDP^{-1}=A$ holds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P * D * P.inv() == A "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course we don't need to go through this process seperately. There is ```diagonalize``` method in SymPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "P, D = A.diagonalize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-1 & -1 & 1\\\\1 & 0 & -1\\\\0 & 1 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡-1  -1  1 ⎤\n",
       "⎢          ⎥\n",
       "⎢1   0   -1⎥\n",
       "⎢          ⎥\n",
       "⎣0   1   1 ⎦"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-2 & 0 & 0\\\\0 & -2 & 0\\\\0 & 0 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡-2  0   0⎤\n",
       "⎢         ⎥\n",
       "⎢0   -2  0⎥\n",
       "⎢         ⎥\n",
       "⎣0   0   1⎦"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We obtain the same results as previous separate steps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes we just want to test if a matrix is diagonalizable, then use ```is_diagonalizable``` in SymPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.is_diagonalizable()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $A$ is symmetric, all of its eigenvectors are orthogonal. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "Av_1 = \\lambda_1v_1 \\quad \\text{and} \\quad Av_2 = \\lambda_2v_2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "v_2 \\cdot Av_1 = \\lambda_1v_1 \\cdot v_2\\\\\n",
    "v_1 \\cdot Av_2 = \\lambda_2v_2 \\cdot v_1\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "v_1 \\cdot Av_2 = v_1 \\cdot A^Tv_2 = v_1 \\cdot \\lambda_2v_2 = \\lambda_2v_1 \\cdot v_2\\\\\n",
    "v_2 \\cdot Av_1 = v_2 \\cdot A^Tv_1 = v_2 \\cdot \\lambda_1v_1 = \\lambda_1v_2 \\cdot v_1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}