3M1 Examples Paper 1 Crib.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3M1 Examples Paper crib (linear algebra)\n",
    "\n",
    "\n",
    "Garth N. Wells (gnw20@cam.ac.uk)\n",
    "\n",
    "(C) 2016-2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import linalg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 1\n",
    "\n",
    "The usual rules for real matrices extend without modification to complex valued matrices.\n",
    "\n",
    "a. $\\det(\\boldsymbol{A}) = (4 + 4i)(4 + i) - (2-i)(-3 + 2i) = 16 + 13 i$\n",
    "\n",
    "b. $\\det(\\boldsymbol{A}^{H}) = (4 - 4i)(4 - i) - (-3 - 2i)(2 + i) = 16 - 13i$\n",
    "\n",
    "c. Inverse:\n",
    "\n",
    "   $$\n",
    "   \\boldsymbol{A}^{-1} = \\frac{1}{16 + 13 i } \n",
    "   \\begin{bmatrix}\n",
    "      4 + i   & -2 + i  \\\\\n",
    "      3 - 2i  & 4 + 4i  \n",
    "   \\end{bmatrix}\n",
    "   $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Extra:** [SymPy](https://www.sympy.org/) can be used for symbolic computations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJEAAAAzCAYAAACXDOrGAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAG2ElEQVR4Ae2d7XHVOBiFTWYLCCkhdMBmOwgdbHY7gA7Y4Vfyb4ftAFoIHcBWkIUOoARICXsejXTH9pV9JVmSdbnWjJGvrI/zHh29kj8Untze3p53nnB3d/foSd6STpgBacKrlTNx8lnHj9Hx9oS52kyfZuB+pBN0c/9EnuirTl5JZZ+my25XNgb2GZBmXir1Bk+0hY2BRQz8sqi0p7DU+VzJeLZXnsvVkw7h0fXXAvVsbbxq/1I4/rIEXSn+zm+lf7Fpq0Qh/JTwRMybF6tY7G/0EB7EvupNhBXQO8Vm8Cn+VZgQz2edX/vNqpZ6kJ+snkgGM6oXBUvapeL3iypS4RA8yvNsaTsZynMjM/DcwoUXYs3BIHiaoY2kKkL4yeaJ1BjTGCN66ajmNpJjUciIZxGOwMJ4m6/CPLabm51zpTPVNRuyiUgW/iljF3uPjEy1hmfONMTyTfxNDcCxuObqqn4ty3Qm45nG3lVHP9HgITy6TqcwVTCVse5YVfxq/2bCFLx7p+tVF9dqL4qfxZ5IDeJqHxV/myCianIgnjfK94+Asd5o8sGq8CEguHV3bDqtFqL4WSwimcUdxaojeUTtLB7bOQ+2DB6gCfGPbOAnAv9gxe65XCYphZ9F05kaZEpImsZUlnIsKMfBPB7Q9cHdis30RelTrr/TtRA8rD3c9PCH6v17DGDt38IHN+CctNVhVB6mnn91xKybbnocuKpcHM1PsogEAlfLnUPSSFY5n0g6pf+uernFZ7oJDsofhEf5HqlUMe1AfEteFFwMhAvFLxQfDMqHPTxXyhJsfZ3iYH6SRSTEdNpvagy32w9mLrfpqLrWnB6L541AM12wnmMwGHH1Dal9Lgx0HE/Pdx5I59jVKU4arAtsCOYnWUQyitvSvZe2SufN7ifFOyIWGBJcNAaP8uKBELsTOIS58+A2c2YUJvAwKMc4EFZVbxnLT7KIZgikgzhaCT48bnQjds75kmG1YDHg0cHDeqgfrpUWNbX3CyeeR/GTTUTWeNO4gGM4pDysQIDhbQ6PrrFAZypj/dEprjrSDcDhPx/1E+4MnuEl8w5tlFT2p/iI4ieniLwL5bLmTtcuImbx6HrV6XYaqRFxC+/vBhBj+DkblGzjx6NgcGzhSBjI5oly2asRsLdYz1X3Vk8ZBlr0RGUs3WotxsAmomLUnk7Fm4hOp6+LWbqJqBi1p1PxJqLT6etilm4iKkbt6VS8ieh0+rqYpZuIilF7OhUXf9ioh4eXotO9guBlKL/frvVQ0eJxb8qvhOW7jslNgsq/yuZGtctb/dU3gYbYX1REAoBo6CAnok7nfNrwUTFf133QebWg9hAwmwR3H3zpnG+s+Vj/hQ7f03KwV8VpCeEFtvsC0yatEh20v/R0xlvpl+ochOOC6yi+4YkOqosvBHxvu0Pq8m4SVEHe1Y0/rjP1qS0+EnOeK6SNxXnUXpZNoAt42tkQYn9pETGS6CAOEwRqd+7SIuNz5edICdcq1PQmQfHDNAZHa/IUxW3p6QyvM9gCLJKcVxp/fBUFPDEzeJ7PCDlVnIlwvMXYdDlYAnhzNZRYVERjO0UOnsBMKTqv/iGY2pz6hojR3+m6WYMoRkxMmVU3N6rdJjaBxtpfejqjbzqBYvRDEIs0Ouo/HU0EsAnIpY7+uidq814OQ4QDDK1sAo2yv4onEkEIx41ypjPuhqrfnU10Ngvq3SZBK6oHmxfPVWuXBbfzfSFPwC2bnGL/rIhUIW4958a4TnWabTqq917nT3V4F5BKZ83E9DcOFyTo+u6xQS/D7ObGXj5zatsYbxKM3rw3rjf2t3AwdSatEa0NOXmKtv+QiOjg5I1xMnCw1uiRy3SG4RzeZzAq6xNJp3Q8WfTmRpUZBNVDx+1tElS6EbVth0FUdO2mdpjGWtoEGm3/rIgGrKf94C/TdiJq0uOkVbuslBXIoU2CPMeqsbkREbW0CdSRG2x/aRGhavZSGXU7dIqv7Ll78Ni7VPZUWPCOs5sElQcPRD63RoFQd67TfEFtwcEeD0pfZRMolsXaX1pEe8QLINMRncRCciwubCgW1B6jnoX0oU2C5OuUf83NjXDEsUaIsr+oiNQJ73Vc6+gvGgE49Z6qNGFBmwSFN2rzXk7QlivTiaoX7hB91U2gsfYXFRHkCpDXXXOtdhCW4E2Cyjv1YLIobLXrvaEo2qin8hj7zzzlW09iCqw6DbZOyAS+ajwV90QTBiYna4TsLUKTK/uJC9bk6Rg90U/c9cdp2iai4+y3plBvImqqO44TjFsT8Rph8ExCv7fF63H2aTHUY42oIaMZ54l4jtP/j/Oa/NvOxdjZKg5lgGdWezr5H4+e42nq4vGBAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}4 + 4 i & 2 - i\\\\-3 + 2 i & 4 + i\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "[4 + 4*I   2 - I]\n",
       "[               ]\n",
       "[-3 + 2*I  4 + I]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import sympy\n",
    "sympy.interactive.printing.init_printing(use_unicode=False, wrap_line=False)\n",
    "A = sympy.Matrix([[4 + sympy.I*4 , 2- sympy.I], [-3 + sympy.I*2 , 4 + sympy.I]])\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAEgAAAAQCAYAAAC1MDndAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAC40lEQVRYCc2X7VEbMRCGTYYC8tEBdACkgjgdhKQCcAdh+Od/DHRAqICPDqCDBDoIHSRxB87zyJIiH3fnc4yxd2a9q9VqtffernTeGI/HvUTD4XAH/RreRR8le1Uyd1qxXWJ7qNhebDgrb+a3SGYQE3qNdHyK/S7apgT2rxi2kYNNflxwAf+G92AX1xK+zgngEXoIjnQsb8MvRl3zjn7mmwDqoX8i0VvkPnxTk7S+wS5AVsq+TugiZxU1kUBYLSXyAvzYtKDNTpw+81vIb21+dXOs6Zr3IesP8ReQBEbK/5i5ZMvb4Jdf9ma2zlBYJOqC96F0xf6xHM+pC668TLL1BVMORM4jOI7aRWeACGPZGThv1B56PWbJ12p5U2aDzZctnU9E8+88AHk+PRLcKvoC/4ItxeuYBOr6E7n2ydJLZoCeWxvdSrYdfab7NDcPQAaQ9lh8NFHDufWH8QH8pJeTzzpI8vPFCs572Lb7AZd07HPB+njWBvBelR5NOosSODvoGfXof4W8KHyawqzUTn4P8BnshXQJWyWh1ZCC9z0m6Hy+dOapINfnhTGY4h62NG3BdDug/iMSsNd9M1V6q4H5fAUXDj5QuF0L27OoxL2BPUs9HjyfPDrSd9xnxidpo04AGQx2jUGbqO37qQ6AHjF9g17zZ01BF7UT2+pwrwRACmmL+dL6zIXjAWk+dkvukk4txgLJ6kitFgyVn7rqqrisZGiF205tuafEwncRvhZE8J8HINukrkp2sRuwtr3SziuUI/ZOLVWm4ZEghbwjIFZbuvoFq1cF6J1GKJwNE3XySwDL8A6Z/4fFoPbsQem7Ar0xb3Lxxr0tcyLv1Epe9QIohZfP2GdU/6lxwz+rGLzWJHvS0rJfbRk/z3M/MtZXgEL5IQXyBFu1vzHPJtYtdAaxvlPe+Plc5YEvAE/+rMZ4AUz08NwBoNmPshwPklgIoOVkNR212mLTs8sfWd6pxJe/23/s8Bedihl+S2jNugAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle 16 + 13 i$"
      ],
      "text/plain": [
       "16 + 13*I"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "det = A.det()\n",
    "sympy.simplify(det)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAEgAAAAOCAYAAACM7Fo2AAAACXBIWXMAAA7EAAAOxAGVKw4bAAACoklEQVRIDc2X21UUQRCGFw8BeMlAMgCMQMxANALYDOTwtm8czECJACEDMAKFDCQDdTNYv6+3q3ccppfZh73UObVVU13d888/Vd2zW5PJZBAyGo128a/QPfxxxNuWsfNW7JLYfSu2ssuncDP+GjDDDOg51utz4rc59p8h/onADna4zY8TLtA/6D7q5E4h1zEJPMFPi2O9VnfQlUlf3DlPvEHQAP89QG+wh+h1B2hzU1yCrJRDk/BlziqqiURYLU3mJfihNmFZ8QVwH4PhmHwJCTIC/yljEStQySsve7tEn3CYJOuS97aZSvxd83oDfVvfIlCTgHmM5qv5pjdBLGPZuXC50fylN2MUvFbLiyYaYr5s5cvU1H8XIcj96YHFraKP6G/UUrzKIHA3X8B6AEoPmSH+10CM71ZhO/pMdzG2CEEuoOwz+WTqpn3rL9dH6KNejpxNsODzxUrOG9S2+4k25dTnQs1xr03kPWtm1HwmBTm7+IX1nP8Ne9HIqS2z1jj47tHPqAfSJWqVpFbDSt6PDNDxcugsUkHOLxPzYpo71NK0BeN0wJ0JACT4OxpEzwbrnkfwUr6tWPcaHXNrtwf3J7eOuNcHrs8CVi+CXAx1jovWZN73k/P2ahOXGQe31THABgFxO1vMdjpgLG0PWCvKl1i6pFeLMUGxOpxck67qquWuMm6F207zsAee9F1ErgWR8hchyCOxq0qsDBfsbK+48xrtmHtHSzVhuCUoCXcmxGqLo1+yBm2CXhlEXk7N7JcFLMNbbPkflhe1Z49mmWvxqrhB44l700QF7mglj3oJVNLL59pn1P9lcMs/qwQ81hR70tKyX20ZP89LP3JtrgSl8sNK5Bmxdn8TXr5w3164yfO50t+pjEoCHv1ZzeslMvHTc/8DB30DQW5HOqgAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle 16 - 13 i$"
      ],
      "text/plain": [
       "16 - 13*I"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "detH = A.H.det()\n",
    "sympy.simplify(detH)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANsAAAAzCAYAAAAelH2vAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAKt0lEQVR4Ae2dX5LUNhDGZ6l9yCNZHvMGNyDhBCw3YMkJAjeA4gneKHIDyAkSuAGbEwS4AeQEUHuD5Ptp1CrZ67Fkez3WaK0qj/7L/XWr1ZJsa45evHhxc9PhXr58edGRHJKUf1vX15DQCqTyW8WriQo3/Dzp4001YK8xkD45+7xL3DlWyiddt1s5bxV/0kprR5+pUeq90dVQOqV/Vloqv93eonGP5VREfNcFrnu6XnksCm6d4q8V+uajF4rDq9g9V4T83+PEQw0LH4PHI11PFP65jUNpT5V2SxeYbyn+rF2m5Ljo/RTjUhi8yBBnevGb0tvGp0/O71SXvhS7c5QN90CNnW+D2b8QQoPtRk1RU/nZN9pTQZToq/jglEc+nehvXT/a/ZXGwATjP+u6qzDxhrIp/aA6m2Hr8j1Gky+dsOGUT6eCZw6z/Iek6TprFCw0IjqRMXKM3WulB0OjMMYEOd+JCyl9p5yV96BV9rHiZzfixIFhmHwUX6r/VnEjNJU/8HazF3+lO/wZ3YXRGivnnHChjB/lY7U33m8w1RWs6AeMurDQjZlLBPGhwjHPGLBRuEuKGdUpIigaGUS6cD32eUYncmfJ1FZKy8/2j7NLXi4IEcGJGOKxtqfyQ90SAqLfKVFECx0pxuBGJ5XDYsP4c4XDbEBhhEf6PYUPYmQXraOd5wP1w/RKaUyrSftFV+ANCSU50chggKwasxJPI8biow9f8lR3tJxHWzbdNIwKCtMxP8iPGd+bfwlFIQnC8FQX04Y38m1KiXC4sGT4dCSmGyjgRj5pd+VjBQ5iZIfuKU5Yke+FLgYf55Rmo39Is7zCfKxX55pa6czOwGUO5WOWhpWfJOfRymaUeAJYPHeOZKl8a6cUX/QihPu6WMcyp8edbL3NF6XBdISBFUch6VjkIyQGHQQTC0tJ1TosOBthN7kUxqLhwvR7Gy3n18vofQ5FKsvggUxtY2iSnCcrmwh5pKs9BVNScKn8ULCUgJhsyoT1wmKZlQ7TC18Gkk/J93FGQRbU18IJMwMsCsfUinWNTcv6+sNg3qhdeOxmEYMrRxXUhhsY5Zs8o9xm0JdlGfGzwm7wpJ4Pj5LzcfMWo2KpG6fyR930KiuJgYzK/+q6r7B1FMdgpTFak4aAzMIpGJwTnG+DTnemMApq7YSCtQXAKUx0QGcp5IOfeLIzD+QF8uGa6mjjjuiztbhr08f/kW84UEosttsAk098Ix9s1Bkl5xs0MtHB8D6Xyu+ru5c8MfBCN/quK+4kMJT0v3ThEFDY+FAd62g2faa8TTEZYK6DgyfMXMw1eGSJpfiSGfJBidwlukzpeJ4aKxqzE5YIDJrImSUD/QM3Ws7H2/qTfumgYXrV0VIqv6PKIkko0nMx95u/+z358RSCNRmbJwiIMjx3sbm8gm7T5FfKKGxCJP1gnbAwojNw0MGYJjrs8m1zgTjpTPHgh3sGKb9452m2RzcsF9jgQ+HYHMN64QenPBtAGVxHyflIr2t9UeWdGxzhbmtg5cBCHFBHZ5MCpTYlX4iScbcV3e6x0VVMI8dRsNZaOXDNOJCcRkor/5ubJ7rHUd89lI9Z59Up/FzHRkVyk6IEfLmAdpW7Cv6UwAfRwFqJKWvbnZCgfJvKxfmsw8JaOs6wsPIX78PQsk4jTSKrXywHpCzrNLJY6ayErRwokAPJaWQXzRppBn2WoPK2q/XBt0f8u9LddmvXPZZMqx1fLm9r5sMS2AYrm4hka7v97Cz1WQLKxY4MdXkU8F7t2NsGipbjaseXy+ma+bAUtkG7kSKSxWv84Ndk99jnWdyev8RKydsZfJLDE/z46wCrs7jvMVSLL5fBNfNhSWzZlk1EshOI8nRZJHaJ+h5sOzl7oLTD5ykXLrGQn9rx5bK5Zj4sjW2IZcN6dT5UVPrOzxIiIT9XGKvBE/g/VIcdppJc7fhyeV0iHxiYr2JwXhRblmXzipG1maGyWD8UKbzKpDQUzN4h3Cj+SnFeh+l9vqb8vTjRA73V4stlYql8EF2h7+RiaZcrAVvSsolINjdO5HetZRqYfFnWa+GdQgoovW3FXFtKZw24qPM0V4svl7k186EUbDmWjTXWpM8SVJ9DYFBAe6PjxHeCpALndpYJ5WrHl8uamvlQBLaksnkFMSXZKI6lYxufzxIukKRP41UbPl2wHUg2TWzXkU2VWLHctE1l4zQV2b8TDWCrFl8uR2vmQynYksoWC0tEo2RjPktA6ezzFU6t2qit3vfZKLNPJ3oYRJjWTsX3k9r4oTR8Q3gp2muW82LYRr8bKYFgmtlhxNFRcXzPZNauN39bvJxf0Y1lBgd+w+Iqz1k+j5mPJfkkKWwAKR6c0tvffIW8EgIpDKl8MBhGj6froNpFofZh8Hk7+60Rbhjld+7AW7kcX22g4GeDLFur4dRbI6n8VnOLR82yYd1ixxQY5WJ6bHkMJJ1O5Wzq3Jm/ZGIKQyof2lWGXeTeg2pLxijasvqlcF65HG9MYAzPLKzz0Uz7rZFU/oRbz1KV8yV2HjqrPD7lYJRrWL1ZKJmp0RSGjHxkXPRBtSkMon+xfjnFsqXeGknlz9SlRjdLRwpOQiOePbqpPAMP1q/mQ1rddEhYmQWAtXFQreKH4Hr75ZxyHG3ZRFTvWyOp/NKkInqDxVKY3dLGobN99Ko800oOh8HyVXlIq8cITjaQ8HnQzJQMBTwYJ3p39luPcTY5jla2mLsiklGODrpr06A3P25r6bBn+NAzWU5EN0KEB0xHL5bGMcP9wYjbdVDtNveAfjv67axynKxsIpgpBVOuxlsjxvNUvpUryGe3MTx3y6FLGE3BmKKwm1mdA6MHFV44V5oNKvHafTB2tXOqa68WUve71G+VNqscp6zZNp5gHmS7Z1M+TroTjI/vzB8slf1UGKUwwsrUik5X8yGtyNUsXCwNU8Q4bUgY3nHtxfX1yznlONqyeYIZxTsPs0zl74Wr425yd1w1p2jsWDLao7A1OmYw4WUEYYVXWIPJLwrvi1kZ/ZIBcxY5TrFsqcMsU/n74u/Q+zBKh6mSVfZCQokQRteBpXS4UYd32j3m9lMYMvJZl/YdVDs3hGT7KQxqINUvZ5Pj6DdIkqjXAisHMjkgBWFj6WAPYU3BFD73yGT0NDJ1gzV/5cDKgSYHktNIaeXiB1yKBhbP6yGtTdkNihUiR9b4TMPbzm26iMautS7rp7BObFckXgK2Lrraaes0ss2RNb53DkhZ1mnk3rm+3nDlQMUcSE4ju7BrJFoPaW1+TsQDUqZAB3EIrcm0ZjmWiG2wsgnEekjr9ts3tpD5TzIcysaOU/GH0EIsrmY5lopt0G6kQLC47Xpb4LHPc4LUDw8/2cqNHxCvh7Qadxb2a5ZjydiyLZtAsCOI8qyHtO5QFi9o+FTcIbRGcs1yLB3bEMuG9er8RFzpOz9bMCHL51N0rCJP6NdDWsWEhVyJcuQVN66prkRsAVOWZZMysTVb7SGmteMzaZeKU3RNfreyVGzGe/ykZRMIFv/VHmJaOz4Tds04DwVbjmVjDbIe0roddHZ9LlTyIbSmbzXL8SCwmbKxcwjBwSl+QUR+1YeY1o7PBFozztKwiZ6GLkkGLm7K5r5LM8HId8e3RXEXVCM8S5p6iGmRh7QCsHZ8ToiV4yxEhu/E5lPjt/fP/wdxe2oajwRGlwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{77}{425} - \\frac{36 i}{425} & - \\frac{19}{425} + \\frac{42 i}{425}\\\\\\frac{22}{425} - \\frac{71 i}{425} & \\frac{116}{425} + \\frac{12 i}{425}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "[ 77   36*I     19   42*I]\n",
       "[--- - ----  - --- + ----]\n",
       "[425   425     425   425 ]\n",
       "[                        ]\n",
       "[ 22   71*I   116   12*I ]\n",
       "[--- - ----   --- + ---- ]\n",
       "[425   425    425   425  ]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ainv = A.inv()\n",
    "sympy.simplify(Ainv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAADUAAAAzCAYAAAA3v0BtAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAC2ElEQVRoBe1a7U3jQBAFdAUg6CB0wEcFQAeEDu4oI/mHKAGoAIUOoAMEHUAHx6WEe8/n0a2dmbG9Cc4QeaTVfu+8tzO73oyyPZlMdrcUmU6nc6U5VBMwqth3gPIV6U8t3YRCb4OZ1XCTx2wblnpH4Qqsn+2536MHHH4B6ZiW2jj50ZYRduEQY2nuI5R7O2/QJUfhN3QfIN2g7QO5KS4pTOZBvEf6RDpGGiH1JtDP836N/JFKSzyvyM+RTGKu+2HiHGmMdIU1H7hwXwKdPB+7yAtC1IsyPYT1W9YtcUlZk3pqH0PPm6LrBW1nIKhe5xwfmdQZ8NHt6yJux35VQpLyrJCw2EvKlWJIUkAogL1b9lu6X2X3lcq+0lY0RbWUdpaEg1iR3y1VQpIqr24C1lxM2uTCWCAWklSJkm/R0QLi/+fNfKtGJsUnGV8xdTlCw1tizXp/p++UHEzx6YXFVtkA0HdY7xP5hayLMl3vEumntGm5+/bjBCzEHaPIx26GNvrzE3Iq/kqhVfiAPUHOi4H5KeraSwNd/6QNKT5X1iIAP4divjs7SeQz1YlIOnggle5G5PJgqcjWSbENlkp3I3J5sFRk66TYNtJSjc8k7gCeK50DiunOLVuG/kOs0TqQ2kgKC2YFFFdAhC/yrECq634glB1QXAGp7ECqSwrAsgOKy5JaZn4TKf6G0oIgEh9gfzgxScH1JMDhge7lV7AHQOszSWGwAOYPNUvaELfmflm7R6qNUolbtBnb2xiPlHaWBJhY0QwoysB15CYpnClxO83FpE0ujHVgN3WapMoZ2QFFU2MPHU2ksgOKPWA3Vbik4ILZAUVTY16HXEhylt1VGt9+mJ0VUHS1tuzEpmYFUhtJlRdG54BiS9zuMOjOCqS67udqDNw5kApsnAq0wVKV7QhckdtvhJtGnj4F3PLWCwy9iJ1UMANsURf343990j8ySqAlNCmA43dsAfdf2JDTgoeVBLwAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0\\\\0 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "[1  0]\n",
       "[    ]\n",
       "[0  1]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.simplify(A.multiply(A.inv()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 2 \n",
    "\n",
    "The diagonal entries must be real since $A_{ii} = \\bar{A}_{ii}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 3\n",
    "\n",
    "\\begin{align}\n",
    "\\det \\left(\\boldsymbol{Q} - \\lambda \\boldsymbol{I} \\right) &=\n",
    "\\det\\left(\n",
    "\\begin{bmatrix}\n",
    "\\cos \\theta - \\lambda & - \\sin \\theta \\\\ \n",
    "\\sin \\theta           & \\cos \\theta - \\lambda\n",
    "\\end{bmatrix} \\right) \n",
    "\\\\\n",
    "&= (\\cos \\theta - \\lambda)^{2} + \\sin^{2}\\theta \\\\\n",
    "&= \\lambda^{2} - (2 \\cos \\theta) \\lambda  + 1 \\\\\n",
    "&= 0\n",
    "\\end{align}\n",
    "\n",
    "Computing roots,\n",
    "\n",
    "\\begin{align}\n",
    "\\lambda &= \\cos \\theta \\pm \\sqrt{\\cos^{2}\\theta - 1} \\\\\n",
    "        &= \\cos \\theta \\pm i \\sin\\theta\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 4\n",
    "\n",
    "a. If $\\boldsymbol{x}$ is an eigenvector and $\\lambda$ is an eigenvalue,\n",
    "\n",
    "   $$\n",
    "   \\boldsymbol{M} \\boldsymbol{x} = \\lambda \\boldsymbol{x}\n",
    "   $$\n",
    "\n",
    "   Premultiplying by $\\boldsymbol{x}^{H}$, \n",
    "\n",
    "   $$\n",
    "   \\boldsymbol{x}^{H} \\boldsymbol{M} \\boldsymbol{x} = \\lambda \\boldsymbol{x}^{H} \\boldsymbol{x}\n",
    "   $$\n",
    "\n",
    "   Taking the complex conjugate of both sides and noting that $\\boldsymbol{x}^{H} \\boldsymbol{x}$ is real,\n",
    "\n",
    "   $$\n",
    "   \\left(\\boldsymbol{x}^{H} \\boldsymbol{M} \\boldsymbol{x}\\right)^{H} = \n",
    "   \\boldsymbol{x}^{H} \\boldsymbol{M}^{H} \\boldsymbol{x} = \n",
    "   \\bar{\\lambda} \\boldsymbol{x}^{H} \\boldsymbol{x}\n",
    "   $$\n",
    "\n",
    "   Since $\\boldsymbol{M}$ is Hermitian ($\\boldsymbol{M} =\\boldsymbol{M}^{H}$), \n",
    "   comparing the above two equations they can hold ony if $\\lambda = \\bar{\\lambda}$, \n",
    "   which is true only if $\\lambda$ is *real*. \n",
    "\n",
    "b. Consider two eigenpairs $(\\lambda_{1}, \\boldsymbol{x}_{1})$ and $( \\lambda_{2}, \\boldsymbol{x}_{2})$, where $\\lambda_{1} \\ne \\lambda_{2}$. We have:\n",
    "\n",
    "   \\begin{align}\n",
    "   \\boldsymbol{M} \\boldsymbol{x}_{1} &= \\lambda_{1} \\boldsymbol{x}_{1} \\\\\n",
    "   \\boldsymbol{M} \\boldsymbol{x}_{2} &= \\lambda_{2} \\boldsymbol{x}_{2}\n",
    "   \\end{align}\n",
    "\n",
    "   Multiplying by $\\boldsymbol{x}_{2}$ and $\\boldsymbol{x}_{1}$, respectively,\n",
    "\n",
    "   \\begin{align}\n",
    "   \\boldsymbol{x}_{2}^{H} \\boldsymbol{M} \\boldsymbol{x}_{1} &= \\lambda_{1} \\boldsymbol{x}_{2}^{H}   \n",
    "   \\boldsymbol{x}_{1} \\\\\n",
    "   \\boldsymbol{x}_{1}^{H} \\boldsymbol{M} \\boldsymbol{x}_{2} &= \\lambda_{2} \\boldsymbol{x}_{1}^{H} \n",
    "   \\boldsymbol{x}_{2}\n",
    "   \\end{align}\n",
    "\n",
    "   Taking the complex conjugate transpose of the first equation and exploiting that $\\boldsymbol{M} = \n",
    "   \\boldsymbol{M}^{H}$,\n",
    "\n",
    "   \\begin{equation}\n",
    "   \\left(\\boldsymbol{x}_{2}^{H} \\boldsymbol{M} \\boldsymbol{x}_{1}\\right)^{H} =  \\boldsymbol{x}_{1}^{H}    \n",
    "   \\boldsymbol{M} \\boldsymbol{x}_{2} = \\lambda_{1} \\boldsymbol{x}_{1}^{H} \\boldsymbol{x}_{2} \n",
    "   \\end{equation}\n",
    "\n",
    "   Comparing this to $\\boldsymbol{x}_{1}^{H}\\boldsymbol{M} \\boldsymbol{x}_{2} = \\lambda_{2} \n",
    "   \\boldsymbol{x}_{1}^{H}\\boldsymbol{x}_{2}$, since $\\lambda_{1} \\ne   \n",
    "   \\lambda_{2}$ both equations can hold only if $\\boldsymbol{x}_{1}^{H} \\boldsymbol{x}_{2} = 0$, i.e. the   eigenvectors are *orthogonal*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 5\n",
    "\n",
    "To show that $\\boldsymbol{A}^{H} \\boldsymbol{A}$ is positive semi-definite:\n",
    "\\begin{equation}\n",
    "\\boldsymbol{x}^{H} \\boldsymbol{A}^{H} \\boldsymbol{A} \\boldsymbol{x}  =\n",
    "(\\boldsymbol{A}\\boldsymbol{x})^{H} (\\boldsymbol{A} \\boldsymbol{x})\n",
    "\\ge 0\n",
    "\\end{equation}\n",
    "If $(\\lambda, \\boldsymbol{x})$ is an eigenpair of $\\boldsymbol{A}^{H} \\boldsymbol{A}$,\n",
    "\n",
    "\\begin{equation}\n",
    "  \\boldsymbol{A}^{H} \\boldsymbol{A} \\boldsymbol{x} \n",
    "  = \\lambda \\boldsymbol{x}\n",
    "\\end{equation}\n",
    "\n",
    "Premultiplying both sides by $\\boldsymbol{x}^{H}$:\n",
    "\n",
    "\\begin{equation}\n",
    "  \\boldsymbol{x}^{H}  \\boldsymbol{A}^{H} \\boldsymbol{A} \\boldsymbol{x} \n",
    "  = \\lambda \\boldsymbol{x}^{H} \\boldsymbol{x}\n",
    "\\end{equation}\n",
    "\n",
    "Since $\\boldsymbol{x}^{H} \\boldsymbol{A}^{H} \\boldsymbol{A} \\boldsymbol{x} \\ge 0$ and $\\boldsymbol{x}^{H} \\boldsymbol{x} \\ge 0$, all eigenvalues $\\lambda$ must be greater than or equal to zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 6\n",
    "\n",
    "a. Eigenvalues of $\\boldsymbol{A}$ satisfy $\\det(\\boldsymbol{A} - \\lambda\\boldsymbol{I}) = 0$. \n",
    "   For the transpose, we have \n",
    "   \\begin{equation}\n",
    "     \\det\\left(\\boldsymbol{A}^{T} - \\lambda\\boldsymbol{I}\\right) = \n",
    "     \\det\\left(\\left(\\boldsymbol{A} - \\lambda \\boldsymbol{I}\\right)^{T}\\right) = \n",
    "     \\det(\\boldsymbol{A} - \\lambda \\boldsymbol{I}),\n",
    "   \\end{equation}\n",
    "   since the $\\det\\boldsymbol{A} = \\det \\boldsymbol{A}^{T}$. Hence eigenvalues of $\\det\\boldsymbol{A}$ and \n",
    "   $\\det \\boldsymbol{A}^{T}$ are the same.\n",
    "  \n",
    "b. Noting that $\\det\\left(\\boldsymbol{A}^{H}\\right) = \\det\\left(\\overline{\\boldsymbol{A}}^{T}\\right) = \\overline{\\det\\left(\\boldsymbol{A}^{T}\\right)} = \\overline{\\det\\left(\\boldsymbol{A}\\right)}$,\n",
    "   \\begin{equation}\n",
    "   \\det\\left(\\boldsymbol{A}^{H} - \\lambda\\boldsymbol{I}\\right) = \n",
    "   \\det\\left(\\left(\\boldsymbol{A} - \\bar{\\lambda} \\boldsymbol{I}\\right)^{H}\\right) = \n",
    "   \\det\\left( \\overline{(\\boldsymbol{A} - \\bar{\\lambda}\\boldsymbol{I})^{T}}\\right) = \n",
    "   \\overline{\\det(\\boldsymbol{A} - \\bar{\\lambda} \\boldsymbol{I}}).\n",
    "   \\end{equation}\n",
    "   \n",
    "c. If $\\lambda$ and $\\boldsymbol{x}$ are an eigenvalue and eigenvector, respectively, of \n",
    "   $\\boldsymbol{A}\\boldsymbol{B}$, then\n",
    "   \\begin{equation}\n",
    "    \\boldsymbol{A}\\boldsymbol{B} \\boldsymbol{x} = \\lambda \\boldsymbol{x}\n",
    "   \\end{equation}\n",
    "   Multiplying $\\boldsymbol{x}$ by $\\boldsymbol{B}$ on both sides,\n",
    "   \\begin{equation}\n",
    "     (\\boldsymbol{B}\\boldsymbol{A}) (\\boldsymbol{B} \\boldsymbol{x}) \n",
    "     = \\lambda (\\boldsymbol{B}\\boldsymbol{x})\n",
    "   \\end{equation}\n",
    "   We see that $\\lambda$ is an eigenvalue of $\\boldsymbol{B}\\boldsymbol{A}$, and the eigenvector is now    \n",
    "   $\\boldsymbol{B} \\boldsymbol{x}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 7\n",
    "\n",
    "Follows directly from the definition of the norms: \n",
    "\n",
    "\\begin{equation}\n",
    "  \\max_{i=1}^{n} |x_{i}| \\le \\sum_{i=1}^{n} |x_{i}| \\le n \\max_{i=1}^{n} |x_{i}|.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 8\n",
    "\n",
    "Recall that from the defintion of a matrix operator norm it follows that $\\| \\boldsymbol{A}\\boldsymbol{x} \\| \\le \\| \\boldsymbol{A} \\| \\| \\boldsymbol{x} \\|$ \n",
    "\n",
    "a. From the definition of the matrix norm:\n",
    "\n",
    "   \\begin{equation}\n",
    "      \\| \\boldsymbol{A} \\boldsymbol{B} \\boldsymbol{x} \\| \n",
    "      \\le \\| \\boldsymbol{A} \\| \\| \\boldsymbol{B} \\boldsymbol{x}\\| \n",
    "      \\le \\| \\boldsymbol{A} \\| \\| \\boldsymbol{B}\\| \\| \\boldsymbol{x}\\|\n",
    "   \\end{equation}\n",
    "   \n",
    "   for all $\\boldsymbol{x}$. Re-arranging,\n",
    "\n",
    "   \\begin{equation}\n",
    "      \\frac{\\| \\boldsymbol{A} \\boldsymbol{B} \\boldsymbol{x} \\|}{\\|\\boldsymbol{x}\\|} \n",
    "      \\le \\| \\boldsymbol{A} \\| \\| \\boldsymbol{B}\\|\n",
    "   \\end{equation}\n",
    "\n",
    "   for all $\\boldsymbol{x} \\ne \\boldsymbol{0}$. From the definition of the norm, $\\|    \n",
    "   \\boldsymbol{A} \\boldsymbol{B} \\|$ is the largest possible value of $\\| \\boldsymbol{A} \n",
    "   \\boldsymbol{B} \\boldsymbol{x} \\| / \\|\\boldsymbol{x}\\|$ (over all\n",
    "   $\\boldsymbol{x}$, exluding the zero vector), hence\n",
    "\n",
    "   \\begin{equation}\n",
    "      \\| \\boldsymbol{A} \\boldsymbol{B} \\|  \\le \\| \\boldsymbol{A} \\| \\| \\boldsymbol{B}\\|\n",
    "   \\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b. Consider\n",
    "   \\begin{equation}\n",
    "      \\| (\\boldsymbol{A} + \\boldsymbol{B}) \\boldsymbol{x} \\| \n",
    "      = \\| \\boldsymbol{A}\\boldsymbol{x} + \\boldsymbol{B} \\boldsymbol{x} \\| \n",
    "   \\end{equation}\n",
    "\n",
    "   From the triagle inequality for vectors,\n",
    "\n",
    "   \\begin{align}\n",
    "      \\| \\boldsymbol{A}\\boldsymbol{x} + \\boldsymbol{B} \\boldsymbol{x} \\| \n",
    "      &\\le\n",
    "      \\| \\boldsymbol{A}\\boldsymbol{x} \\| + \\| \\boldsymbol{B} \\boldsymbol{x} \\|  \\\\\n",
    "      &\\le\n",
    "      \\| \\boldsymbol{A} \\| \\| \\boldsymbol{x} \\| + \\| \\boldsymbol{B} \\| |\\boldsymbol{x} \\| \n",
    "   \\end{align}\n",
    "\n",
    "   Re-arranging\n",
    "\n",
    "   \\begin{equation}\n",
    "      \\frac{\\| (\\boldsymbol{A} + \\boldsymbol{B}) \\boldsymbol{x} \\|}{\\| \\boldsymbol{x} \\|}\n",
    "      \\le\n",
    "      \\| \\boldsymbol{A} \\| + \\| \\boldsymbol{B} \\|  \n",
    "   \\end{equation}\n",
    "\n",
    "   Using same argument as from part a), we get $\\| \\boldsymbol{A} + \\boldsymbol{B} \\| \\le  \n",
    "   \\| \\boldsymbol{A} \\| + \\| \\boldsymbol{B} \\|$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 9\n",
    "\n",
    "Recall that the smallest possible value of $C$ for which $\\|{\\boldsymbol{A}\\boldsymbol{x}}\\| \\le C \\| \\boldsymbol{x}\\|$ holds is the operator norm $\\|\\boldsymbol{A}\\|$. The task is to find $C$ for the different norms. It simplifies the proofs if we consider all vectors for which $\\|\\boldsymbol{x}\\| = 1$ (this is possible since $\\|\\alpha \\boldsymbol{x}\\| = |\\alpha| \\|\\boldsymbol{x}\\|$), in which case the task is to find the smallest possible $C$ such that $\\|\\boldsymbol{A}\\boldsymbol{x}\\| \\le C$.\n",
    "\n",
    "What we need to do is (i) prove the inequality, and then (ii) find an equality to show that is is a weak inequality. \n",
    "\n",
    "a. For the 1-norm, denoting the $j$th column of $\\boldsymbol{A}$ by\n",
    "   $\\boldsymbol{a}_{j}$, and for a vector $\\|\\boldsymbol{x}\\|_{1} = 1$:\n",
    "\n",
    "   \\begin{align*}\n",
    "        \\| \\boldsymbol{A} \\boldsymbol{x}\\|_{1}\n",
    "        = \\|\\sum_{j=1}^{n} \\boldsymbol{a}_{j} x_{j}\\|_{1}\n",
    "        \\le\n",
    "        \\sum_{j=1}^{n} \\|\\boldsymbol{a}_{j}\\|_{1} |x_{j}|\n",
    "        \\le \\max_{j=1}^{n} \\|\\boldsymbol{a}_{j}\\|_{1}\n",
    "   \\end{align*}\n",
    "\n",
    "   This is the maximum column sum.\n",
    "\n",
    "   The term on the right could possibly be larger than $\\| \\boldsymbol{A} \\|_{1}$, whereas the norm is the smallest \n",
    "   possible value for the RHS that still satisfies the inequalities.  If we can show a case for which equality is\n",
    "   reached, $\\|\\boldsymbol{A} \\boldsymbol{x}\\|_{1} = \\max_{j=1}^{n} \\|\\boldsymbol{a}_{j}\\|_{1}$, we have the norm.\n",
    "   For the vector the $\\boldsymbol{e}_{j}$ with $e_{j} = 1$, where $j$ is the column with the greatest 1-norm,\n",
    "   and $e_{i \\ne j} = 0$, we have\n",
    "   \\begin{equation}\n",
    "        \\|\\boldsymbol{A} \\boldsymbol{e}_{j}\\|_{1} = \\max_{j=1}^{n} \\|\\boldsymbol{a}_{j}\\|_{1}\n",
    "   \\end{equation}\n",
    "   Therefore,\n",
    "   \\begin{equation}\n",
    "      \\|\\boldsymbol{A}\\|_{1} = \\max_{j=1}^{n} \\|\\boldsymbol{a}_{j}\\|_{1}\n",
    "   \\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b. For a vector $\\|\\boldsymbol{x}\\|_{\\infty} = 1$:\n",
    "\n",
    "   \\begin{align*}\n",
    "        \\|\\boldsymbol{A} \\boldsymbol{x}\\|_{\\infty}\n",
    "          = \\max_{i=1}^{m} |\\sum_{j=1}^{n} a_{ij} x_{j}|\n",
    "          &\\le \\max_{i=1}^{m} \\sum_{j=1}^{n} |a_{ij}| |x_{j}|\n",
    "          \\\\\n",
    "          &\\le \\max_{i=1}^{m} \\sum_{j=1}^{n} |a_{ij}|\n",
    "   \\end{align*}\n",
    "\n",
    "   The RHS is the maximum row sum.\n",
    "\n",
    "   As before, we need to find a case with equality.  If the row\n",
    "   with the maximum sum is row $k$, we choose a vector $\\boldsymbol{x}$\n",
    "   where $x_{i} = \\pm 1$ such that the sign of $x_{i}$ is the\n",
    "   same as the sign of the entry $a_{ki}$. We then have\n",
    "\n",
    "   \\begin{equation}\n",
    "      \\|\\boldsymbol{A}\\|_{\\infty}\n",
    "      = \\max_{i=1}^{m} \\sum_{j=1}^{n} |a_{ij}|\n",
    "   \\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 10\n",
    "\n",
    "Recall that the $l_{2}$ norm is the square root of the maximum eigenvalue of $\\boldsymbol{A}^{T}\\boldsymbol{A}$ (matrix is real in this question). We have\n",
    "$$\n",
    "\\boldsymbol{A}^{T}\\boldsymbol{A} =\n",
    "\\begin{bmatrix}\n",
    "26 & 5 \\\\ 5 & 1\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "For this matrix, $\\lambda = (27 \\pm \\sqrt{27^{2} -4})/2$. Hence $\\|\\boldsymbol{A}\\|_{2} \\approx 5.1926$.\n",
    "\n",
    "Checking for a vector of length $1$:\n",
    "$$\n",
    "\\boldsymbol{A}\n",
    "\\begin{bmatrix} \\cos \\theta \\\\ \\sin \\theta   \\end{bmatrix} =    \n",
    "\\begin{bmatrix} \\cos \\theta \\\\ 5\\cos \\theta + \\sin \\theta  \\end{bmatrix}\n",
    "$$\n",
    "Comptuing the norm,\n",
    "$$\n",
    "F(\\theta) = \\|\\boldsymbol{A}\\boldsymbol{x}\\|_{2}^{2} = 25 \\cos^{2} \\theta + 10 \\cos\\theta\\sin\\theta  + 1\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "To find the extreme points of the function, we differentiate $F$ with respect to $\\theta$  and set the derivative equal to zero:\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\frac{d F}{d \\theta} &= -50 \\cos\\theta \\sin\\theta + 10 (\\cos^{2} \\theta - \\sin^{2} \\theta) \\\\\n",
    "&= -25 \\sin 2\\theta + 10 \\cos 2\\theta \\\\\n",
    "&= 0\n",
    "\\end{align*}\n",
    "$$\n",
    "The norm is maximised when $\\theta = \\arctan(2/5) /2$. Evaluating:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "θ = 0.19025318855618245\n",
      "F = 5.192582403567252\n"
     ]
    }
   ],
   "source": [
    "θ = np.arctan(2/5)/2\n",
    "print(\"θ =\", θ)\n",
    "\n",
    "F = 25*np.cos(θ)*np.cos(θ) + 10*np.cos(θ)*np.sin(θ) + 1\n",
    "print(\"F =\", np.sqrt(F))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 11 (condition number)\n",
    "\n",
    "Recall that $\\kappa(\\boldsymbol{A}) = \\| \\boldsymbol{A} \\| \\| \\boldsymbol{A}^{-1} \\|$.\n",
    "\n",
    "We have $\\| \\boldsymbol{A}\\|_{1} = 3$ and $\\| \\boldsymbol{A}\\|_{\\infty} \\approx 2$.\n",
    "Note that\n",
    "$$\n",
    "\\boldsymbol{A}^{-1} \\approx -\\frac{1}{2}\n",
    "\\begin{bmatrix}\n",
    "1 & -2 \\\\\n",
    "-1 & 10^{-4} \n",
    "\\end{bmatrix}\n",
    "$$ \n",
    "so $\\| \\boldsymbol{A}^{-1} \\|_{1} \\approx 1$ and $\\| \\boldsymbol{A}^{-1} \\|_{\\infty} = 3/2$.\n",
    "Therefore $\\kappa_{1}(\\boldsymbol{A}) \\approx 3$ and $\\kappa_{\\infty}(\\boldsymbol{A}) \\approx 3$.\n",
    "\n",
    "Recall that $\\kappa_{2}(\\boldsymbol{A}) = \\sqrt{\\lambda_{\\max}(\\boldsymbol{A}^{T}\\boldsymbol{A}})/\\sqrt{\\lambda_{\\min}(\\boldsymbol{A}^{T}\\boldsymbol{A}})$. As an approximation we ignore the $(1, 1)$ entry:\n",
    "$$\n",
    "  \\boldsymbol{A}^{T}\\boldsymbol{A} \\approx\n",
    "  \\begin{bmatrix}\n",
    "     1 & 1 \\\\ 1 & 5\n",
    "  \\end{bmatrix}\n",
    "$$\n",
    "Computing the eigenvalues, $\\lambda \\approx 3 \\pm \\sqrt{5}$, therefore $\\kappa_{2} \\approx \\sqrt{3 + \\sqrt{5}}/\\sqrt{3 - \\sqrt{5}} \\approx 2.618$. \n",
    "\n",
    "The matrix is very well conditioned, however LU factorisation may require pivoting. This is an issue with LU factorisation rather than a pathological problem with the matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 12 (least squares)\n",
    "\n",
    "A projection matrix $\\boldsymbol{P}$  has the property $\\boldsymbol{P} = \\boldsymbol{P}\\boldsymbol{P}$, and $\\boldsymbol{P}^{H} = \\boldsymbol{P}$.\n",
    " \n",
    "a. The solution to the least squares problem is $\\hat{\\boldsymbol{x}} = \\boldsymbol{A}\n",
    "   (\\boldsymbol{A}^{H} \\boldsymbol{A})^{-1} \\boldsymbol{A}^{H} \\boldsymbol{b}$.\n",
    "   Therefore $\\boldsymbol{r} = \\boldsymbol{A} \\hat{\\boldsymbol{x}} - \\boldsymbol{b}\n",
    "   = \\boldsymbol{A}(\\boldsymbol{A}^{H} \\boldsymbol{A})^{-1} \\boldsymbol{A}^{H}\n",
    "   \\boldsymbol{b} - \\boldsymbol{b}$. Insert this expression for $\\boldsymbol{r}$ into the expression in the\n",
    "   question and re-arrange. to show the result.\n",
    "\n",
    "   Vectors $\\boldsymbol{A}\\boldsymbol{z}$ lie in the *column space* of $\\boldsymbol{A}$,\n",
    "   hence the expression says that the least-squares residual is *orthogonal* to the column\n",
    "   space of $\\boldsymbol{A}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b.\n",
    "   $$\n",
    "     \\boldsymbol{P}\\boldsymbol{P}\n",
    "        = \\boldsymbol{A} (\\boldsymbol{A}^{H}\\boldsymbol{A})^{-1}\\boldsymbol{A}^{H}\n",
    "    \\boldsymbol{A}(\\boldsymbol{A}^{H}\\boldsymbol{A})^{-1}\\boldsymbol{A}^{H} \n",
    "     = \\boldsymbol{A} (\\boldsymbol{A}^{H}\\boldsymbol{A})^{-1}\\boldsymbol{A}^{H}\n",
    "    = \\boldsymbol{P}\n",
    "   $$\n",
    "   and\n",
    "   $$\n",
    "        \\boldsymbol{P}^{H}\n",
    "        = \\boldsymbol{A}(\\boldsymbol{A}^{H}\\boldsymbol{A})^{-H}\\boldsymbol{A}^{H}\n",
    "        = \\boldsymbol{A}(\\boldsymbol{A}^{H}\\boldsymbol{A})^{-1}\\boldsymbol{A}^{H}\n",
    "   $$\n",
    "   by $\\boldsymbol{A}^{H}\\boldsymbol{A}$ being Hermitian\n",
    "\n",
    "c. We can phrase a least squares problem as\n",
    "   $$\n",
    "        \\boldsymbol{A} \\hat{\\boldsymbol{x}}\n",
    "        = \\boldsymbol{A}(\\boldsymbol{A}^{H}\\boldsymbol{A})^{-1}\\boldsymbol{A}^{H} \n",
    "        \\boldsymbol{b}\n",
    "        = \\boldsymbol{P} \\boldsymbol{b}\n",
    "   $$\n",
    "   which says that $\\boldsymbol{P}$ projects $\\boldsymbol{b}$ into the column space of $\\boldsymbol{A}$. If\n",
    "   $\\boldsymbol{b}^{\\prime} = \\boldsymbol{P}\\boldsymbol{b}$ is in the column space of $\\boldsymbol{A}$, then\n",
    "   $\\boldsymbol{A} \\hat{\\boldsymbol{x}} = \\boldsymbol{b}^{\\prime}$ has a solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 13 (pseudo inverse)\n",
    "\n",
    "a.\n",
    "\n",
    "   - Firstly, if the $m \\times n$ matrix $\\boldsymbol{A}$ has linearly independent\n",
    "     rows, then the rank of $\\boldsymbol{A}$ is $m$ and the column space of\n",
    "     $\\boldsymbol{A}$ spans $\\mathbb{C}^{m}$, and the nullspace space\n",
    "     of $\\boldsymbol{A}^{H}$ contains the zero vector only.\n",
    "\n",
    "   - Now, consider the nullspace of $\\boldsymbol{A}\\boldsymbol{A}^{H}$:\n",
    "\n",
    "     \\begin{equation}\n",
    "      \\boldsymbol{A}\\boldsymbol{A}^{H} \\boldsymbol{x} = \\boldsymbol{0} \\ \\rightarrow\n",
    "      \\ \\boldsymbol{x}^{H}\\boldsymbol{A} \\boldsymbol{A}^{H} \\boldsymbol{x} = 0 \\ \\rightarrow\n",
    "      \\ (\\boldsymbol{A}^{H} \\boldsymbol{x})^{H}\\boldsymbol{A}^{H} \\boldsymbol{x} = 0\n",
    "     \\end{equation}\n",
    "\n",
    "     The above holds only if $\\boldsymbol{A}^{H} \\boldsymbol{x} = \\boldsymbol{0}$, which\n",
    "     says that $\\boldsymbol{x}$ must come from the nullspace of\n",
    "     $\\boldsymbol{A}^{H}$.  We have already determined that the nullspace of\n",
    "     $\\boldsymbol{A}^{H}$ contains only the zero vector, therefore\n",
    "     $\\boldsymbol{A}\\boldsymbol{A}^{H}$ is full rank (the nullspace of\n",
    "     $\\boldsymbol{A}\\boldsymbol{A}^{H}$ contains the zero vector only) and can be inverted.\n",
    "\n",
    "b. Since $\\boldsymbol{A}\\boldsymbol{A}^{H}$ is square and full rank, it can be inverted,\n",
    "   $$\n",
    "     \\boldsymbol{A} \\boldsymbol{A}^{+} = \\boldsymbol{A} \\boldsymbol{A}^{H}\n",
    "     (\\boldsymbol{A}\\boldsymbol{A}^{H})^{-1} = \\boldsymbol{I}.\n",
    "   $$\n",
    "   hence $\\boldsymbol{A}^{+}$ is a right-inverse of $\\boldsymbol{A}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 14 (stationary iterative methods)\n",
    "\n",
    "Define matrix $\\boldsymbol{A}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2 -1]\n",
      " [-1  2]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[2, -1], [-1, 2]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We split $\\boldsymbol{A}$ such that $\\boldsymbol{A} = \\boldsymbol{N} - \\boldsymbol{P}$. A method will converge if the largest absolute eigenvalue of $\\boldsymbol{N}^{-1}\\boldsymbol{P}$ is less the one.\n",
    "\n",
    "For the Richardson method $\\boldsymbol{N} = \\boldsymbol{I}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0. -2.]\n"
     ]
    }
   ],
   "source": [
    "# Richardson\n",
    "N = np.identity(2)\n",
    "P = N - A\n",
    "M = np.linalg.inv(N).dot(P) \n",
    "print(np.linalg.eigvals(M))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Largest eigenvalue (absolute value) is greater than 1, therefore method will not converge.\n",
    "\n",
    "For the Jacobi method, $\\boldsymbol{N} = \\text{diag}(\\boldsymbol{A})$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.5 -0.5]\n"
     ]
    }
   ],
   "source": [
    "# Jacobi\n",
    "N = np.diag(np.diag(A))\n",
    "P = N - A\n",
    "M = np.linalg.inv(N).dot(P) \n",
    "print(np.linalg.eigvals(M))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Largest eigenvalue (absolute value) is less than 1, therefore method will converge.\n",
    "\n",
    "For Gauss-Seidel, $\\boldsymbol{N}$ is the lower triangular part of $\\boldsymbol{A}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.   0.25]\n"
     ]
    }
   ],
   "source": [
    "# Gauss-Seidel\n",
    "N = np.tril(A)\n",
    "P = N - A\n",
    "M = np.linalg.inv(N).dot(P) \n",
    "print(np.linalg.eigvals(M))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gauss-Seidel will converge because largest eigenvalue is less than one."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 15 (SVD)\n",
    "\n",
    "\\begin{align}\n",
    "  \\boldsymbol{A}^{-1} &= \\left(\\boldsymbol{U} \\boldsymbol{\\Sigma} \\boldsymbol{V}^{H}\\right)^{-1} \\\\\n",
    "  &=\\boldsymbol{V}^{-H} \\boldsymbol{\\Sigma}^{-1} \\boldsymbol{U}^{-1} \\\\\n",
    "  &=\\boldsymbol{V} \\boldsymbol{\\Sigma}^{-1} \\boldsymbol{U}^{H}\n",
    "\\end{align}\n",
    "Non-singular matrix cannot have any zero singular values. In fact, smallest singular values is a measure of the 'distance' to a singular matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 16 (SVD)\n",
    "\n",
    "Define matrix $\\boldsymbol{A}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 3]\n",
      " [2 2]\n",
      " [3 1]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1, 3], [2, 2], [3, 1]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the reduced SVD (recall that NumPy uses $\\boldsymbol{A} = \\boldsymbol{U} \\boldsymbol{\\Sigma} \\boldsymbol{V}$ rather than  $\\boldsymbol{A} = \\boldsymbol{U} \\boldsymbol{\\Sigma} \\boldsymbol{V}^{T}$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[4.89897949 2.        ]\n",
      "[[-5.77350269e-01  7.07106781e-01]\n",
      " [-5.77350269e-01  8.98662938e-17]\n",
      " [-5.77350269e-01 -7.07106781e-01]]\n",
      "[[-0.70710678 -0.70710678]\n",
      " [-0.70710678  0.70710678]]\n"
     ]
    }
   ],
   "source": [
    "U, s, V = np.linalg.svd(A, full_matrices=False)\n",
    "print(s)\n",
    "print(U)\n",
    "print(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute pseudoinverse by creating $\\boldsymbol{\\Sigma}^{+} =  \\boldsymbol{\\Sigma}_{1}^{-1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.16666667  0.08333333  0.33333333]\n",
      " [ 0.33333333  0.08333333 -0.16666667]]\n"
     ]
    }
   ],
   "source": [
    "# Pseudoinverse\n",
    "Sigma_p = np.diag(1.0/s)\n",
    "Ap = (V.T).dot(Sigma_p.dot(U.T))\n",
    "print(Ap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute $\\boldsymbol{A}^{+}\\boldsymbol{A}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.00000000e+00 -3.88578059e-16]\n",
      " [ 2.49800181e-16  1.00000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "# Check that A^{+}A = I\n",
    "print(Ap.dot(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which is the identity. Compute now $\\boldsymbol{A}\\boldsymbol{A}^{+}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.83333333  0.33333333 -0.16666667]\n",
      " [ 0.33333333  0.33333333  0.33333333]\n",
      " [-0.16666667  0.33333333  0.83333333]]\n"
     ]
    }
   ],
   "source": [
    "# Check that AA^{+} \\ne I\n",
    "print(A.dot(Ap))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which is clearly not the identity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(e) Recall that from $\\boldsymbol{A}^{T} \\boldsymbol{A} \\hat{\\boldsymbol{x}} = \\boldsymbol{A}^{T} \\boldsymbol{b}$ we have \n",
    "\\begin{align}\n",
    "\\hat{\\boldsymbol{x}} &= (\\boldsymbol{A}^{T} \\boldsymbol{A})^{-1}\\boldsymbol{A}^{T} \\boldsymbol{b} \\\\\n",
    "&= \\boldsymbol{A}^{+} \\boldsymbol{b}\n",
    "\\end{align}\n",
    "Multiplying both sides by $\\boldsymbol{A}$, \n",
    "$$\n",
    "\\boldsymbol{A}\\hat{\\boldsymbol{x}} = \\underbrace{\\boldsymbol{A} \\boldsymbol{A}^{+}}_{\\boldsymbol{P}} \\boldsymbol{b}\n",
    "$$\n",
    "where $\\boldsymbol{P}$ is the projection matrix from an earlier questions. Recall that $\\boldsymbol{P}$ projects a vector into the column space of $\\boldsymbol{A}$.\n",
    "Since $\\boldsymbol{P}$ is a projector, it does nothing if $\\boldsymbol{b}$ is already in column space. Therefore any vector $\\boldsymbol{b}$ in column space of $\\boldsymbol{A}$ is a solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 17 (pseudo inverse)\n",
    "\n",
    "Define matrix $\\boldsymbol{A}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 3 0]\n",
      " [2 0 0]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[0, 3, 0], [2, 0, 0]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute SVD:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 0.]\n",
      " [0. 1.]]\n",
      "[3. 2.]\n",
      "[[0. 1. 0.]\n",
      " [1. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "U, s, V = np.linalg.svd(A, full_matrices=False)\n",
    "print(U)\n",
    "print(s)\n",
    "print(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute $\\boldsymbol{\\Sigma}^{+} =  \\boldsymbol{\\Sigma}_{1}^{-1}$, and then $\\boldsymbol{A}^{+}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.         0.5       ]\n",
      " [0.33333333 0.        ]\n",
      " [0.         0.        ]]\n"
     ]
    }
   ],
   "source": [
    "# \\Sigma^{+}\n",
    "Sigma_p = np.diag(1.0/s)\n",
    "\n",
    "# A^{+}\n",
    "Ap = (V.T).dot(Sigma_p).dot(U)\n",
    "print(Ap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute $\\boldsymbol{A}^{+}\\boldsymbol{A}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 0. 0.]\n",
      " [0. 1. 0.]\n",
      " [0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "ApA = Ap.dot(A)\n",
    "print(ApA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this matrix, any $\\boldsymbol{x} = [x_{1} \\ \\ x_{2} \\ \\ 0]$ (which is from the row space of $\\boldsymbol{A}$) satisfies $\\boldsymbol{A}^{+} \\boldsymbol{A} \\boldsymbol{x} = \\boldsymbol{x}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 18 (rank deficient least squares)\n",
    "\n",
    "Define matrix and RHS vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 0 0]\n",
      " [1 0 0]\n",
      " [1 1 1]]\n",
      "[0 2 2]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1, 0, 0], [1, 0, 0], [1, 1, 1]])\n",
    "print(A)\n",
    "\n",
    "b = np.array([0, 2, 2])\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the SVD and print singular values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2. 1. 0.]\n"
     ]
    }
   ],
   "source": [
    "# Compute SVD\n",
    "U, s, V = np.linalg.svd(A)\n",
    "print(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is one zero singular value, so we need to 'trim' the last column from $\\boldsymbol{U}$ and the last row from $\\boldsymbol{U}$, and compute $\\boldsymbol{\\sigma}^{+}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.40824829 -0.57735027]\n",
      " [-0.40824829 -0.57735027]\n",
      " [-0.81649658  0.57735027]]\n",
      "[[-0.81649658 -0.40824829 -0.40824829]\n",
      " [-0.57735027  0.57735027  0.57735027]]\n",
      "[[0.5 0. ]\n",
      " [0.  1. ]]\n"
     ]
    }
   ],
   "source": [
    "# Create view of U with last columns \n",
    "U1 = U[:, :2]\n",
    "print(U1)\n",
    "V1 = V[:2,:]\n",
    "print(V1)\n",
    "\n",
    "# Create Sigma^{+}\n",
    "S1 = np.diag(1.0/s[:-1])\n",
    "print(S1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve problem $\\boldsymbol{x} = \\boldsymbol{V}_{1} \\boldsymbol{\\Sigma}_{1}^{-1} \\boldsymbol{U}_{1}^{T}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.  0.5 0.5]\n"
     ]
    }
   ],
   "source": [
    "x = np.transpose(V1).dot(S1.dot(U1.T).dot(b))\n",
    "print(x)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}