LA4DS_ch08.ipynb

{
 "nbformat": 4,
 "nbformat_minor": 0,
 "metadata": {
  "colab": {
   "name": "LA4DS_ch08.ipynb",
   "provenance": [],
   "collapsed_sections": [],
   "authorship_tag": "ABX9TyPOR7ncTzmPaCV3gCoDByVT"
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Practical Linear Algebra for Data Science\n",
    "## Mike X Cohen (sincxpress.com)\n",
    "### https://www.oreilly.com/library/view/practical-linear-algebra/9781098120603/\n",
    "\n",
    "#### Code for chapter 8"
   ],
   "metadata": {
    "id": "SbGFWGzkd44U"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# NOTE: these lines define global figure properties used for publication.\n",
    "import matplotlib_inline.backend_inline\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats('svg') # display figures in vector format\n",
    "plt.rcParams.update({'font.size':14}) # set global font size"
   ],
   "metadata": {
    "id": "KqvHgWXHlZjY"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "CgEFEzu47xLN"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# The matrix inverse"
   ],
   "metadata": {
    "id": "0vCbZcd07yYd"
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "KyivubDplVat"
   },
   "outputs": [],
   "source": [
    "# a matrix\n",
    "A = np.array([ [1,4],[2,7] ])\n",
    "\n",
    "# its inverse\n",
    "Ainv = np.linalg.inv(A)\n",
    "\n",
    "# confirm that it produces the identity matrix\n",
    "A@Ainv"
   ]
  },
  {
   "cell_type": "code",
   "source": [
    "# The matrices visualized\n",
    "\n",
    "fig,axs = plt.subplots(1,3,figsize=(10,6))\n",
    "\n",
    "\n",
    "# the matrix\n",
    "axs[0].imshow(A,cmap='gray')\n",
    "axs[0].set_title('The matrix')\n",
    "for (j,i),num in np.ndenumerate(A):\n",
    "  axs[0].text(i,j,num,color=[.8,.8,.8],ha='center',va='center',fontsize=28)\n",
    "\n",
    "# its inverse\n",
    "axs[1].imshow(Ainv,cmap='gray')\n",
    "axs[1].set_title('Its inverse')\n",
    "for (j,i),num in np.ndenumerate(Ainv):\n",
    "  axs[1].text(i,j,num,color=[.8,.8,.8],ha='center',va='center',fontsize=28)\n",
    "\n",
    "# their product\n",
    "AAi = A@Ainv\n",
    "axs[2].imshow(AAi,cmap='gray')\n",
    "axs[2].set_title('Their product')\n",
    "for (j,i),num in np.ndenumerate(AAi):\n",
    "  axs[2].text(i,j,num,color=[.8,.8,.8],ha='center',va='center',fontsize=28)\n",
    "\n",
    "\n",
    "# common properties\n",
    "for i in range(3):\n",
    "  axs[i].set_xticks([])\n",
    "  axs[i].set_yticks([])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('Figure_08_01.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "o8efOyor732_"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "KLWQjOGC75ux"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# the full inverse is two-sided\n",
    "print( A@Ainv ), print(' ')\n",
    "print( Ainv@A )"
   ],
   "metadata": {
    "id": "qbollo2olhGu"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# reminder to use the correct operator:\n",
    "A*Ainv # Hadamard multiplication!"
   ],
   "metadata": {
    "id": "Z3elO0gElovq"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# try again with a singular matrix\n",
    "A = np.array([ [1,4],[2,8] ])\n",
    "\n",
    "# its inverse\n",
    "Ainv = np.linalg.inv(A)\n",
    "\n",
    "# does it produce the identity matrix?\n",
    "A@Ainv"
   ],
   "metadata": {
    "id": "2eIRa7zSobRs"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "064zh1iKeiHm"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Inverse of a diagonal matrix"
   ],
   "metadata": {
    "id": "mJRYL0ceeiSw"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "D = np.diag( np.arange(1,6) )\n",
    "Dinv = np.linalg.inv(D)\n",
    "\n",
    "print('The diagonal matrix:')\n",
    "print(D), print(' ')\n",
    "\n",
    "print('Its inverse:')\n",
    "print(Dinv), print(' ')\n",
    "\n",
    "print('Their product:')\n",
    "print(D@Dinv)"
   ],
   "metadata": {
    "id": "NDkTHmqpekWu"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "wO0dB11weiVe"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# The left-inverse"
   ],
   "metadata": {
    "id": "ho1--WmbeiYC"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# making an invertible square matrix from a tall full column-rank matrix\n",
    "\n",
    "# here's a tall matrix.\n",
    "T = np.random.randint(-10,11,size=(40,4))\n",
    "\n",
    "# confirm that it has its maximum possible rank (full column-rank)\n",
    "print( f'This matrix has rank={np.linalg.matrix_rank(T)}\\n\\n' )\n",
    "\n",
    "# next, create a square full-rank matrix\n",
    "TtT = T.T@T\n",
    "\n",
    "# check whether it has an inverse\n",
    "TtT_inv = np.linalg.inv(TtT)\n",
    "print( np.round(TtT_inv@TtT,4) )"
   ],
   "metadata": {
    "id": "vTutSKdlqm1a"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# finish creating the left-inverse\n",
    "\n",
    "# our left-inverse\n",
    "L = TtT_inv @ T.T\n",
    "\n",
    "# confirm that it works\n",
    "print( np.round( L@T,6 ) ), print(' ')\n",
    "\n",
    "# but it's one-sided!\n",
    "print( np.round( T@L,6 ) )\n"
   ],
   "metadata": {
    "id": "r2LBAGIrZY4Y"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# visualize! of course :)\n",
    "\n",
    "fig,axs = plt.subplots(2,2,figsize=(10,10))\n",
    "\n",
    "axs[0,0].imshow(T,cmap='gray')\n",
    "axs[0,0].set_title('Tall matrix')\n",
    "\n",
    "axs[0,1].imshow(L,cmap='gray')\n",
    "axs[0,1].set_title('Left inverse')\n",
    "\n",
    "axs[1,0].imshow(L@T,cmap='gray')\n",
    "axs[1,0].set_title('L@T')\n",
    "\n",
    "axs[1,1].imshow(T@L,cmap='gray')\n",
    "axs[1,1].set_title('T@L')\n",
    "\n",
    "for a in axs.flatten():\n",
    "  a.set_xticks([])\n",
    "  a.set_yticks([])\n",
    "  \n",
    "plt.tight_layout()\n",
    "plt.savefig('Figure_08_04.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "5aOfDBtjXlgD"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "zHfYnGAAQRnc"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "pU8VXS0ZQRsw"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# MP pseudoinverse"
   ],
   "metadata": {
    "id": "yqiDEgEwXl2I"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# The same singular matrix as before\n",
    "A = np.array([ [1,4],[2,8] ])\n",
    "\n",
    "# its inverse\n",
    "Apinv = np.linalg.pinv(A)\n",
    "print(Apinv*85), print(' ')\n",
    "\n",
    "# does it produce the identity matrix?\n",
    "A@Apinv"
   ],
   "metadata": {
    "id": "yRmmDFGxXqfS"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# an exmple with random numbers\n",
    "A = np.random.randn(7,5) @ np.random.randn(5,7)\n",
    "print(f'The rank of this matrix is {np.linalg.matrix_rank(A)}.\\n')\n",
    "\n",
    "Apinv = np.linalg.pinv(A)\n",
    "plt.imshow(A@Apinv)\n",
    "plt.title('The matrix times its pinv')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "ZaoQSjckXl5j"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "ry7dnwRcXl8s"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 1"
   ],
   "metadata": {
    "id": "rkRbbrzyNnFi"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "n = 5\n",
    "\n",
    "# the matrix\n",
    "A = np.random.randn(n,n)\n",
    "\n",
    "# its inverse, and its inverse's inverse\n",
    "Ai  = np.linalg.inv(A)\n",
    "Aii = np.linalg.inv(Ai)\n",
    "\n",
    "# equal the original matrix within tolerance\n",
    "np.round( A-Aii ,10)"
   ],
   "metadata": {
    "id": "neUW7OpCNpJT"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "bc2RFoe_NnKY"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 2"
   ],
   "metadata": {
    "id": "WLso1YT5lFR2"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# create matrix\n",
    "m = 4\n",
    "A = np.random.randn(m,m)\n",
    "\n",
    "# initialize\n",
    "M = np.zeros((m,m)) # minors matrix\n",
    "G = np.zeros((m,m)) # grid matrix\n",
    "\n",
    "# compute minors matrix\n",
    "for i in range(m):\n",
    "  for j in range(m):\n",
    "    \n",
    "    # select rows and cols\n",
    "    rows = [True]*m\n",
    "    rows[i] = False\n",
    "    \n",
    "    cols = [True]*m\n",
    "    cols[j] = False\n",
    "    \n",
    "    # compute the minors\n",
    "    M[i,j]=np.linalg.det(A[rows,:][:,cols])\n",
    "    \n",
    "    # compute Grid\n",
    "    G[i,j] = (-1)**(i+j)\n",
    "\n",
    "        \n",
    "# compute cofactors matrix\n",
    "C = M * G\n",
    "\n",
    "# compute adjugate matrix\n",
    "Ainv = C.T / np.linalg.det(A)\n",
    "\n",
    "# 'regular' inverse function\n",
    "AinvI = np.linalg.inv(A)\n",
    "\n",
    "# compare against inv()\n",
    "np.round( AinvI-Ainv ,8)"
   ],
   "metadata": {
    "id": "xygHMRZNldLE"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# plot them\n",
    "\n",
    "fig,axs = plt.subplots(2,3,figsize=(14,7))\n",
    "\n",
    "axs[0,0].imshow(M,cmap='gray')\n",
    "axs[0,0].set_title('Minors matrix')\n",
    "\n",
    "axs[0,1].imshow(G,cmap='gray')\n",
    "axs[0,1].set_title('Grid matrix')\n",
    "\n",
    "axs[0,2].imshow(C,cmap='gray')\n",
    "axs[0,2].set_title('Cofactors matrix')\n",
    "\n",
    "axs[1,0].imshow(Ainv,cmap='gray')\n",
    "axs[1,0].set_title('Adjugate (inverse)')\n",
    "\n",
    "axs[1,1].imshow(AinvI,cmap='gray')\n",
    "axs[1,1].set_title('np.linalg.inv')\n",
    "\n",
    "axs[1,2].imshow(A@Ainv,cmap='gray')\n",
    "axs[1,2].set_title('A@Ainv')\n",
    "\n",
    "for a in axs.flatten():\n",
    "  a.set_xticks([])\n",
    "  a.set_yticks([])\n",
    "\n",
    "\n",
    "plt.savefig('Figure_08_03.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "HJYFfC0HljWm"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "r_Hs9rf9lFU3"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 4"
   ],
   "metadata": {
    "id": "DxQiunFXVMUp"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Start from the code for the left-inverse, and swap as necessary.\n",
    "\n",
    "# here's a wide matrix.\n",
    "W = np.random.randint(-10,11,size=(4,40))\n",
    "\n",
    "# confirm that it has its maximum possible rank (full row-rank)\n",
    "print( f'This matrix has rank={np.linalg.matrix_rank(W)}\\n\\n' )\n",
    "\n",
    "# next, create a square full-rank matrix\n",
    "WWt = W@W.T\n",
    "\n",
    "# check whether it has an inverse\n",
    "WWt_inv = np.linalg.inv(WWt)\n",
    "print( np.round(WWt_inv@WWt,4) )"
   ],
   "metadata": {
    "id": "agRugmk2VMU2"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# finish creating the right-inverse\n",
    "\n",
    "# our right-inverse\n",
    "R = W.T @ WWt_inv\n",
    "\n",
    "# confirm that it works\n",
    "print( np.round( W@R,6 ) ), print(' ')\n",
    "\n",
    "# but it's one-sided!\n",
    "print( np.round( R@W,6 ) )\n"
   ],
   "metadata": {
    "id": "hG7w2nkzVMU2"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# visualize! of course :)\n",
    "\n",
    "fig,axs = plt.subplots(2,2,figsize=(10,10))\n",
    "\n",
    "axs[0,0].imshow(W,cmap='gray')\n",
    "axs[0,0].set_title('Wide matrix')\n",
    "\n",
    "axs[0,1].imshow(R,cmap='gray')\n",
    "axs[0,1].set_title('Right inverse')\n",
    "\n",
    "axs[1,0].imshow(R@W,cmap='gray')\n",
    "axs[1,0].set_title('R@W')\n",
    "\n",
    "axs[1,1].imshow(W@R,cmap='gray')\n",
    "axs[1,1].set_title('W@R')\n",
    "\n",
    "for a in axs.flatten():\n",
    "  a.set_xticks([])\n",
    "  a.set_yticks([])\n",
    "  \n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "HgvTTg7RVMU3"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "9SsY7JsYWPaN"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 5"
   ],
   "metadata": {
    "id": "xuM_RQvvlFXf"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Full inverse case\n",
    "M = 4\n",
    "\n",
    "A = np.random.randn(M,M)\n",
    "\n",
    "Ainv  = np.linalg.inv(A)\n",
    "Apinv = np.linalg.pinv(A)\n",
    "\n",
    "np.round( Ainv-Apinv,10 )"
   ],
   "metadata": {
    "id": "Gs3OIFGGPhTy"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# left inverse case\n",
    "M,N = 14,4\n",
    "\n",
    "A = np.random.randn(M,N)\n",
    "\n",
    "ALeft = np.linalg.inv(A.T@A) @ A.T\n",
    "Apinv = np.linalg.pinv(A)\n",
    "\n",
    "np.round( ALeft-Apinv,10 )"
   ],
   "metadata": {
    "id": "lwZ4x2pFPfKS"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# right inverse case\n",
    "M,N = 4,14\n",
    "\n",
    "A = np.random.randn(M,N)\n",
    "\n",
    "ARight = A.T @ np.linalg.inv(A@A.T)\n",
    "Apinv  = np.linalg.pinv(A)\n",
    "\n",
    "np.round( ARight-Apinv,10 )"
   ],
   "metadata": {
    "id": "zBlbKvnxPfM8"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "CiZfkIrZPfP5"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 6"
   ],
   "metadata": {
    "id": "RcGvtRAyPfTE"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# create the matrices\n",
    "N = 4\n",
    "A = np.random.randn(N,N)\n",
    "B = np.random.randn(N,N)\n",
    "\n",
    "# compute the three specified options\n",
    "op1 = np.linalg.inv(A@B)\n",
    "op2 = np.linalg.inv(A) @ np.linalg.inv(B)\n",
    "op3 = np.linalg.inv(B) @ np.linalg.inv(A)\n",
    "\n",
    "# compute distances\n",
    "dist12 = np.sqrt(np.sum( (op1-op2)**2 ))\n",
    "dist13 = np.sqrt(np.sum( (op1-op3)**2 ))\n",
    "\n",
    "# print results!\n",
    "print(f'Distance between (AB)^-1 and (A^-1)(B^-1) is {dist12:.8f}')\n",
    "print(f'Distance between (AB)^-1 and (B^-1)(A^-1) is {dist13:.8f}')"
   ],
   "metadata": {
    "id": "cRbW0CT7Y3MB"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "0k8WxCOkp0Pa"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 7"
   ],
   "metadata": {
    "id": "Brhf9aasp0Tc"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# create the matrices\n",
    "M,N = 14,4\n",
    "T = np.random.randn(M,N)\n",
    "\n",
    "# compute the three specified options\n",
    "op1 = np.linalg.inv(T.T@T)\n",
    "op2 = np.linalg.inv(T) @ np.linalg.inv(T.T)\n",
    "\n",
    "# The answer is No, it doesn't work, because a tall matrix has no inverse."
   ],
   "metadata": {
    "id": "cvdW1BpNao0v"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "7hQ-lzldao3j"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 8"
   ],
   "metadata": {
    "id": "-443e5NZwaXH"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Transformation matrix\n",
    "T = np.array([ \n",
    "              [1,.5],\n",
    "              [0,.5]\n",
    "            ])\n",
    "\n",
    "# Compute its inverse\n",
    "Ti = np.linalg.inv(T)\n",
    "\n",
    "\n",
    "# define the set of points (a circle)\n",
    "theta = np.linspace(0,2*np.pi-2*np.pi/20,20)\n",
    "origPoints = np.vstack( (np.cos(theta),np.sin(theta)) )\n",
    "\n",
    "# apply transformation\n",
    "transformedPoints = T @ origPoints\n",
    "\n",
    "# undo the transformation via the inverse of the transform\n",
    "backTransformed   = Ti @ transformedPoints\n",
    "\n",
    "\n",
    "# plot the points\n",
    "plt.figure(figsize=(6,6))\n",
    "plt.plot(origPoints[0,:],origPoints[1,:],'ko',label='Original')\n",
    "plt.plot(transformedPoints[0,:],transformedPoints[1,:],'s',\n",
    "         color=[.7,.7,.7],label='Transformed')\n",
    "plt.plot(backTransformed[0,:],backTransformed[1,:],'rx',markersize=15,\n",
    "         color=[.7,.7,.7],label='Inverse-transformed')\n",
    "\n",
    "plt.axis('square')\n",
    "plt.xlim([-2,2])\n",
    "plt.ylim([-2,2])\n",
    "plt.legend()\n",
    "plt.savefig('Figure_08_06.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "CYv5RDb5p0bX"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "s-EDdvRNwj_s"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 9"
   ],
   "metadata": {
    "id": "y-rB90yRXl_o"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# a function to create a Hilbert matrix\n",
    "def hilbmat(k):\n",
    "  H = np.zeros((k,k))\n",
    "  for i in range(k):\n",
    "    for j in range(k):\n",
    "\n",
    "      # note: the math formula has denominator: i+j-1\n",
    "      #   but with 0-based indexing, this is: (i+1)+(j+1)-1\n",
    "      #   which can be shortened to: i+j+1\n",
    "\n",
    "      H[i,j] = 1 / (i+j+1)\n",
    "  return H\n",
    "\n",
    "  \n",
    "\n",
    "# The double for-loop above is a direct implementation of the math.\n",
    "# The function below gives the same result but without the loops.\n",
    "def hilbmat(k):\n",
    "  k = np.arange(1,k+1).reshape(1,-1) # reshape to a row vector (instead of a 1D array)\n",
    "  return 1 / (k.T+k-1) # outer product and element-wise division"
   ],
   "metadata": {
    "id": "Ucq_gWr5lE94"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "print( hilbmat(5) ), print(' ')\n",
    "\n",
    "# you can confirm the accuracy of your function against the scipy Hilbert-matrix function:\n",
    "from scipy.linalg import hilbert\n",
    "print( hilbert(5) )"
   ],
   "metadata": {
    "id": "m0_hMHAtljdU"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# create a 5x5 Hilbert matrix and show it, its inverse, and their product\n",
    "H = hilbmat(5)\n",
    "Hi = np.linalg.inv(H)\n",
    "\n",
    "fig,axs = plt.subplots(1,3,figsize=(12,6))\n",
    "h = [0,0,0]\n",
    "\n",
    "# the matrix\n",
    "h[0] = axs[0].imshow(H,cmap='gray')\n",
    "axs[0].set_title('Hilbert matrix')\n",
    "\n",
    "# its inverse\n",
    "h[1] = axs[1].imshow(Hi,cmap='gray')\n",
    "axs[1].set_title('Inv(hilbert)')\n",
    "\n",
    "# their product\n",
    "h[2] = axs[2].imshow(H@Hi,cmap='gray')\n",
    "axs[2].set_title('Their product')\n",
    "\n",
    "\n",
    "for i in range(2):\n",
    "  fig.colorbar(h[i],ax=axs[i],fraction=.045)\n",
    "  axs[i].set_xticks([])\n",
    "  axs[i].set_yticks([])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('Figure_08_05.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "nBhTOHNRlFIa"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "2kRscLQojJOl"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 10"
   ],
   "metadata": {
    "id": "aVTBbgoUjJjR"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "\n",
    "matSizes = np.arange(3,13)\n",
    "\n",
    "identityError = np.zeros((len(matSizes),2))\n",
    "condNumbers   = np.zeros((len(matSizes),2))\n",
    "\n",
    "\n",
    "for i,k in enumerate(matSizes):\n",
    "    \n",
    "    ### for the Hilbert matrix\n",
    "    H   = hilbmat(k)       # the matrix\n",
    "    Hi  = np.linalg.inv(H) # its inverse\n",
    "    HHi = H@Hi             # should be identity\n",
    "    err = HHi - np.eye(k)  # difference from true identity\n",
    "    identityError[i,0] = np.sqrt(np.sum(err**2))  # Euclidean distance\n",
    "    condNumbers[i,0] = np.linalg.cond(H) # condition number\n",
    "    \n",
    "\n",
    "    ### repeat for a random matrix\n",
    "    H = np.random.randn(k,k) # the matrix\n",
    "    Hi  = np.linalg.inv(H)   # its inverse\n",
    "    HHi = H@Hi               # should be identity\n",
    "    err = HHi - np.eye(k)    # difference from true identity\n",
    "    identityError[i,1] = np.sqrt(np.sum(err**2))  # Euclidean distance\n",
    "    condNumbers[i,1] = np.linalg.cond(H) # condition number\n",
    "    \n",
    "\n",
    "\n",
    "# now plot\n",
    "fig,axs = plt.subplots(1,2,figsize=(14,5))\n",
    "\n",
    "## plot the Euclidean distance to the identity matrix\n",
    "h = axs[0].plot(matSizes,np.log(identityError),'s-',markersize=12)\n",
    "h[0].set_color('k') # adjust the individual line colors and shapes\n",
    "h[0].set_marker('o')\n",
    "h[1].set_color('gray')\n",
    "\n",
    "axs[0].legend(['Hilbert','Random'])\n",
    "axs[0].set_xlabel('Matrix size')\n",
    "axs[0].set_ylabel('Log Euclidan distance')\n",
    "axs[0].set_title('Distance to identity matrix')\n",
    "\n",
    "\n",
    "\n",
    "## plot the condition numbers\n",
    "h = axs[1].plot(matSizes,np.log(condNumbers),'s-',markersize=12)\n",
    "h[0].set_color('k') # adjust the individual line colors and shapes\n",
    "h[0].set_marker('o')\n",
    "h[1].set_color('gray')\n",
    "\n",
    "axs[1].legend(['Hilbert','Random'])\n",
    "axs[1].set_xlabel('Matrix size')\n",
    "axs[1].set_ylabel('Log Kappa')\n",
    "axs[1].set_title('Matrix condition number')\n",
    "\n",
    "plt.savefig('Figure_08_07.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "ZtaAeWj2lFLx"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "## interesting to see the \"identity\" matrix\n",
    "H   = hilbmat(k)\n",
    "Hi  = np.linalg.inv(H)\n",
    "HHi = H@Hi \n",
    "\n",
    "plt.imshow(HHi,vmin=0,vmax=.1)\n",
    "plt.title('Should be identity')\n",
    "plt.colorbar();"
   ],
   "metadata": {
    "id": "rSkb1Tscp0XV"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "qemw6B5FwgTm"
   },
   "execution_count": null,
   "outputs": []
  }
 ]
}