phys512/matrices.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "96c50950",
   "metadata": {},
   "source": [
    "# Matrices in numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97abd76c",
   "metadata": {},
   "source": [
    "We're going to be dealing with matrices quite a bit in this section, so let's have a review of how numpy handles matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "076056cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f14947d5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.57487069 0.9823841  0.09112579]\n",
      " [0.85884055 0.89054154 0.77408934]\n",
      " [0.786177   0.65692552 0.72699611]\n",
      " [0.66648673 0.92769796 0.33097987]]\n",
      "[0.85884055 0.89054154 0.77408934]\n",
      "0.7740893379034613\n",
      "0.7740893379034613\n"
     ]
    }
   ],
   "source": [
    "# Matrices are 2D arrays in numpy\n",
    "# A[i][j] = ith row and jth column\n",
    "A = np.random.rand(12).reshape(4,3)\n",
    "print(A)\n",
    "print(A[1])\n",
    "print(A[1][2])\n",
    "print(A[1,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "78a723e8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.57487069, 0.85884055, 0.786177  , 0.66648673],\n",
       "       [0.9823841 , 0.89054154, 0.65692552, 0.92769796],\n",
       "       [0.09112579, 0.77408934, 0.72699611, 0.33097987]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Can also use np.transpose(A)  or  A.T\n",
    "A.transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "01bcc062",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.81301628, 4.96219165, 4.28101637, 3.51482228])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Matrix multiplication uses the @ operator\n",
    "y = np.array((1,2,3))\n",
    "A@y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c5653f87",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 0, 0, 0],\n",
       "       [0, 2, 0, 0],\n",
       "       [0, 0, 3, 0],\n",
       "       [0, 0, 0, 4]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Diagonal matrix\n",
    "np.diag((1,2,3,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7f4e7ee1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.        , 0.        , 0.        , 0.        ],\n",
       "       [0.        , 0.5       , 0.        , 0.        ],\n",
       "       [0.        , 0.        , 0.33333333, 0.        ],\n",
       "       [0.        , 0.        , 0.        , 0.25      ]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Inverse\n",
    "np.linalg.inv((np.diag((1,2,3,4))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "11a0999c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.00000000e+00, -6.01445136e-16],\n",
       "       [-6.60539330e-16,  1.00000000e+00]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.random.rand(4).reshape(2,2)\n",
    "B = np.linalg.inv(A)\n",
    "B@A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7f3959e8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.10191999972462222\n",
      "-0.10191999972462229\n"
     ]
    }
   ],
   "source": [
    "# Determinant\n",
    "print(np.linalg.det(A))\n",
    "print(A[0,0]*A[1,1] - A[0,1]*A[1,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "826cced3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "23.999999999999993"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# should be 1 * 2 * 3 * 4 = 24\n",
    "np.linalg.det(np.diag((1,2,3,4)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e619e0da",
   "metadata": {},
   "source": [
    "# Matrix decompositions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "aa6ee144",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We need scipy.linalg for eigenvalues\n",
    "import scipy.linalg\n",
    "\n",
    "# Visualize matrices as a color map\n",
    "def plot_matrices(A,titles=[]):\n",
    "    plt.clf()\n",
    "    n = len(A)\n",
    "    if titles==[]:\n",
    "        titles = [\"\"]*n\n",
    "    plt.figure(figsize=(n*4,4))\n",
    "    for i,AA in enumerate(A):\n",
    "        plt.subplot(1, n, i+1)\n",
    "        plt.imshow(AA)\n",
    "        plt.colorbar()\n",
    "        plt.title(titles[i])\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e2c140b",
   "metadata": {},
   "source": [
    "## LU decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f0f93712",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9oAAAFWCAYAAACW+uXNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABSt0lEQVR4nO3dfXhU1aHH+98kIRMQAiLkDVODVgUECSaSG7BqDzlGpZxyjseiUsFUsQqpSDzWoJCgVKKtYlSQCJWiVyj0DWsLDaWxaC2pSDCt9PBSRCEHnQCXSiAIgZm5f9CZYZwZyN7zmpnv53n28zQ7e629BvDXtWatvbbF6XQ6BQAAAAAAQiIp2g0AAAAAACCeMNAGAAAAACCEGGgDAAAAABBCDLQBAAAAAAghBtoAAAAAAIQQA20AAAAAAEKIgTYAAAAAACHEQBsAAAAAgBBioA0AAAAAQAgx0AYAAAAAIIQYaAOImnfeeUfjxo1TTk6OLBaL3njjjXOW2bBhg6666ipZrVZ99atf1bJly8LeTgCINPIRALo2BtoAoqa9vV3Dhw/XwoULO3X9xx9/rLFjx+rrX/+6mpub9eCDD+qee+7RunXrwtxSAIgs8hEAujaL0+l0RrsRAGCxWLR69WqNHz8+4DWPPPKI1qxZo61bt7rP3Xbbbfr8889VX18fgVYCQOSRjwDQ9TCjDaDLaGxsVElJide50tJSNTY2RqlFABAbyEcAiC0p0W4AgNhy/PhxdXR0mCrrdDplsVi8zlmtVlmt1lA0TTabTZmZmV7nMjMz1dbWpi+++ELdu3cPyX0AwJ9g8lEKb0aSjwCiLZiMTE1NVVpaWohbFF0MtAG4HT9+XAMv6inbfrup8j179tTRo0e9zlVXV2vOnDkhaB0ARE+w+SiRkQDiV7AZmZWVpY8//jiuBtsMtAG4dXR0yLbfro+bLlJ6L2NPlrQdcWhgwR61tLQoPT3dfT5Us9nS6RBubW31Otfa2qr09HRmawCEVTD5KIU/I8lHANEUij5kR0cHA20A8e28nqcPI+z/2lYxPT3dqxMZSsXFxVq7dq3XufXr16u4uDgs9wOALzOTj1L4M5J8BBALgulDxhs2QwPgwyGnqcOoo0ePqrm5Wc3NzZJOv56mublZe/fulSTNnDlTkyZNcl9/3333affu3fr+97+v7du366WXXtLPfvYzzZgxIySfGwDOxWw+Gs1I8hFAVxSpPmRXwIw2AB8OOeQwUcaozZs36+tf/7r754qKCknS5MmTtWzZMn322WfuTqUkDRw4UGvWrNGMGTP0/PPP68ILL9SPf/xjlZaWGr43AJhhJh9d5YwgHwF0RZHqQ3YFvEcbgFtbW5t69+6tT3dcaOr5mpzL/0+HDx8O29JxAIiWYPJRIiMBxDf6kL6Y0Qbgw+50ym7wOzij1wNAV2QmH13lACDe0Yf0YKANwIeZ52Xi9fkaADiT2ecJyUgAiYA+pAcDbQA+HHLKTkgCgA8z+egqBwDxjj6kBwNtAD74NhIA/GNGGwACow/pwUAbgA+erwEA/3hGGwACow/pwXu0AQAAAAAIIWa0Afhw/OswWgYA4p2ZfHSVA4B4Rx/Sg4E2AB92ExtZmNkcCAC6GjP56CoHAPGOPqQHA20APuzO04fRMgAQ78zko6scAMQ7+pAeDLQB+GDZDwD4x9JxAAiMPqQHA20APhyyyC6L4TIAEO/M5KOrHADEO/qQHuw6DgAAAABACDGjDcCHw3n6MFoGAOKdmXx0lQOAeEcf0oOBNgAfdhPLfswspQSArsZMPrrKAUC8ow/pwUAbgA9CEgD8Y6ANAIHRh/RgoA3Ah8NpkcNpcCMLg9cDQFdkJh9d5QAg3tGH9GCgDcAH30YCgH/MaANAYPQhPdh1HAEtW7ZMFovFfaSlpemyyy5TeXm5Wltbo908AIgaVz5u3rw52k0BgJgxZ84cWSwWHTx40O/vhw4dquuvvz6yjQKihBltnNMTTzyhgQMH6vjx43r33Xe1aNEirV27Vlu3blWPHj2i3TyEgV1Jshv8Hs4eprYAQCwxk4+nywFA/KMP6cFAG+d00003qbCwUJJ0zz336IILLtD8+fP161//WrfffnuUW4dwcJp4vsYZp8/XAMCZzOSjqxwAxDv6kB4sHYdh//Zv/yZJ+vjjj6PcEoSL6/kaowcAxDuz+UhGAkgE5KMHM9ow7KOPPpIkXXDBBVFuCcLF7kyS3Wlw2Y8zTI0BgBhiJh9PlwtDYwAgxtCH9GCgjXM6fPiwDh48qOPHj+vPf/6znnjiCXXv3l3f+MY3ot00hIlDFjkMLnhxKE5TEgDOYCYfT5cjIwHEP/qQHgy0cU4lJSVeP1900UVavny5BgwYEKUWAQAAAEDsYqCNc1q4cKEuu+wypaSkKDMzU5dffrmSkni8P57xDkQA8I/3aAPBsVj4byGe0Yf0YLSEcxo5cqRKSkp0/fXXa/DgwQyyE4Dr+RqjBwDEO7P5aCYjFy5cqLy8PKWlpamoqEibNm066/W1tbW6/PLL1b17d+Xm5mrGjBk6fvy42Y8KGJaWliZJ+uKLL/z+/tixY+5rEJ/oQ3rE56cCEJTTz9cYPwAg3pnNR6MZuWrVKlVUVKi6ulpbtmzR8OHDVVpaqv379/u9fsWKFaqsrFR1dbW2bdumV155RatWrdKjjz4aio8NdMpFF10kSdqxY4fP744dO6aWlhb3NYhP9CE9GGgD8OFQkuwGDzObAwFAV2MmH81k5Pz58zVlyhSVlZVpyJAhqqurU48ePbR06VK/12/cuFGjR4/WHXfcoby8PN1www26/fbbzzkLDoTSmDFjlJqaqkWLFsnhcHj9bvHixTp16pRuuummKLUOkUAf0oNntAH4MPdqhvjcMRI4m6VLl6q+vt7n/PTp09WrV68otAjhZv71Xp3PyI6ODjU1NWnmzJnuc0lJSSopKVFjY6PfMqNGjdLrr7+uTZs2aeTIkdq9e7fWrl2rO++803BbAbMyMjJUVVWlWbNm6dprr9V//Md/qEePHtq4caN++tOf6oYbbtC4ceOi3UyEEX1IDwbaAACYtGjRIr/n77rrLgba8Kutrc3rZ6vVKqvV6nXu4MGDstvtyszM9DqfmZmp7du3+633jjvu0MGDB3XNNdfI6XTq1KlTuu+++1g6joh77LHHlJeXpwULFuiJJ57QqVOnNHDgQD3++ON65JFH2OsHCYN/6QjorrvuktPpVGFhYbSbgghz/GsZj9EDSBSufAx0XHjhhdFuIsLEbD66MjI3N1e9e/d2HzU1NSFp14YNGzRv3jy99NJL2rJli371q19pzZo1mjt3bkjqB4yYOHGiGhsbdfToUR0/flzbtm1TVVWVz5dKiD/0IT2Y0Qbgw+60yO40+GoGg9cDQFdkJh9d5SSppaVF6enp7vP+Bh79+vVTcnKyWltbvc63trYqKyvLb/2zZ8/WnXfeqXvuuUeSNGzYMLW3t+vee+/VY489xiwigIigD+lB6gLwYWajHztxAiABmM1HV0amp6d7Hf4G2qmpqSooKFBDQ4P7nMPhUENDg4qLi/2269ixYz6D6eTkZEmSM06ffwQQe+hDejCjDcCHw5kkh8GNLBx05AAkADP5eLqcsYysqKjQ5MmTVVhYqJEjR6q2tlbt7e0qKyuTJE2aNEkDBgxwLz0fN26c5s+frxEjRqioqEi7du3S7NmzNW7cOPeAGwDCjT6kBwNtAD7MfLtoV3yGJACcyezsi9GMnDBhgg4cOKCqqirZbDbl5+ervr7evUHa3r17vWawZ82aJYvFolmzZmnfvn3q37+/xo0bpyeffNJwWwHALPqQHgy0AQAAYlB5ebnKy8v9/m7Dhg1eP6ekpKi6ulrV1dURaBkA4FwYaAPw4ZDxjSkc4WkKAMQUM/noKgcA8Y4+pEfEB9oOh0OffvqpevXqJYslPneYA6LN6XTqyJEjysnJMbXTrJlXLcTrqxkiiXwEIiOYjDT7KhoyMjjkIxAZ9CFDJ+ID7U8//VS5ubmRvi2QkFpaWky9z9fuTJLd4EYWRq+HL/IRiCwzGWkmH13lYB75CERWLPc
      "text/plain": [
       "<Figure size 1200x400 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAFXCAYAAACPydStAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAArc0lEQVR4nO3df1xUZb4H8M+A8kNjxlRgoEZF2wR/Af6IwNZw5Qqu65Wtl1e77Esll34s3EK8lnQTt0xHd1O5lSvaptheSet2tbLWljDh5YqpKKmtopgJoQO6m4PgAjZz7h/GGScB58xwmIfh8369ntc2x/Oc5znZfvnyPc88RyNJkgQiIhKOl7snQEREbWOAJiISFAM0EZGgGKCJiATFAE1EJCgGaCIiQTFAExEJigGaiEhQDNBERIJigCYiEhQDNBHdpqSkBDNmzEBoaCg0Gg127dolxHinTp3Cv/7rv0Kn06Fv376YMGECqqqqVJ2bOzFAE9FtGhsbERkZifXr1wsz3rlz5/DQQw8hPDwc+/btw/Hjx7F06VL4+fl1yRzdQcPNkoioIxqNBjt37kRycrJ8rLm5Gf/1X/+Fd955B1evXsWoUaOwevVqxMfHqzIeAMyZMwe9e/fGn/70J5fH6C6YQRORYhkZGSgtLcX27dtx/PhxzJo1C0lJSTh79qwq41mtVnz88ce4//77kZiYiKCgIMTExKheenE3BmgiUqSqqgpbtmzBe++9h5/+9KcYNmwY/vM//xMPPfQQtmzZosqYdXV1aGhowKpVq5CUlIS//OUv+OUvf4lHHnkExcXFqowpgl7ungARdS8nTpyAxWLB/fffb3e8ubkZAwYMAACcPn0aERERHV7n+eefx6pVqxwa02q1AgBmzpyJhQsXAgCioqJw4MAB5OXl4eGHH1Z6G90CAzQRKdLQ0ABvb2+UlZXB29vb7s/uuusuAMDQoUNx6tSpDq/TGswdMXDgQPTq1QsjRoywOx4REYH9+/c7fJ3uhgGaiBSJjo6GxWJBXV0dfvrTn7Z5jo+PD8LDwzttTB8fH0yYMAEVFRV2x8+cOYPBgwd32jiiYYAmots0NDSgsrJS/nz+/HmUl5ejf//+uP/++5GSkoK5c+dizZo1iI6OxuXLl1FUVIQxY8Zg+vTpnTreoEGDAACLFy/G7NmzMWnSJEyePBl79uzBRx99hH379rl8v8KSiIh+5PPPP5cA3NbmzZsnSZIktbS0SDk5OdKQIUOk3r17SyEhIdIvf/lL6fjx46qM1+qtt96S7rvvPsnPz0+KjIyUdu3a5eKdio3roImIBMVldkREgmKAJiISFB8SEpEwmpqa0NLS4lRfHx8fj9uXgwGaiITQ1NSEsMF3wVRncaq/Xq/H+fPnPSpIM0ATkRBaWlpgqrPgfNlgaAOUVV/rr1kRNu4CWlpaGKCJiNTS966bTQmLh65F40NCIiJBMYMmIqFYIcEKZSmx0vO7CwZoIhKKFVZYnejjiRigiUgoFkmCReEXnJWe310wQBORUFjisGGAJiKhWCHBwgANgKs4iIiExQyaiITCEocNM2giEkrrQ0KlTQmj0YgJEyYgICAAQUFBSE5Ovu1tLW157733EB4eDj8/P4wePRqffPKJs7fpEAZoIhKK1cmmRHFxMdLT03Hw4EEUFhbixo0bmDp1KhobG9vtc+DAATz22GNYsGABjh07huTkZCQnJ+PkyZOK79FR3LCfiIRQX18PnU6Hr04FIUDhXhzXrlkxMqIOZrMZWq1W8diXL19GUFAQiouLMWnSpDbPmT17NhobG7F792752IMPPoioqCjk5eUpHtMRzKCJSCgWybkG3Azyt7bm5maHxjSbzQCA/v37t3tOaWkpEhIS7I4lJiaitLTUuRt1AAM0EXkMg8EAnU4nN6PReMc+VqsVmZmZmDhxIkaNGtXueSaTCcHBwXbHgoODYTKZXJ53e7iKg4iE4kxNufX86upquxKHr6/vHfump6fj5MmT2L9/v8JR1ccATURCsUIDCzSK+wCAVqtVVIPOyMjA7t27UVJSgnvvvbfDc/V6PWpra+2O1dbWQq/XK5qrEixxEJFQrJJzTQlJkpCRkYGdO3di7969CAsLu2Of2NhYFBUV2R0rLCxEbGysssEVYAZNREKxOJFBKz0/PT0dBQUF+OCDDxAQECDXkXU6Hfz9/QEAc+fOxT333CPXsZ999lk8/PDDWLNmDaZPn47t27fjyJEj2LRpk6KxlWAGTURCaQ3QSpsSGzZsgNlsRnx8PEJCQuS2Y8cO+ZyqqipcunRJ/hwXF4eCggJs2rQJkZGR+N///V/s2rWrwweLruI6aCISQus66ANfheAuheugG65ZETfyktProEXFEgcRCcUqaWCVFD4kVHh+d8EATURC6YoadHfBAE1EQrHACxaFj8csKs3F3RigiUgokhMlDoklDiIi9bHEYcNldkREgmIGTURCsUhesEgKa9Aeuli4ywO01WrFxYsXERAQAI3GM38tIeqpJEnCtWvXEBoaCi8v535Bt0IDq8Jf7j31lVddHqAvXrwIg8HQ1cMSUReqrq6+4+ZD7WEN2qbLA3RAQAAA4N7fvggvP7+uHp7a0Lda3UcRkY98per1AeDA4QjVx+gKQ5ccVn2Mr1dNUO3a1qYmfPvbV+T/nzvDuRIHM+hO0VrW8PLzY4AWhLevugHa5y4fVa8PwGP+W+ql6a36GF3x78qV8uXNEodz2416Gq7iICISFFdxEJFQrE58k5APCYmIugBr0DYM0EQkFCu8uMzuBwzQRCQUi6SBReHeGkrP7y6ceki4fv16DBkyBH5+foiJicGhQ4c6e15ERD2e4gC9Y8cOZGVlYdmyZTh69CgiIyORmJiIuro6NeZHRD1M63ajSpsnUnxXa9euRVpaGlJTUzFixAjk5eWhT58+2Lx5sxrzI6Iexip5OdU8kaIadEtLC8rKypCdnS0f8/LyQkJCAkpLS9vs09zcjObmZvlzfX29k1Mlop7AuQ37PfMhoaJ/C1euXIHFYkFwcLDd8eDgYPm15T9mNBqh0+nkxn04iKgjVtgeFDrarO6etEpU/70gOzsbZrNZbtXV1WoPSUTdWOsyO6XNEykqcQwcOBDe3t6ora21O15bWwu9Xt9mH19fX/j6+jo/QyKiHkrRjx0fHx+MGzcORUVF8jGr1YqioiLExsZ2+uSIqOdp/Sah0uaJFN9VVlYW3nzzTWzduhWnTp3C008/jcbGRqSmpqoxPyLqYVp3s1PalCopKcGMGTMQGhoKjUaDXbt2dXj+vn37oNFobmvtPX/rDIq/STh79mxcvnwZOTk5MJlMiIqKwp49e257cEhE5Azn9uJQnkE3NjYiMjISjz/+OB555BGH+1VUVECr1cqfg4KCFI/tKKe+6p2RkYGMjIzOngsRkZPL7JQH6GnTpmHatGmK+wUFBaFfv36K+znDMws3RNRtWSWNUw24+T2LW9ut38HoLFFRUQgJCcG//Mu/4K9//WunX/9WDNBE5DEMBoPd9y6MRmOnXTskJAR5eXl4//338f7778NgMCA+Ph5Hjx7ttDF+jLvZEZFQnNuw/+b51dXVdvXhzlziO3z4cAwfPlz+HBcXh3PnzmHdunX405/+1Gnj3IoBmoiE4szeGq3na7VauwCttgceeAD79+9X7foM0EQkFAs0sChcNqf0/M5SXl6OkJAQ1a7PAO0CH7P6JfyBX1pUH+PiJHV3MigpHanq9T1J5boHVR8jtES9jYW+vyGhysVruJJBK9HQ0IDKykr58/nz51FeXo7+/ftj0KBByM7ORk1NDd5++20AQG5uLsLCwjBy5Eg0NTXhj3/8I/bu3Yu//OUvisd2FAM0EQnFAuUZsTNpzJEjRzB58mT5c1ZWFgBg3rx5yM/Px6VLl1BVZftx09LSgkWLFqGmpgZ9+vTBmDFj8Nlnn9ldo7MxQBNRjxQfHw+pg5fN5ufn231+7rnn8Nxzz6k8K3sM0EQklK4qcXQHDNBEJJSu+qp3d8AATURCkZzY/Ehy0yoOtTFAE5FQmEHbMEATkVBu3VtDSR9P5Jk/doiIPAAzaCI
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = np.random.rand(100).reshape(10,10)\n",
    "P, L, U = scipy.linalg.lu(A)\n",
    "plot_matrices([P, L, U], titles=[\"P\", \"L\", \"U\"])\n",
    "plot_matrices([P@L@U - A])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d0e2812",
   "metadata": {},
   "source": [
    "## Eigen-decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19159d23",
   "metadata": {},
   "source": [
    "The eigenvectors of an $n\\times n$ matrix $\\mathbf{A}$ satisfy\n",
    "\n",
    "$$\\mathbf{A}\\mathbf{v_i} = \\lambda_i \\mathbf{v_i}$$\n",
    "\n",
    "with eigenvalues $\\lambda_i$. Now consider the matrix $Q$ whose columns correspond to the $n$ eigenvectors $\\mathbf{v_i}$: you can easily show that this satisfies\n",
    "\n",
    "$$\\mathbf{A}\\mathbf{Q} = \\mathbf{Q}\\mathbf{\\Lambda},$$\n",
    "\n",
    "where $\\mathbf{\\Lambda} = \\mathrm{diag}(\\lambda_1, \\lambda_2\\dots \\lambda_n)$.\n",
    "\n",
    "Therefore \n",
    "\n",
    "$$\\mathbf{A} = \\mathbf{Q}\\mathbf{\\Lambda}\\mathbf{Q^{-1}}$$\n",
    "\n",
    "which is the *eigendecomposition* of the matrix.\n",
    "\n",
    "Since $(\\mathbf{A}\\mathbf{B}\\mathbf{C})^{-1} =\\mathbf{C}^{-1}\\mathbf{B}^{-1}\\mathbf{A}^{-1}$, the inverse of the matrix is\n",
    "\n",
    "$$\\mathbf{A^{-1}} = \\mathbf{Q}\\mathbf{\\Lambda}^{-1}\\mathbf{Q^{-1}}$$\n",
    "\n",
    "where \n",
    "\n",
    "$$\\mathbf{\\Lambda}^{-1} = \\mathrm{diag}(1/\\lambda_1, 1/\\lambda_2\\dots 1/\\lambda_n).$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "43a8a74f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues= [11.57481875  0.11814652  8.81628573  1.38184428  2.236823    3.37700902\n",
      "  7.11855695  6.48539964  5.20294079  4.33568771]\n",
      "Eigenvector check:\n",
      "0 1.0436096431476471e-14\n",
      "1 8.259872278336411e-15\n",
      "2 5.440092820663267e-15\n",
      "3 2.345346139520643e-15\n",
      "4 4.6629367034256575e-15\n",
      "5 3.3306690738754696e-15\n",
      "6 2.6645352591003757e-15\n",
      "7 5.218048215738236e-15\n",
      "8 5.329070518200751e-15\n",
      "9 1.7763568394002505e-15\n"
     ]
    }
   ],
   "source": [
    "# real symmetric matrix\n",
    "n = 10\n",
    "A = np.random.rand(n*n).reshape(n,n)\n",
    "A = 0.5*(A + A.T)\n",
    "A = A + np.diag(np.arange(n))\n",
    "lam, Q = scipy.linalg.eig(A)\n",
    "lam = np.real(lam)\n",
    "Q = np.real(Q)\n",
    "\n",
    "print(\"Eigenvalues=\", lam)\n",
    "\n",
    "print(\"Eigenvector check:\")\n",
    "for i in range(n):\n",
    "    # The eigenvectors are given by the columns of Q, ie. Q[:,i]\n",
    "    print(i, np.max(np.abs(A@Q[:,i]-lam[i]*Q[:,i])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "e3adf1fc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVsAAAFXCAYAAAAMDe2WAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAkwElEQVR4nO3de3BUVbo28KfTkE6QdCaEXIcmJHjhDoFASjJj4ZDiUsA3WBZqVZgJcU5ETYQYS010AEckTUbkpAQmXEqRKrn5lRP14IDDCQLDDBSQCELJRUSgJ5gEHExDkA7pvc8fmNY2JGR1eu1e6Ty/qlVltnuvtTYkL2/evfZqk67rOoiISKqQQE+AiKg7YLAlIjIAgy0RkQEYbImIDMBgS0RkAAZbIiIDMNgSERmAwZaIyAAMtkREBmCwJSIyAIMtEbWyd+9ezJgxA4mJiTCZTPjggw8CPt6cOXNgMpm82pQpU6TOy58YbImolcbGRowcORKrVq1SarwpU6bgm2++8bTNmzcbMj9/6BHoCRCReqZOnYqpU6e2+f9dLhdefvllbN68Gd999x2GDRuG0tJSTJgwQcp4LSwWC+Lj430aI9CY2RKRsPz8fOzfvx9btmzB559/jlmzZmHKlCn48ssvpY67e/duxMbG4r777sNTTz2Fb7/9Vup4/mTiFotE1B6TyYSKigrMnDkTAHDhwgWkpKTgwoULSExM9JyXmZmJcePGoaSkxK/jtdiyZQt69eqF5ORkfPXVV3jppZfQu3dv7N+/H2azuVNjGoFlBCIScuzYMbjdbtx7771ex10uF6KjowEAJ0+exODBg9vt58UXX8TSpUs7PO5jjz3m+e/hw4djxIgRGDhwIHbv3o2JEycK3EFgMNgSkZBr167BbDajqqqqVUbZu3dvAEBKSgpOnDjRbj8tgdlXKSkp6Nu3L86cOcNgS0TBJzU1FW63G/X19fj1r39923NCQ0MxaNAgqfP497//jW+//RYJCQlSx/EXBlsiauXatWs4c+aM5+uvv/4aR44cQZ8+fXDvvfciKysLv//97/HGG28gNTUVly5dQmVlJUaMGIFp06b5dbz+/fvj2rVr+NOf/oSHH34Y8fHx+Oqrr/DCCy/g7rvvxuTJk/1yz9LpREQ/8+mnn+oAWrXs7Gxd13W9qalJX7hwoT5gwAC9Z8+eekJCgv7QQw/pn3/+uZTxrl+/rk+aNEmPiYnRe/bsqSclJem5ubl6bW2tn+5YPq5GICIyANfZEhEZgMGWiLodt9uNBQsWIDk5GeHh4Rg4cCAWL14Mmb/o8wEZESnjxo0baGpq8una0NBQhIWFdejc0tJSlJeXY8OGDRg6dCgOHz6MnJwcREZGYt68eT6Nfyes2RKREm7cuIHkpN6orXf7dH18fDy+/vrrDgXc6dOnIy4uDm+99Zbn2MMPP4zw8HC8++67Po1/J8xsiUgJTU1NqK134+uqJFgjxCqczqsaksecx+XLl2G1Wj3HLRYLLBZLq/PHjx+PtWvX4vTp07j33ntx9OhR7Nu3D8uXL+/0fbSFwZaIlHJX71tNhPuH389tNpvX8UWLFuGVV15pdX5RURGcTicGDRoEs9kMt9uNJUuWICsry8dZ3xmDLREFDYfD0SqzvZ333nsPGzduxKZNmzB06FAcOXIEBQUFSExMRHZ2tpS5MdgSkVI06NAg9iip5Xyr1eoVbNvy/PPPo6ioyLO5zfDhw3H+/HnY7XYGWyLqHjRo0Hy4RsT169cREuJdFzabzdA00ZE7jsGWiJTi1nW4BRdJiZ4/Y8YMLFmyBP3798fQoUPx2WefYfny5Xj88ceF+hHBYEtESulMGaGjVqxYgQULFuDpp59GfX09EhMTMXfuXCxcuFCoHxFcZ0tESnA6nYiMjMTXJxMQIbj06+pVDcmDvkFDQ0OHaraBwNd1iYgMwDICESnFiDJCIDDYEpFSjHhAFggMtkSkFO2HJnqN6hhsiUgpbuhwC5YFRM8PBAZbIlKKW/9xrwORa1TH1QhERAZgZktESmHNlojIABpMcMMkfI3qGGyJSCmafquJXqM6BlsiUorbh8xW9PxAYLAlIqUEa7DlagQiIgMwsyUipWi6CZou+IBM8PxAYLAlIqUEaxmBwZaIlOJGCNyCFU63pLn4E4MtESlF96GMoLOMQEQkJljLCFyNQERkAGa2RKQUtx4Cty5Ys+UbZK1pmoaLFy8iIiICJpP6qT8RdZyu67h69SoSExMREuLbL84aTNAEf+nmx+LcxsWLF2Gz2YwelogM5HA40K9fP5+uDdaareHBNiIiAgDwy/8uQki4Rdo4MZ/2lNZ3i8Z4+SXv3hflbx7XHCb3G7X3v5uk9g8Auln+D5srUv6PS9Shi9LHaL5QI69v3MQ+/M3zc+4L38oIzGxbaSkdhIRbEBIeJm0cc6j8YGu2yA+25p7yg60eKjdQ9egh/8/JiGDr7in/x6VHiLwExMMk8Wfjh5jXmRLhrTJC8G2xyNUIREQG4GoEIlKK5sMbZHxARkQkiDVbIiIDaAjh0i8iItncugluwb0ORM8PBJ8ekK1atQoDBgxAWFgY0tPTcfDgQX/Pi4goqAgH261bt6KwsBCLFi1CdXU1Ro4cicmTJ6O+vl7G/Iiom2nZYlG0iaqpqcHs2bMRHR2N8PBwDB8+HIcPH5ZwR7cIz3D58uXIzc1FTk4OhgwZgtWrV6NXr154++23ZcyPiLoZTQ/xqYm4cuUKMjIy0LNnT2zfvh1ffPEF3njjDURFRUm6K8GabVNTE6qqqlBcXOw5FhISgszMTOzfv/+217hcLrhcLs/XTqfTx6kSUXfg2+bhYg/ISktLYbPZsH79es+x5ORkoT5ECd3R5cuX4Xa7ERcX53U8Li4OtbW1t73GbrcjMjLS07gvAhG1R8OPD8k62lres3Q6nV7tp4neT3300UdIS0vDrFmzEBsbi9TUVKxbt07qfUl/g6y4uBgNDQ2e5nA4ZA9JRF1Yy9Iv0QYANpvNK7mz2+23HePs2bMoLy/HPffcg08++QRPPfUU5s2bhw0bNki7L6EyQt++fWE2m1FXV+d1vK6uDvHx8be9xmKxwGIx4H1vIur2HA4HrFar5+u2Yo+maUhLS0NJSQkAIDU1FcePH8fq1auRnZ0tZW5CmW1oaCjGjBmDyspKzzFN01BZWYn777/f75Mjou6n5Q0y0QYAVqvVq7UVbBMSEjBkyBCvY4MHD8aFCxek3ZfwSw2FhYXIzs5GWloaxo0bh7KyMjQ2NiInJ0fG/IiomzFi16+MjAycOnXK69jp06eRlJQk1I8I4WD76KOP4tKlS1i4cCFqa2sxatQo7Nixo9VDMyIiX/i2N4LY+c8++yzGjx+PkpISPPLIIzh48CDWrl2LtWvXCvUjwqfXdfPz85Gfn+/vuRAR+bj0S+z8sWPHoqKiAsXFxXj11VeRnJyMsrIyZGVlCfUjgnsjEJFSNN0ETXCvA9HzAWD69OmYPn268HW+4ubhREQGYGZLRErxbfNw9fNGBlsiUoovex2Inh8IDLZEpBR+lLmfRf8jFObQUGn916fL/1TaXr9skD5GncP3j4TuqLB6ud+o38fIf4PQFS1/p37dLH8M5wD5e4ccn/8/0vp2XtUQdW/n+mBmS0RkADfEM1W3nKn4lfr/HBARBQFmtkSkFJYRiIgMYMTruoHAYEtEStF92IhG52oEIiIxzGyJiAxg1N4IRlP/nwMioiDAzJaIlGLEFouBwGBLREoJ1jICgy0RKeWnn5Yrco3qGGyJSClu3QS3YKYqen4gMNgSkVKCtYygfu5NRBQEmNkSkVJ0H/ZG0PlSAxGRGG4eTkRkAE0Xr8Fq8vd17zQGWyJSCrdYJCIygObDrl+i5weC+v8cEBEFAWa2RKQUvtRARGQA1myJiAygwYc3yLpAzZbBloiUwo/F8bPGeBPMFnl/QGEJjdL6bpGeeF76GFWfjpA+hi75u+BqilvuAAAsl83Sx3CluKSPoV0Jkz7Gb+b8l7S+m5tvAHilU31wbwQioiC1dOlSmEwmFBQUSBuDZQQiUorRD8gOHTqENWvWYMQIub9
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAFXCAYAAACPydStAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAnPklEQVR4nO3df1RUZf4H8PcMPwZEQFHkxzYK2g/8laKYgbulK0dpLbU8lrt0VHK1NSgN1w3cxNKU7JiyqfnrmLonSetslOt3szXyx5ooipHaCoaVEDaglYzgyuDc+/3DuNMo4NwZLvMwvF/nPOfE9T73ea7Yhw+f+9xndLIsyyAiIuHo3T0BIiJqGgM0EZGgGKCJiATFAE1EJCgGaCIiQTFAExEJigGaiEhQDNBERIJigCYiEhQDNBGRoBigiegWBw8exCOPPILIyEjodDp88MEHbh/vpZdeQkxMDAICAtC1a1ckJibi6NGjms7L3RigiegWdXV1GDRoENauXSvMeHfffTfWrFmDU6dO4dChQ4iKisKYMWNw8eLFNpmjO+i4WRIRtUSn0yEvLw8TJ05UjtXX1+Ovf/0r3nnnHVy+fBkDBgzA8uXLMXLkSE3Ga4rZbEZwcDA++eQTjB492uVxRcQMmohUS0tLQ0FBAXbs2IGTJ09i8uTJSEpKwldffdUm41ssFmzcuBHBwcEYNGhQm4zpDt7ungARtS/l5eXYsmULysvLERkZCQD485//jD179mDLli1YtmyZZmPv3r0bU6ZMwdWrVxEREYG9e/eie/fumo3nbsygiUiVU6dOwWq14u6770bnzp2VduDAAZw7dw4AUFJSAp1O12LLyMhQPfaoUaNQXFyMw4cPIykpCY8//jiqq6tb+xaFwQyaiFSpra2Fl5cXioqK4OXlZfdnnTt3BgD07t0bZ86cafE63bp1Uz12QEAA7rzzTtx55524//77cdddd2Hz5s3IzMxUfa32gAGaiFSJjY2F1WpFdXU1fvOb3zR5jq+vL2JiYjSfiyRJqK+v13wcd2GAJqJb1NbWoqysTPn6m2++QXFxMUJCQnD33XcjOTkZU6dOxeuvv47Y2FhcvHgR+fn5uPfeezFu3LhWHa9nz56oq6vD0qVLMX78eERERODSpUtYu3YtKisrMXny5Fa5ZyHJREQ32bdvnwzgljZt2jRZlmXZYrHIWVlZclRUlOzj4yNHRETIjz76qHzy5ElNxvvf//4nP/roo3JkZKTs6+srR0REyOPHj5cLCwtb6Y7FxHXQRESC4ioOIiJBMUATEQmKDwmJSBjXrl2DxWJxqq+vry/8/PxaeUbuxQBNREK4du0aont1hqna6lT/8PBwfPPNNx4VpBmgiUgIFosFpmorvinqhaBAddVX8xUJ0UPPw2KxMEATEWkloPONpobVQ9ei8SEhEZGgmEETkVAkyJCgLiVWe357wQBNREKRIEFyoo8nYoAmIqFYZRlWlS84qz2/vWCAJiKhsMRhwwBNREKRIMPKAA2AqziIiITFDJqIhMIShw0DNBEJhQ8JbRigiUgo0s9NbR9PxABNREKxOvGQUO357QUDNBEJxSqr31uDe3EQEVGbYgZNREJhDdqGAZqIhCJBByt0qvt4IgZoIhKKJN9oavt4IgZoIhKK1YkMWu357QUDNBEJhQHahqs4iIgExQyaiIQiyTpIssqHhCrPby8YoIlIKCxx2DBAE5FQrNDDqrL6atVoLu7GAE1EQpGdKHHILHEQEWmPJQ4bruIgIhIUM2giEopV1sMqq6xB803C1iFJEi5cuIDAwEDodJ75awlRRyXLMq5cuYLIyEjo9c79gi5BB0nlL/f8yKtWcuHCBRiNxrYelojaUEVFBe644w6n+rIGbdPmATowMBAA8Gv8Dt7w0WycvLOnNLt2oxer7tV8jIr/ddV8jEv/C9D0+vd3/1bT6wNAWV2o5mMY/X/SfIz/1oRrPkaw4Zpm126os+D/JuYq/587w7kSBzPoVtFY1vCGD7x12gXooEDtn38a6rSbfyMfva/mY3jrDZpe39C5Df6eoP3fk6GT9vfhfV3b7wUA+Bi03z3ZlfLljRIHtxsFuIqDiEhYXMVBREKRnHiTkA8JiYjaAGvQNgzQRCQUCXous/sZa9BEJBSrrHOqOWPt2rWIioqCn58fhg8fjsLCwmbP3bp1K3Q6nV3z8/Nz9jYd4lSAVnNTRERqNO5mp7aptXPnTqSnp2PRokU4ceIEBg0ahLFjx6K6urrZPkFBQfj++++Vdv78eVdu9bZU35UzN0VEJJqVK1di5syZSElJQb9+/bB+/Xp06tQJb731VrN9dDodwsPDlRYWFqbpHFUHaGduiojIUZKsd6qpYbFYUFRUhMTEROWYXq9HYmIiCgoKmu1XW1uLXr16wWg0YsKECfjyyy+dvk9HqLorZ26qvr4eZrPZrhERNceVEsfNsaa+vr7JMS5dugSr1XpLBhwWFgaTydRkn3vuuQdvvfUWPvzwQ7z99tuQJAkJCQn47rvvWvcv4BdUBWhnbio7OxvBwcFK4z4cRNQSCeofFDa+G2k0Gu3iTXZ2dqvNKz4+HlOnTsXgwYPx4IMP4v3330doaCg2bNjQamPcTPNldpmZmUhPT1e+NpvNDNJE1CznltndOL+iogJBQUHKcYOh6Vfnu3fvDi8vL1RVVdkdr6qqQni4Y/uh+Pj4IDY2FmVlZarmqoaqvwVnbspgMCAoKMiuERFp4eZY01yA9vX1xdChQ5Gfn68ckyQJ+fn5iI+Pd2gsq9WKU6dOISIiolXm3hRVAbo1boqIqCWNbxKqbWqlp6dj06ZN2LZtG86cOYPZs2ejrq4OKSkpAICpU6ciMzNTOX/x4sX497//ja+//honTpzAk08+ifPnz+OPf/xjq937zVSXONLT0zFt2jTExcXhvvvuQ05Ojt1NERG5oq12s3viiSdw8eJFZGVlwWQyYfDgwdizZ4/yjK28vNzuQwd++uknzJw5EyaTCV27dsXQoUNx+PBh9OvXT/XYjlIdoG93U0RErnBuLw7nXopOS0tDWlpak3+2f/9+u69XrVqFVatWOTWOs5x6SNjSTRERucKZNwOdeZOwPeBmSUQkFEnWQVK5t4ba89sLz/yxQ0TkAZhBE5FQnNuw3zNzTQZoIhKKM3trqD2/vWCAJiKhWKGDVeWyObXntxduC9B5Z09p+snb9/xnqmbXbhQZUqP5GAE+Fs3HCPWv1fT6h6r7aHp9AOjs2/SmOK3p+A89NR+jq+Gq5mOYLdptMn/d4nqgZAZtwwyaiIRihfqM2KrNVNzOM3/sEBF5AGbQRCQUljhsGKCJSCht+aq36BigiUgoshObJclcxUFEpD1m0DYM0EQkFO7FYeOZP3aIiDwAM2giEgq3G7VhgCYiobDEYcMATURCceVTvT0NAzQRCcUq62BVmRGrPb+9YIAmIqGwxGHjmb8XEBF5AGbQRCQU2Ym9OGS+qEJEpD1u2G/DAE1EQpFk9TVlSdZoMm7GAE1EQuF2ozYM0EQkFMmJ3ezUnt9eeOaPHSIiD8AMmoiEwhdVbBigiUgorEHbMEATkVAkOPEmoYfWoBmgiUgo/MgrG7cF6KzqgTBc9dHs+pEhNZpdu5G/d4PmY+h12i/wtFi1/WfQycei6fWBtvl7CmiD+2iLX9W99ZJm15b1rn8fuBeHjWcWboiIPABLHEQkFD4ktPHMuyKidquxxKG2OWPt2rWIioqCn58fhg8fjsLCwhbPf++99xATEwM/Pz8MHDgQ//rXv5wa11EM0EQklMY3CdU2tXbu3In09HQsWrQIJ06cwKBBgzB27FhUV1c3ef7hw4fx+9//HjNmzMDnn3+OiRMnYuLEiTh9+rSrt9wsBmgiEkpbZdArV67EzJkzkZKSgn79+mH9+vXo1KkT3nrrrSbP/9vf/oakpCTMnz8fffv2xZIlSzBkyBCsWbPG1VtuFgM0EQmlLQK0xWJBUVEREhMTlWN6vR6JiYkoKChosk9BQYHd+QAwduzYZs9vDXxISEQew2w2231tMBhgMBhuOe/SpUuwWq0ICwuzOx4WFoaSkpImr20
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Check eigendecomposition\n",
    "plot_matrices([Q@np.diag(lam)@np.linalg.inv(Q) - A])\n",
    "\n",
    "# and the inverse\n",
    "plot_matrices([Q@np.diag(1/lam)@np.linalg.inv(Q) - np.linalg.inv(A)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "3bb677c2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAFXCAYAAACPydStAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAly0lEQVR4nO3df1RUdf4/8OfMKAMY4E8YyFHQMvz9CyV0t+wjJ3TNk/vt6+ouHYkK+wGV0ZbQJm6akm0am3ok3WO45/irPZvl6Qd9jFKPK4VClJY/0lxhzQEtmRGUQefezx/GnWYF4w68Z94Mz8c573NivPe+36Px5MX7vu97DKqqqiAiIukY/T0AIiJqGQOaiEhSDGgiIkkxoImIJMWAJiKSFAOaiEhSDGgiIkkxoImIJMWAJiKSFAOaiEhSDGgius7evXsxc+ZMxMTEwGAw4J133vF7fw888AAMBoNHmzZtmtBx+RsDmoiu09DQgNGjR2Pt2rVS9Tdt2jScPXtWa1u3bvXJ+Pylm78HQETymT59OqZPn97qnzudTvzpT3/C1q1bUVdXhxEjRmDFihWYMmWKkP6amc1mWCwWr/rojFhBE5FuWVlZKC0txbZt2/DVV19h9uzZmDZtGr799luh/e7evRuRkZG47bbb8Nhjj+GHH34Q2p+/GbjdKBHdiMFgwI4dOzBr1iwAQFVVFQYNGoSqqirExMRoxyUnJ2PixIlYvnx5h/bXbNu2bQgNDUVcXBxOnjyJ559/HjfddBNKS0thMpna1aesOMVBRLocOnQILpcLQ4YM8Xjd6XSiT58+AICjR49i6NChN7zOwoUL8fLLL7e537lz52r/PXLkSIwaNQqDBw/G7t27MXXqVB3voPNgQBORLvX19TCZTCgvL7+ucr3pppsAAIMGDcKRI0dueJ3mMPfWoEGD0LdvX5w4cYIBTUQEAGPHjoXL5UJtbS1+/etft3hMUFAQ4uPjhY7jP//5D3744QdER0cL7cefGNBEdJ36+nqcOHFC+/rUqVOorKxE7969MWTIEKSmpmLevHlYuXIlxo4di3PnzqGkpASjRo3CjBkzOrS/AQMGoL6+Hi+++CLuu+8+WCwWnDx5Es899xxuueUWpKSkdMh7lpJKRPRfPv30UxXAdS0tLU1VVVVtampS8/Ly1NjYWLV79+5qdHS0+tvf/lb96quvhPR36dIl9e6771b79eundu/eXR04cKCakZGh2my2DnrHcuIqDiIiSXEdNBGRpBjQRESS4k1CIpJGY2MjmpqavDo3KCgIwcHBHTwi/2JAE5EUGhsbETfwJthqXV6db7FYcOrUqYAKaQY0EUmhqakJtloXTpUPRHiYvtlXx0UFceNPo6mpiQFNRCRKj5uuNT1cAboWjTcJiYgkxQqaiKSiQIUCfSWx3uM7CwY0EUlFgQLFi3MCEQOaiKTiUlW4dD7grPf4zoIBTURS4RSHGwOaiKSiQIWLAQ2AqziIiKTFCpqIpMIpDjcGNBFJhTcJ3RjQRCQV5aem95xAxIAmIqm4vLhJqPf4zoIBTURScan699bgXhxERORTrKCJSCqcg3ZjQBORVBQY4IJB9zmBiAFNRFJR1GtN7zmBiAFNRFJxeVFB6z2+s2BAE5FUGNBuXMVBRF3S3r17MXPmTMTExMBgMOCdd975xXN2796NcePGwWw245ZbbkFRUZHQMTKgiUgqimrwqunV0NCA0aNHY+3atW06/tSpU5gxYwbuuusuVFZWYsGCBXj44Yfx0Ucf6e67rTjFQURS8dUUx/Tp0zF9+vQ2H19YWIi4uDisXLkSADB06FDs27cPr732GlJSUnT33xasoIlIKi4YvWoA4HA4PJrT6eywcZWWliI5OdnjtZSUFJSWlnZYH/+NAU1EUlG9mN5Qf5risFqtiIiI0Fp+fn6HjctmsyEqKsrjtaioKDgcDly+fLnD+vk5TnEQkVTaM8VRXV2N8PBw7XWz2dyhY/M1BjQRBYzw8HCPgO5IFosFNTU1Hq/V1NQgPDwcISEhQvpkQBORVFyqES5V3+yrL3azS0pKwgcffODx2q5du5CUlCSsT58HtKIo+P777xEWFgaDITAXlxN1Vaqq4uLFi4iJiYHR6N0tLgUGKDpvj3nzkVf19fU4ceKE9vWpU6dQWVmJ3r17Y8CAAcjNzcWZM2fw97//HQDw6KOPYs2aNXjuuefw4IMP4pNPPsFbb72F999/X3ffbeXzgP7+++9htVp93S0R+VB1dTX69+/v1bm+WmZ38OBB3HXXXdrX2dnZAIC0tDQUFRXh7NmzqKqq0v48Li4O77//Pp5++mn89a9/Rf/+/fG3v/1N2BI7ADCoqm8/zMtut6Nnz54YmLMIRnOwsH7UAWLuqv5c6IFQ4X04RjUJ70M0U11gzKS5Qn2wqWWQ+D6M3V3Crq1cduI/T76Curo6RERE6DrX4XAgIiICO768FT3CTLrObbjowm9Hfwu73S5sDtoffP6d0zytYTQHwxgsMKBDxf/cMQn8AdPMGNL5V0IaGwMjoNWQAAnoIHEB3aw905fXpji43SjAddBERNIKjNKGiAKG8rMnA9t+TmBuCM2AJiKpeLfMjgFNRCScAqNPltl1BgxoIpKKSzXApXP7UL3HdxZe3SRcu3YtYmNjERwcjMTERJSVlXX0uIiIujzdAb19+3ZkZ2dj8eLFqKiowOjRo5GSkoLa2loR4yOiLqY9240GGt3vatWqVcjIyEB6ejqGDRuGwsJChIaGYuPGjSLGR0RdjKIavWqBSNccdFNTE8rLy5Gbm6u9ZjQakZyc3Oqm1U6n02PTbIfD4eVQiagr8KYidgXoTUJdfwvnz5+Hy+VqcdNqm83W4jn5+fkeG2hzHw4iuhEF7huFbW0+eMbTL4T/XpCbmwu73a616upq0V0SUSfWvMxObwtEuqY4+vbtC5PJ1OKm1RaLpcVzzGZzp/9UAyIif9D1YycoKAjjx49HSUmJ9pqiKCgpKRG6aTURdR3NTxLqbYFI94Mq2dnZSEtLQ0JCAiZOnIiCggI0NDQgPT1dxPiIqIvhbnZuugN6zpw5OHfuHPLy8mCz2TBmzBgUFxdfd+OQiMgb3u3FwQpak5WVhaysrI4eCxGRl8vsGNBERMIpqgGKzr019B7fWQTmjx0iogDACpqIpOLdhv2BWWsyoIlIKt7srcG9OIiIfMAFA1w6l83pPb6z8FtAK2YVMIvb4ET5QfzTi/aRV4T3EXw6SHgfqknsRjPdLon/5rkc6YNPw74svkrrcVz8t2RjP3H/3obG9l+DFbQbK2gikooL+itil5ih+F1g/tghIgoArKCJSCqc4nBjQBORVPiotxsDmoikonqxWZLKVRxEROKxgnZjQBORVLgXh1tg/tghIgoArKCJSCrcbtSNAU1EUuEUhxsDmoik4s2ndHM3OyIiH3CpBrh0VsR6j+8sGNBEJBVOcbgF5u8FREQBgBU0EUlF9WIvDpUPqhARiccN+90Y0EQkFUXVP6esiP3MCb9hQBORVLjdqBsDmoikonixm53e4zuLwPyxQ0QUAFhBE5FU+KCKGwOaiKTCOWg3BjQRSUWBF08SBugcNAOaiKTCj7xy81tAm88ZYTKL+7XE1Cjs0ppuyReE91EX0kN4H+GlIUKvfzW5Tuj1AQCnw4V30f9jRXgf/RaeEN5HZdktwq6tdMCCZO7F4RaYEzdERAGAUxxEJBXeJHRjQBORVDjF4caAJiKp8ElCNwY0EUmFFbQbA5qIpMKAdgvMmXUiojZYu3YtYmNjERwcjMTERJSVlbV6bFFREQwGg0cLDg4WOj4GNBFJpbmC1tv02r59O7Kzs7F48WJUVFRg9OjRSElJQW1tbavnhIeH4+zZs1o7ffp0e97qL9IV0Pn5+ZgwYQLCwsIQGRmJWbNm4dixY6LGRkRdkK8CetWqVcjIyEB6ejqGDRuGwsJChIaGYuPGja2eYzAYYLFYtBYVFdWet/qLdAX0nj17kJmZic8++wy7du3ClStXcPfdd6OhoUHU+Iioi1HhXsnR1qb3+cWmpiaUl5cjOTlZe81oNCI5ORmlpaWtnldfX4+
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_matrices([np.linalg.inv(Q) - Q.T])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f725b10f",
   "metadata": {},
   "source": [
    "## Singular value decomposition (SVD)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c8d2cc0",
   "metadata": {},
   "source": [
    "The singular value decomposition (SVD) of a $m\\times n$ matrix $\\mathbf{A}$ is\n",
    "\n",
    "$$\\mathbf{A} = \\mathbf{U}\\ \\mathbf{S}\\ \\mathbf{V^T},$$\n",
    "\n",
    "where \n",
    "\n",
    "- $\\mathbf{U}$ is an $m\\times n$ matrix with orthonormal columns\n",
    "- $\\mathbf{S}$ is an $n\\times n$ diagonal matrix whose diagonal elements are the **singular values** of the matrix $s_i$\n",
    "- $\\mathbf{V}$ is an $n\\times n$ orthogonal matrix ($V^TV = VV^T = 1$)\n",
    "\n",
    "If we have a square matrix ($m=n$), we can write the inverse as\n",
    "\n",
    "$$\\mathbf{A^{-1}} = \\mathbf{V}\\,\\mathbf{diag}(1/s_i)\\,\\mathbf{U^T}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "66806267",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Singular values are:  [1.76883682 0.58748002 0.34937691]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/sAAAFlCAYAAAC0i/ggAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAByuUlEQVR4nO3dfVhUZcIG8HsGBXR1QFK+FEXTFFNBIdhBt3QlQb3YaF0zc8NYxa11WpW2glJR2aTyCzOSSpHcZDUr7UNflDB0TdJESW2VFr8gY1CXBMEEZeb9g2Vy5AxwzjBzhpn7917n2ubM88x5hre99zzzfByFXq/Xg4iIiIiIiIjshlLuBhARERERERFR+2Jnn4iIiIiIiMjOsLNPREREREREZGfY2SciIiIiIiKyM+zsExEREREREdkZdvaJiIiIiIiI7Aw7+0RERERERER2hp19IiIiIiIiIjvDzj4RERERERGRnWFnn4iIiIiIiMjOsLNPslm0aBEUCgUuXLgg+H6fPn0wZswY6zaKiIiIiIjIDrCzT7L59ttv4e7uDn9//2bv/fe//8WlS5cQFBRk9XYRERERERF1dOzsk2y+/fZbjBgxwuR7ABAYGGjNJhEREREREdkFdvZJFteuXUNpaanJzjw7+0RERERERNKxs0+yOHHiBADTnflvv/0WSqUSw4YNs2aziIiIiIiI7AI7+ySL1kbuv/32WwwaNAhdu3a1ZrOIiIiIiIjsAjv7JItvv/0WTk5OgiP3t27dwr///W9uzkdERERERCQRO/skixMnTuC+++6Dq6trs/fOnDmD+vp6rtcnIiIiIiKSiJ19ksW///1vBAQECL536NAhAEBYWJg1m0RERERERGQ32Nknq7t9+zZqa2uhUCgE39+6dSt69OiB0aNHW7llRERERERE9oGdfbK6Tp06wc/PDwcPHkR1dbXRe5s2bUJ+fj40Gg1cXFxkaiEREREREVHHptDr9Xq5G0GOJz09HRqNBoMGDcL06dPRuXNnfPXVV8jJycFDDz2EnJwcwfX8ZB9u3ryJ+vp6SXWdnZ357wYR2S3mIxGRMHPyEXDMjGRnn2Tz8ccfY8WKFfj3v/+NW7duYciQIZg+fTrmz5+Pzp07y908spCbN2+if79u0F5ukFTf29sb58+fd7iwJiL7x3wkIhJmbj4CjpmR7OwTkVVVV1fDzc0NFwv9oeoubiVR9XUd+gVfQFVVFVQqlYVaSEQkD+YjEZEwc/IRcNyM7CR3A4jIMXXrrkC37sKbNJqig7jyREQdEfORiEiYlHwEHDcj2dknIlk06HVoEDmvqEGvs0xjiIhsCPORiEiYlHxsqueI2NknIlnooIcO4tJabHkioo6I+UhEJExKPjbVc0QWe/ReZWUlZsyYAZVKBXd3d8yaNQs1NTUt1hk7diwUCoXR8fTTT1uqiUQkI53E/yMisnfMRyIiYVLz0VEz0mIj+zNmzEB5eTlyc3Nx69YtxMXFYc6cOcjOzm6xXnx8PJYtW2Z43bVrV0s1kYhk1KDXo0Hk/qBiyxMRdUTMRyIiYVLysameI7JIZ//06dPIycnBN998g5CQEADAunXrMGnSJKxcuRK+vr4m63bt2hXe3t6WaBYRERERERGRQ7BIZ7+goADu7u6Gjj4AREREQKlU4vDhw3j00UdN1t2yZQvef/99eHt7Izo6GosWLWpxdL+urg51dXWG1zqdDpWVlbjnnnugUDjmrotE1qDX63H9+nX4+vpCqRS/IohrUomIhDEfiYiEcc2+OBbp7Gu1Wnh6ehpfqFMneHh4QKvVmqz3xBNPoF+/fvD19cWJEyfw4osvori4GB9//LHJOqmpqVi6dGm7tZ2IxCkrK0OfPn1E19NBjwbezBIRNcN8JCISJiUfm+o5IlGd/cTERLz22mstljl9+rTkxsyZM8fwz8OHD4ePjw/Gjx+Ps2fP4t577xWsk5SUhISEBMPrqqoq9O3bF/c/sQhOzq6S22LvroXflLsJNu/bsVvkboJNq67Rod+oC+jevbuk+tYcuUpPT8eKFSug1WoRGBiIdevWITQ01GT5tLQ0rF+/HqWlpejZsyf+8Ic/IDU1Fa6uzBQisjyO7BMRCePIvjiiOvvPPfccnnrqqRbLDBgwAN7e3rh8+bLR+du3b6OyslLUevywsDAAQElJicnOvouLC1xcXJqdd3J2ZWe/BUrue9gqVXeLPazCrkhdLmOtDai2bduGhIQEZGRkICwsDGlpaYiMjERxcXGzGUgAkJ2djcTERGRmZiI8PBzff/89nnrqKSgUCqxevVr09YmIxOIGfUREwrhBnziiOvu9evVCr169Wi2nVqtx7do1FBYWIjg4GACwb98+6HQ6Qwe+LYqKigAAPj4+YppJRB2A7n+H2DpirV69GvHx8YiLiwMAZGRkYNeuXcjMzERiYmKz8ocOHcLo0aPxxBNPAAD8/f0xffp0HD58WMLViYjEs1Y+EhF1NFLysameI7LI0GVAQACioqIQHx+PI0eO4KuvvoJGo8Hjjz9u2In/0qVLGDJkCI4cOQIAOHv2LFJSUlBYWIgLFy7g008/RWxsLB588EGMGDHCEs0kog6qurra6Lhzk8471dfXo7CwEBEREYZzSqUSERERKCgoEKwTHh6OwsJCQzadO3cOu3fvxqRJk9r/ixARERERWYhFNugDGnfV12g0GD9+PJRKJaZMmYI33njD8P6tW7dQXFyMGzduAACcnZ3xxRdfIC0tDbW1tfDz88OUKVOwcOFCSzWRiGTUIGGDlabyfn5+RueTk5OxZMmSZuWvXr2KhoYGeHl5GZ338vLCmTNnBK/xxBNP4OrVqxgzZgz0ej1u376Np59+Gi+99JKothIRSWVOPhIR2TMp+dhUzxFZrLPv4eGB7Oxsk+/7+/tDf8faCT8/P+zfv99SzSEiG9OgbzzE1gEanwCgUqkM54X27ZAqPz8fy5cvx1tvvYWwsDCUlJRg3rx5SElJwaJFi9rtOkREppiTj0RE9kxKPjbVc0QW6+wTEbXEnDWpKpXKqLNvSs+ePeHk5ISKigqj8xUVFSY3C120aBGefPJJzJ49G0Djk0Fqa2sxZ84cvPzyy1AquXEjEVkW1+wTEQnjmn1xeNdKRLLQQYEGkYcO4nb+d3Z2RnBwMPLy8n65rk6HvLw8qNVqwTo3btxo1qF3cnICAKPZSERElmKNfCQi6oik5KMjZyRH9olIFjp94yG2jlgJCQmYOXMmQkJCEBoaatgXpGl3/tjYWPTu3RupqakAgOjoaKxevRojR440TONftGgRoqOjDZ1+IiJLslY+EhF1NFLysameI2Jnn4js2rRp03DlyhUsXrwYWq0WQUFByMnJMWzaV1paajSSv3DhQigUCixcuBCXLl1Cr169EB0djVdeeUWur0BEREREJBo7+0Qki6ZpVWLrSKHRaKDRaATfy8/PN3rdqVMnJCcnIzk5WdK1iIjMZc18JLIXGzduxLx58wA0PnpXr9cbNvANDw/H3r175WwetRMp+dhUzxGxs09EsuDNLBGRMOYjkXizZs3CrFmzAABz5sxBt27dsHr1aplbRe2NnX1xuEEfEclCp1dIOoiI7B3zkcg8J06cwIgRI+RuBlmA1Hx01IxkZ5+IZCFlJ1VH/VWWiBwL85GoUVxcHFxdXdHQ0GCyzMSJE9G1a1f88MMPABqfnHPq1Cl29u2U1Hx01IxkZ5+IZNEApaSDyJGcOHECM2bMQO/eveHs7Axvb2+MGTMGS5culbtpZEHMR6JGAQEBqKurw/nz5wXf/9e//oWcnBz89a9/RZ8+fQAAZ8+eRV1dHe6//35rNpWsRGo+OmpGOua3JiIisnEff/wxHnjgARw5cgR//vOf8dZbb2HOnDnQ6XR4//335W4eEZHFDR06FABw5swZwfeTkpLg4eGBxMREw7lvv/0W9913n2FzPiKp0tPT4e/vD1dXV4SFheHIkSMmy44dOxYKhaLZMXnyZEOZp556qtn7UVFRFv0O7OwTkSz0EtZa6R10vRU5np9++gl/+tOf8MADD+C7777D4sWLMXv2bCxbtgyHDh3
      "text/plain": [
       "<Figure size 1200x400 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArAAAAFWCAYAAACRsa60AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0mElEQVR4nO3df3RV1Z3//9cNkhsykABiEn6khlb5JZAgSAy0CqsZUmUxTdvpouhqMFX8akkHTWeEWCQqnzFtVZp2RDP+QGZmyYLqFOgSioOhgUFiKWBQrOCASBiGG2AoCcSahNz7/YNy7Z38MOcmJ3efzfOx1lmr99x9ztlJ4eWbffbZxxcKhUICAAAAPCIu1h0AAAAAnKCABQAAgKdQwAIAAMBTKGABAADgKRSwAAAA8BQKWAAAAHgKBSwAAAA8hQIWAAAAnkIBCwAAAE+hgAUAAICnUMACuCLt2LFDc+bM0bBhw+Tz+bRhw4bPPaaqqko33nij/H6/rrvuOq1evdr1fgKAyWKVpRSwAK5IjY2NyszM1MqVK7vU/ujRo5o9e7ZmzpypmpoaPfDAA7rnnnv0xhtvuNxTADBXrLLUFwqFQtF0GABs4fP5tH79euXn53fYZvHixdq0aZMOHDgQ3ved73xH586d05YtW3qhlwBgtt7MUkZgAaALqqurlZubG7EvLy9P1dXVMeoRAHhPT2XpVT3ZKQDork8//VTNzc1RHRsKheTz+SL2+f1++f3+bvcrEAgoNTU1Yl9qaqoaGhr0pz/9Sf369ev2NQCgJ5iao1LPZSkFLABjfPrppxp5bX8FTrVGdXz//v114cKFiH2lpaV69NFHe6B3AGC+KyVHKWABGKO5uVmBU606uvdaJQ1wNsOp4XxQIycf0/Hjx5WUlBTe31OjBmlpaaqrq4vYV1dXp6SkJEZfARjD5ByVei5LKWABGOev+l/anGj98+OoSUlJEcHbU3JycrR58+aIfVu3blVOTk6PXwsAusvEHJV6Lkt5iAuAcYIKRbU5ceHCBdXU1KimpkbSpaVdampqVFtbK0kqKSlRQUFBuP19992njz76SA899JAOHjyoZ599Vr/85S/14IMP9tjPDQA9pTdyVIpdljICC8A4QQUVjOIYJ/bs2aOZM2eGPxcXF0uS5s+fr9WrV+vkyZPhAJakkSNHatOmTXrwwQf185//XCNGjNCLL76ovLw8hz0FAPf1Ro5KsctS1oEFYIyGhgYlJyfrfw6NiGru1rDR/636+nrXbn0BgOmulBxlBBaAcVpDIbU6/Le10/YAYDPbc5QCFoBxopmLFc3cLQCwle05SgELwDhBhdRqcfACgNtsz1EKWADGsX3kAADcZnuOUsACMI7tc7cAwG225yjrwAIAAMBTGIEFYJzgnzenxwAALrE9RylgARinNYqHD5y2BwCb2Z6jFLAAjNMa+uyd3E6OAQBcYnuOMgcWxnjppZfUv39/9e/fX/Hx8erbt2/486xZs2LdPfSiYJQbcCUjQ/GXbM9RXiULI917773q37+/VqxYEeuuoBddfgXivj+kqr/DVyBeOB/UjePqPPEKRMBtZOiV60rJUUZgYaR3331XEydOjHU3AMCTyFDYjgIWrnrkkUfk8/n08ccft/v9iBEj9OUvfzliXygU0oEDBzoM32jOCW8JhqLbANsUFhYqISFBra2tHba57bbblJiYqP/+7/+W1HmGRnM+eJPtOUoBC1ft379fAwcOVEZGRpvv/vd//1cnTpxQVlZWxP4jR46oqalJN9xwQ4+dE97SKl9UG2CbsWPHqqmpSUePHm33+//8z//Uli1b9Hd/93caMWKEpM4zNJrzwZtsz1EKWLhq//79HY6k7t+/X5KUmZnZZv+oUaPk9/t77JzwFtuDF+iqcePGSZIOHjzY7vclJSUaPHiwlixZEt7XWYZGcz54k+05SgEL15w7d061tbUdFpMdFZudzd2K9pzwlmDIF9UG2OZywfnBBx+0+W7Tpk166623VFJSooEDB4b3d5ah0ZwP3mR7jrIOLFzz7rvvSuq4mNy/f7/i4uI0fvz4NsdNnTq1R88Jb4lmJMBLIwdAV2VkZKhfv35tRkxDoZCWLl2q9PR0FRUVRXzXWYZGcz54k+05yggsXPN5o6H79+/X9ddfr8TExIj9nY0eRHtOAPCiuLg4jRkzps2I6bp161RTU6PHH39cCQkJEd91lqHRnA8wESOwcM3+/fvVp0+fdkdDW1pa9Ic//EHf+MY32nx35MiRHj8nvKVVcWp1+O/rjp+pBrxt3Lhx+s1vfhP+fPHiRS1btkwTJkxQQUFBm/adZWg054M32Z6jjMDCNe+++65GjRrV7r/mDx48qObmZsdzVd04J8wTimLeVshDc7cAJ8aOHauzZ8/q1KlTkqSXX35Z//Vf/6WysjLFxTn/z3hPnw9msj1H+ZMK1/zhD3/Q2LFj2/1u165dkqTs7OyYnxPmsf3pWcCJv3zwqqmpScuXL9ctt9yi2bNnG3E+mMn2HGUKAVxx8eJFNTY2yudr/y/D2rVrNWjQIE2fPj2m54SZWkNxag05vPXloQW4ASf+cumrffv26fjx43r11VeNOR/MZHuOUsDCFVdddZXS09O1c+dONTQ0RLxT+eWXX1ZVVZUeeeSRDtd67a1zwkxB+RR0eIMoKA8lL+DAddddp/j4eP3+97/Xr3/9a33zm9/s1p2mnj4fzGR7jlLAwjWLFy9WUVGRpkyZonnz5qlv37566623tGXLFt166616+OGHjTgnAJisT58+GjVqlFavXi2fz6cnnnjCqPMBscAcWLhm4cKF+vd//3ddffXVKi8v1xNPPKG6ujr99Kc/1datW6NaqsWNc8I8ts/dApwaN26cWltb9b3vfU+jR4827nwwj+056guFQt4ZLwZgtYaGBiUnJ2v9/uv1VwP6ODq28XyrvpH5X6qvr4+YXgIAV5IrJUeZQgDAOJfmbjkbCXDaHgBsZnuOUsACME4wigW4vfTwAQC4zfYcpYAFYJzoln/xTvACgNtsz1Ee4gIAAICnuFbAnj17VnfeeaeSkpI0cOBA3X333bpw4UKnx8yYMUM+ny9iu++++9zqIgBDBRUX1WYjshRANGzPUdemENx55506efKktm7dqpaWFhUWFuree+/VmjVrOj1uwYIFevzxx8OfExMT3eoiAEO1hnxqdfhObqftvYIsBRAN23PUlQL2gw8+0JYtW/T73/9eU6ZMkST90z/9k26//XY99dRTGjZsWIfHJiYmKi0tzY1uAfCI1igePmj10MMHXUWWAoiW7TnqylhxdXW1Bg4cGA5cScrNzVVcXJx+97vfdXrsK6+8oiFDhmj8+PEqKSnRJ5984kYXARgsGIqLarMNWQogWrbnqCsjsIFAQCkpKZEXuuoqDR48WIFAoMPj7rjjDl177bUaNmyY3n33XS1evFiHDh3Sr371qw6PaWpqUlNTU/hzMBjU2bNndfXVV8vn885QOGCTUCik8+fPa9iwYYqLcx6Ito8cdFVvZSk5CpiHHO2cowJ2yZIl+slPftJpmw8++CDqztx7773h/z1hwgQNHTpUX/3qV3XkyBF96UtfaveYsrIyPfbYY1FfE4B7jh8/rhEjRsS6G8YxLUvJUcBc5Gj7HBWwP/zhD3XXXXd12uaLX/yi0tLSdOrUqYj9Fy9e1NmzZx3NycrOzpYkHT58uMMCtqSkRMXFxeHP9fX1+sIXvqBj+zKU1N87Q+Gx8o1RE2LdBVjoolq0U5s1YMCAqI4PyvnDBMGorhQbpmUpOdo95CjcQI52zlEBe8011+iaa6753HY5OTk6d+6c9u7dq8mTJ0uStm3bpmAwGA7SrqipqZEkDR06tMM2fr9ffr+/zf6k/nFKGkDwfp6rfH1j3QXY6M93oaK9/RzNci5eWv7FtCwlR7uHHIUryNFOudLTsWPH6mtf+5oWLFig3bt366233lJRUZG+853vhJ+aPXHihMaMGaPdu3dLko4cOaLly5dr7969+vjjj/XrX/9aBQUFuuWWWzRx4kQ3ugnAUJffION0sw1ZCiBatueoa+vAvvLKKyo
      "text/plain": [
       "<Figure size 800x400 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArAAAAFWCAYAAACRsa60AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABINElEQVR4nO3dfVRU17k/8O+AMkgR1BpeVAymvkEUUAhk8CZqQ8SXa0Ob+rM2LUiVrKRMq5k0ieQaidI6aX3DVitRQ0jbsLB6q8mNBkM0aI34BpJoouRqVKhxQK+RURIHnJnfH5SJEwbkzMyZObP9ftbaa3XO7DPnGRofHvbZZ2+V1Wq1goiIiIjIR/h5OwAiIiIiIilYwBIRERGRT2EBS0REREQ+hQUsEREREfkUFrBERERE5FNYwBIRERGRT2EBS0REREQ+hQUsEREREfkUFrBERERE5FNYwBIRERGRT2EBS0RERERO2b9/P2bOnIlBgwZBpVJhx44ddzynsrIS48ePh1qtxvDhw1FSUiL5uixgiYiIiMgpLS0tiI+Px/r163vU/9y5c5gxYwYmT56M2tpaLFy4EPPnz8fu3bslXVdltVqtzgRMRERERNRBpVJh+/btyMjI6LLPCy+8gJ07d+LkyZO2Yz/5yU9w7do1lJeX9/haHIElIiIiIo+oqqpCWlqa3bH09HRUVVVJ+pxe7gyKiMhVN2/eRGtrq1PnBgQEIDAw0M0RERH5FlfyqNVqhUqlsjumVquhVqvdERoMBgPCw8PtjoWHh8NoNOLrr79Gnz59evQ5LGCJSDFu3ryJYfcGw9Bkdur8iIgInDt3jkUsEd21XM2jwcHBuHHjht2x/Px8vPzyy26Izn1YwBKRYrS2tsLQZMa56nsR0lfaDCfjdQuGJV5Aa2srC1giumu5I482NDQgJCTEdtxdo69A+0BDY2Oj3bHGxkaEhIT0ePQVYAFLRAr0neD2JoWZj6MSEdm4kkdDQkLsClh30mg02LVrl92xiooKaDQaSZ/Dh7iISHEssDrViIionafy6I0bN1BbW4va2loA7ctk1dbWor6+HgCQl5eHzMxMW/+nnnoKn3/+OZ5//nmcPn0af/7zn/H3v/8dzzzzjKTrcgSWiBTHAgssTpxDRETtPJVHjx07hsmTJ9te63Q6AEBWVhZKSkpw6dIlWzELAMOGDcPOnTvxzDPPYO3atRgyZAg2b96M9PR0SddlAUtERERETpk0aRK621LA0S5bkyZNwvHjx126LgtYIlIcs9UKs8Q9VqT2JyISmeh5lAUsESmOM3OxOAeWiOgboudRFrBEpDgWWGEWOPESEclN9DzKApaIFEf0kQMiIrmJnkdZwBKR4og+d4uISG6i51GuA0tEREREPoUjsESkOJZ/N6nnEBFRO9HzKAtYIlIcsxMPH0jtT0QkMtHzKAtYIlIcs/WbPbmlnENERO1Ez6OcA0s+6bXXXkNwcDCCg4MREBCA3r17215PmTLF2+GRiyxONiLqGeZQ8YmeR1XW7vb/IvIBTz75JIKDg7F69Wpvh0IuMhqNCA0NRc2n4QjuK+3v6xvXLRgf24jm5maEhITIFCGReJhDxXK35FGOwJLP+/jjjxEXF+ftMIiIfBJzKPkiFrAkm5deegkqlQrnz593+P6QIUPwH//xH073BwCr1YqTJ08y+QrGYnWuEYkmOzsbgYGBMJvNXfaZNm0agoKCJPX917/+BYA5VGSi51EWsCSbjz76CP369UN0dHSn9/7v//4PFy9eREJCgtP9AeDs2bMwmUy4//773Rs8eZUZKqcakWhiYmJgMplw7tw5h+//85//RHl5OX79619L6jtkyBAAzKEiEz2PsoAl2Xz00Udd/lX/0UcfAQDi4+Od7t9xfOTIkVCr1e4ImRRC9MRL1FOxsbEAgNOnTzt8Py8vDwMGDMCiRYsk9e3AHCou0fMoC1iSxbVr11BfX9+p4Ozw7YJUav8OnLslJotV5VQjEk1HUXrq1KlO7+3cuRMffvgh8vLy0K9fP0l9OzCHikv0PMp1YEkWH3/8MYDOBWeHjz76CH5+fhgzZoxT/W+/TnJysrvCJoVwZiTAl0YOiHoqOjoaffr06TSqarVasXjxYkRFRUGr1Uru24E5VFyi51GOwJIsuhoxvf39ESNGICgoyKn+HTh6QEQi8/Pzw+jRozuNqm7ZsgW1tbVYtmwZAgMDJfftwBxKvoojsCSLjz76CP7+/p1GTAGgra0Nn376KX74wx863b/D2bNn3Rs4KYIZfjBL/Pu66+euiXxbbGws3n33XdvrW7duYcmSJRg7diwyMzOd7gswh4pM9DzKEViSxccff4yRI0d2+msfaH/AoLW11W60VWp/EpvViXlbVh+au0UkRUxMDK5evYqmpiYAwOuvv47//d//hV6vh5+fn9N9SWyi51H+10yy+PTTTxETE+PwvYMHDwIAUlJSnO5PYhP96VkiKW5/OMtkMqGgoAAPP/wwZsyY4VJfEpvoeZRTCMjtbt26hZaWFqhUjv8hlJWVoX///pgwYYJT/Ul8ZqsfzFaJt758aAFuIiluXx6rpqYGDQ0N2Lp1q8t9SWyi51EWsOR2vXr1QlRUFA4cOACj0Wi3n/Lrr7+OyspKvPTSS7Z1B6X2J/FZoIJF4g0iC3wo8xJJMHz4cAQEBODo0aN4++238aMf/ajLO1JS+pLYRM+jLGBJFi+88AK0Wi2SkpIwZ84c9O7dGx9++CHKy8sxceJEvPjiiy71JyK6W/j7+2PkyJEoKSmBSqXC8uXL3dKXyJdxDizJIjc3F//93/+N7373uygsLMTy5cvR2NiIP/zhD6ioqOj0sJbU/iQ20eduEUkVGxsLs9mMX/ziFxg1apTb+pK4RM+jKqvV6jvjxUQkNKPRiNDQUGz/aAS+09df0rkt1834Yfz/orm52W4aChHR3eRuyaOcQkBEitM+d0vaSIDU/kREIhM9j7KAJSLFsTixALcvPXxARCQ30fMoC1giUhznln/xncRLRCQ30fMoH+IiIiIiIp8iWwF79epVPPHEEwgJCUG/fv0wb9483Lhxo9tzJk2aBJVKZdeeeuopuUIkIoWywM+pJiLmUiJyhuh5VLYpBE888QQuXbqEiooKtLW1ITs7G08++SRKS0u7PS8nJwfLli2zvQ4KCpIrRCJSKLNVBbPEPbml9vcVzKVE5AzR86gsBeypU6dQXl6Oo0ePIikpCQDwpz/9CdOnT8fKlSsxaNCgLs8NCgpCRESEHGERkY8wO/HwgdmHHj7oKeZSInKW6HlUlrHiqqoq9OvXz5ZwASAtLQ1+fn44fPhwt+e++eabGDhwIMaMGYO8vDx89dVXcoRIRApmsfo51UTDXEpEzhI9j8oyAmswGBAWFmZ/oV69MGDAABgMhi7P++lPf4p7770XgwYNwscff4wXXngBdXV1+Mc//tHlOSaTCSaTyfbaYrHg6tWr+O53vwuVyneGwolEYrVacf36dQwaNAh+ftITougjBz3lqVzKPEqkPMyj3ZNUwC5atAi///3vu+1z6tQpp4N58sknbf977NixiIyMxCOPPIKzZ8/ie9/7nsNz9Ho9li5d6vQ1iUg+DQ0NGDJkiLfDUByl5VLmUSLlYh51TFIB++yzz2Lu3Lnd9rnvvvsQERGBpqYmu+O3bt3C1atXJc3JSklJAQCcOXOmywI2Ly8POp3O9rq5uRlDhw7FGwdGIChY2hZqvmLdk497OwTZGVK/4+0QZBVZ2P3tX193C204gF3o27evU+dbIP1hAotTV/IOpeXSrvLohZpohAT7zi1Fb/nhyLHeDoEExDzaPUkF7D333IN77rnnjv00Gg2uXbuG6upqJCYmAgD27t0Li8ViS6Q9UVtbCwCIjIzsso9arYZare50PCjYH0ES9wD2Fb16BXo7BNn5q8X+jr1Uvb0dgrz+fRfK2dvPzizn4kvLvygtl3aVR0OC/RDS13d+rt4i/L9n8g7m0W7JEmlMTAymTp2KnJwcHDlyBB9++CG0Wi1+8pOf2J6avXjxIkaPHo0jR44AAM6ePYuCggJUV1fj/PnzePvtt5G
      "text/plain": [
       "<Figure size 800x400 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAAF4CAYAAACB05i5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAsOklEQVR4nO3dfVRU17k/8O+AziDCYIjCoIwvKb6ACiKKQXMVGxSN1yVtlzHe3AUhvvQFcmNJjZKV6FJ/6SQ1vt1o1TQa0heq11RMq6mWoEiNGCuKilVaco0QLwOm6hCnOuDM/v1hmHRkQM4whzln/H7W2mv1nNnnnGem+rjznH320QghBIiIyKcCfB0AERExGRMRKQKTMRGRAjAZExEpAJMxEZECMBkTESkAkzERkQIwGRMRKQCTMSnWjh07EBISgpCQEGi1WvTs2dO5PX36dF+HR+RVGj6BR2qwePFihISEYP369b4OhUgWHBmTKpw7dw7x8fG+DoNINkzG5NZrr70GjUaDzz//3O3n0dHReOKJJ1z2nTt3Ds8++ywGDBgArVYLg8GAJ554AqtWrQIAZGdnIygoCHa7vd3rzpw5E8HBwfjiiy+c+4QQqKqq8kkyvnjxInr06AGNRoPr1693+/Xp4cFkTG6dPXsWffr0weDBg9t89o9//ANXr17FmDFjnPv27t2L8ePH4+TJk/j+97+Pn//851i8eDEcDgd+/etfAwBiY2Nhs9lw+fJlt9f885//jIMHD+K//uu/EB0d7dz/2WefwWazYeTIkV79jp2Rn58Po9EIAKiqqur269PDo4evAyBlOnv2bLsj0bNnzwIAEhISAAA3btzA888/j/Hjx+Pw4cPQarXOvqtXr8b//d//AQDi4uIAAJcuXUJMTEyb8+bn5yM8PBzLly9vc71hw4ZBp9N1/YtJcPz4cRw4cADFxcWYOnUqqqqqMHny5G6NgR4eHBlTGzdv3kRtba0z2d7v/mT88ccfw2KxYMGCBS6JuFX//v0BfJOML1682KbPgQMH8MknnyA/Px99+vRx+cxX9eJly5bhueeeQ2pqKsLCwjgyJllxZExtnDt3DgA6TMYBAQEYNWoUAMBqtboc157BgwejV69euHTpkst+IQReffVVGI1G5Obmuo0nOTlZ8vfoit///vc4deoUfvvb3wK4V2I5f/58t8ZADxeOjKmN+0e+7j4fOnQogoODAQBTp05FcHAwNm7ciGHDhmHZsmU4cuRImxt1AQEBGDFiRJuR8e7du1FZWYnVq1cjKCiozfW6e2Rst9vxyiuv4Ec/+pGzdh0bG4sLFy50Wwz0EBJE91mwYIEIDAwUt2/fbvNZc3Oz0Gq1Yt68eS77z549K55++mnRu3dvAUAAEAMHDhTFxcUu/Z599lkRHh7u3G5paRFDhw4Vo0ePFna7vUtx2+12cfv27U41h8PR7nl27NghQkNDxbVr15z7fvaznwkA4osvvuhSjETt4ciY2jh37hyGDRvmdpR66dIlNDc3txk1x8fHY/fu3bh+/TqOHDmC7Oxs1NXV4ZlnnsE///lPZ7/Y2Fhcv34djY2NAID33nsPf//732EymRAQ0LU/jmVlZejVq1enWnV1tdtz3LlzBytXrsQzzzyDmzdvoqamBjU1NdDr9QDAUgXJhjVjauOvf/0r0tPT3X52/PhxAMCECRPcfq7VapGamorU1FRcu3YN+/fvR3V1NRITEwG43sQLCwvDmjVrMHnyZMyaNavLcY8YMQLvvfdep/pGRUW53b9p0yZ88cUX+MUvfoFf/OIXbT6vqqrCjBkzuhSnPygrK8PatWtRUVGB+vp6FBUVISMjw+fXu3jxIpYtW4ajR4/i7t27iIuLw+9+9zsMHDhQtti8hcmYXNy9exdWqxUajcbt57t27cIjjzyCSZMmPfBcrVPRwsLCnPv+dXrb6dOnUVdXhz179nghcsBgMOC5557z+PgbN27gjTfewOLFizFt2rQ2n2dlZXFGxdesVisSEhLw/PPP47vf/a4irvfZZ5/hiSeewIIFC7Bq1Sro9XpcuHDB7X/hKZKv6ySkPEajUURGRgqLxeKyf+fOnQKAeO2115z7/vznP4t//vOfbc5x9uxZ0bt3b5GYmOiy/+7du0Kr1YoFCxaIfv36ie9+97vyfAkP/OQnPxGhoaHixo0bbj8fNWqUGDt2bPcGpQIARFFRkcu+O3fuiJdeekn0799fBAcHi+TkZHHkyBHZrieEEPPmzRP/+Z//6ZVr+AKTMbWxefNmAUAMHTpUrFixQqxZs0bMmDFDABBTpkxxubE3adIk0a9fP/GjH/1IbN++XWzZskUsWrRIBAUFicjISFFVVdXm/KNGjRKBgYGiR48e4tKlS9351dpVW1srgoKCxE9+8pN2+8yZM0f06tWryzca/Y275Lhw4UIxceJEUVZWJmpqasTatWuFTqcTf/vb32S5nt1uFyEhIWL16tVi+vTpol+/fiI5Odlt0lYqJmNy63e/+514/PHHhV6vF7169RKJiYniZz/7mWhubnbpt3fvXjF//nwRExMjevfuLYKCgkRsbKxYunSpaGxsdHvup59+WgAQixcv7o6v0inPPfec0Gq14urVq+32ycvLEwC8klD8yf3J8cqVKyIwMLDNb/nkk0+K/Px8r19PCCHq6+sFABEcHCzWr18vzpw5I0wmk9BoNKK0tLTL1+wOTMZE1CX3J8f9+/cLAKJ3794urUePHuLpp58WQghx8eJF5xTI9tqyZcs6dT0hhLh69aoAIObPn++yf/bs2eKZZ57x6veVC2/gEZFX3bp1C4GBgaioqEBgYKDLZyEhIQCAxx57zO1j8f/q0Ucf7fQ1+/btix49ejhvELeKjY3FsWPHOn0eX2IyJiKvSkxMhN1uR2NjI/7t3/7NbR+tVosRI0Z47ZparRbjx49vM3/8b3/7GwYNGuS168iJyZiIJLt16xZqamqc25cvX0ZlZSXCw8MxbNgwPPvss8jMzMS6deuQmJiIa9euoaSkBPHx8R7NKe/oeq1ziJcuXYp58+Zh8uTJmDp1Kg4ePIg//OEPKC0t7fL37Ra+rpMQkfocOXLEbZ03KytLCHHvsfkVK1aIwYMHi549e4qoqCjxne98R5w7d06W67XasWOHiImJEUFBQSIhIUHs27evi9+0+8j2Drzr16/jhRdewB/+8AcEBATge9/7HjZt2uSsGbmTmpqKo0ePuuz7/ve/j23btskRIhGRYsiWjGfOnIn6+nps374dLS0tyM7Oxvjx41FYWNjuMampqRg2bBhWr17t3BccHOxcF4CIyF/JUjO+ePEiDh48iL/85S8YN24cAODtt9/GU089hbfeesu52Lg7wcHBMBgMcoRFRKRYsiTj8vJy9OnTx5mIASAtLQ0BAQH49NNP8Z3vfKfdY3/zm9/g17/+NQwGA2bPno3XXnvNuW6uOzabDTabzbntcDhw/fp1PProo+2ur0BE3U8Iga+++gr9+/f3eIW+O3fuoLm52aNjtVqtotepkCUZm81mREREuF6oRw+Eh4fDbDa3e9x//Md/YNCgQejfvz/OnTuHZcuWobq6Gnv37m33GJPJ5Hz7MBEpX11dncsLZzvrzp07GDIoBObG9t8u3hGDwYDLly8rNiFLSsbLly/Hm2++2WGfB03k7sjixYud/3v06NGIiorCk08+ic8++wzf+ta33B6Tn5+PvLw857bFYsHAgQMx/PkVCNQq80dXok/y3vV1CKoz48fP+ToEVbnbcgcVf/opQkNDPTq+ubkZ5kY7LlcMgj5U2si66SsHhiRdQXNzs38k45deeumBSxQ+9thjMBgMzsXDW929exfXr1+XVA9uXTO3pqam3WSs0+ncvjU4UBuEQJ0yf3QlkvqHm4AePfnnyxNdLR/qQwP88s+rpGTcr18/9OvX74H9UlJScPPmTVRUVCApKQkAcPjwYTgcjnYXJXensrISQPsLgRPRw8cuHLBLnANmFw55gvEiWf55iY2NxYwZM7Bo0SKcPHkSn3zyCXJzc/HMM884Z1JcvXoVI0aMwMmTJwHcWxh6zZo1qKiowOeff47
      "text/plain": [
       "<Figure size 400x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m=4; n=3\n",
    "\n",
    "A = np.random.rand(m*n).reshape(m,n)\n",
    "\n",
    "U, Sdiag, VT = np.linalg.svd(A, full_matrices=False)\n",
    "print(\"Singular values are: \", Sdiag)\n",
    "\n",
    "S = np.diag(Sdiag)\n",
    "plot_matrices([U, S, VT], titles=[r\"$U$\", r\"$S$\", r\"$V^T$\"])\n",
    "\n",
    "# U and V should be orthogonal\n",
    "plot_matrices([U.transpose()@U,VT@VT.transpose()], titles=[r\"$U^TU$\",r\"$V^TV$\"])\n",
    "plot_matrices([U@U.transpose(),VT.transpose()@VT], titles=[r\"$UU^T$\",r\"$VV^T$\"])\n",
    "\n",
    "# Check the reconstructed matrix\n",
    "plot_matrices([U@S@VT - A], titles=[r\"$USV^T - A$\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4a9ab62",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}