LA4DS_ch13.ipynb

{
 "nbformat": 4,
 "nbformat_minor": 0,
 "metadata": {
  "colab": {
   "name": "LA4DS_ch13.ipynb",
   "provenance": [],
   "collapsed_sections": [],
   "authorship_tag": "ABX9TyMI+fm1cZnhp9T4aQtJBwtP"
  },
  "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 13"
   ],
   "metadata": {
    "id": "SbGFWGzkd44U"
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "D9pnXvaOVrxh"
   },
   "outputs": [],
   "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"
   ]
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "ECeT7dTbxV2D"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Geometry of eigenvectors"
   ],
   "metadata": {
    "id": "LJ7QkMH5xVu1"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# in 2D of course, for visualization\n",
    "\n",
    "# the matrix\n",
    "M = np.array([ [-1,1],\n",
    "               [-1,2] ])\n",
    "\n",
    "# its eigenvalues and eigenvectors\n",
    "eigenvalues,eigenvectors = np.linalg.eig(M)\n",
    "print(eigenvalues)\n",
    "\n",
    "# some random vectors\n",
    "notEigenvectors = np.random.randn(2,2)\n",
    "\n",
    "# multipy to create new vectors\n",
    "Mv = M @ eigenvectors\n",
    "Mw = M @ notEigenvectors\n",
    "\n",
    "\n",
    "\n",
    "## and now plot\n",
    "_,axs = plt.subplots(1,2,figsize=(10,6))\n",
    "\n",
    "# the two eigenvectors\n",
    "axs[0].plot([0,eigenvectors[0,0]],[0,eigenvectors[1,0]],'k',linewidth=2,label='$v_1$')\n",
    "axs[0].plot([0,Mv[0,0]],[0,Mv[1,0]],'k--',linewidth=2,label='$Mv_1$')\n",
    "\n",
    "axs[0].plot([0,eigenvectors[0,1]],[0,eigenvectors[1,1]],'r',linewidth=2,label='$v_2$')\n",
    "axs[0].plot([0,Mv[0,1]],[0,Mv[1,1]],'r--',linewidth=2,label='$Mv_2$')\n",
    "\n",
    "# the two non-eigenvectors\n",
    "axs[1].plot([0,notEigenvectors[0,0]],[0,notEigenvectors[1,0]],'k',linewidth=2,label='$w_1$')\n",
    "axs[1].plot([0,Mw[0,0]],[0,Mw[1,0]],'k--',linewidth=2,label='$Mw_1$')\n",
    "\n",
    "axs[1].plot([0,notEigenvectors[0,1]],[0,notEigenvectors[1,1]],'r',linewidth=2,label='$w_2$')\n",
    "axs[1].plot([0,Mw[0,1]],[0,Mw[1,1]],'r--',linewidth=2,label='$Mw_2$')\n",
    "\n",
    "\n",
    "# adjust the graphs a bit\n",
    "for i in range(2):\n",
    "  axs[i].axis('square')\n",
    "  axs[i].set_xlim([-1.5,1.5])\n",
    "  axs[i].set_ylim([-1.5,1.5])\n",
    "  axs[i].grid()\n",
    "  axs[i].legend()\n",
    "\n",
    "plt.savefig('Figure_13_01.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "mtyEVYsTxY3O"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "qPBU7z76TVsN"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Finding eigenvalues"
   ],
   "metadata": {
    "id": "-rmMzuh3cQAh"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "matrix = np.array([\n",
    "             [1,2],\n",
    "             [3,4]\n",
    "             ])\n",
    "\n",
    "# get the eigenvalues\n",
    "evals = np.linalg.eig(matrix)[0]\n",
    "evals"
   ],
   "metadata": {
    "id": "wlA3lzFaYVS5"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# Finding eigenvectors\n",
    "\n",
    "evals,evecs = np.linalg.eig(matrix)\n",
    "print(evals), print(' ')\n",
    "print(evecs)"
   ],
   "metadata": {
    "id": "7VxNv7eEZx-y"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "eWzvQjVSxNmd"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Finding eigenvectors"
   ],
   "metadata": {
    "id": "4ql4otn-H7cU"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# same matrix as above\n",
    "evals,evecs = np.linalg.eig(matrix)\n",
    "\n",
    "print('List of eigenvalues:')\n",
    "print(evals)\n",
    "\n",
    "print(f'\\nMatrix of eigenvectors (in the columns!):')\n",
    "print(evecs)"
   ],
   "metadata": {
    "id": "3zygvzknH9og"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "KRJSm0YnxOHQ"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Diagonalizing a matrix"
   ],
   "metadata": {
    "id": "qXHhZO4Ip4nw"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# using variables created above\n",
    "D = np.diag(evals)\n",
    "D"
   ],
   "metadata": {
    "id": "4P3qgVdBp52j"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# confirm the matrix eigenvalue equation:\n",
    "LHS = matrix @ evecs\n",
    "RHS = evecs @ D\n",
    "\n",
    "\n",
    "# print out the two sides of the equation\n",
    "print('Left-hand side:')\n",
    "print(LHS)\n",
    "\n",
    "print(f'\\nRight-hand side:')\n",
    "print(RHS)"
   ],
   "metadata": {
    "id": "qYMtFhCKp55c"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "SPnIOStPp58n"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Special properties of symmetric matrices"
   ],
   "metadata": {
    "id": "0LFLh86Sp5_e"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# just some random matrix\n",
    "A = np.random.randint(-3,4,(3,3))\n",
    "\n",
    "# and make it symmetric\n",
    "A = A.T@A\n",
    "\n",
    "# its eigendecomposition\n",
    "L,V = np.linalg.eig(A)\n",
    "\n",
    "# all pairwise dot products\n",
    "print( np.dot(V[:,0],V[:,1]) )\n",
    "print( np.dot(V[:,0],V[:,2]) )\n",
    "print( np.dot(V[:,1],V[:,2]) )"
   ],
   "metadata": {
    "id": "i69Glp4VGnGh"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# show that V'V=I\n",
    "np.round( V.T@V ,10) # rounded for visibility (precision errors...)"
   ],
   "metadata": {
    "id": "8Uc7WZqhHEig"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# real-valued matrix with complex-valued eigenvalues\n",
    "\n",
    "# a matrix\n",
    "A = np.array([[-3, -3, 0],\n",
    "              [ 3, -2, 3],\n",
    "              [ 0,  1, 2]])\n",
    "\n",
    "\n",
    "# btw, random matrices often have complex eigenvalues (though this is not guaranteed):\n",
    "#A = np.random.randint(-3,4,(3,3))\n",
    "\n",
    "# its eigendecomposition\n",
    "L,V = np.linalg.eig(A)\n",
    "L.reshape(-1,1) # print as column vector"
   ],
   "metadata": {
    "id": "x38WoIfRp6Cl"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# repeat for symmetric matrices\n",
    "\n",
    "# a matrix\n",
    "A = np.array([[-3, -3, 0],\n",
    "              [-3, -2, 1],\n",
    "              [ 0,  1, 2]])\n",
    "\n",
    "\n",
    "# you can also demonstrate this with random symmetric matrices\n",
    "#A = np.random.randint(-3,4,(3,3))\n",
    "#A = A.T@A\n",
    "\n",
    "# its eigendecomposition\n",
    "L,V = np.linalg.eig(A)\n",
    "L.reshape(-1,1) # print as column vector"
   ],
   "metadata": {
    "id": "eu5mTjE_6FT7"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "vh_5iTTuoiP3"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Eigendecomposition of singular matrices"
   ],
   "metadata": {
    "id": "q34eID_-p6IU"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# a singular matrix\n",
    "A = np.array([[1,4,7],\n",
    "              [2,5,8],\n",
    "              [3,6,9]])\n",
    "\n",
    "# its eigendecomposition\n",
    "L,V = np.linalg.eig(A)\n",
    "\n",
    "\n",
    "# print its rank...\n",
    "print( f'Rank = {np.linalg.matrix_rank(A)}\\n' )\n",
    "\n",
    "# ... and its eigendecomposition\n",
    "print('Eigenvalues: ')\n",
    "print(L.round(2)), print(' ')\n",
    "\n",
    "print('Eigenvectors:')\n",
    "print(V.round(2))"
   ],
   "metadata": {
    "id": "1Lgb_vxnLxgU"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# FYI, random singular matrix\n",
    "M = np.random.randn(5,3) @ np.random.randn(3,5)\n",
    "M = M.T@M # make it symmetric for real-valued eigenvalues\n",
    "\n",
    "# print its eigenvalues (rounded and columnized for clarity)\n",
    "np.linalg.eig(M)[0].reshape(-1,1).round(3)"
   ],
   "metadata": {
    "id": "HMCtlmP4p6LU"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "dbtzVjDup6Fd"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Quadratic form"
   ],
   "metadata": {
    "id": "hCc5qgUyoiDd"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# a matrix with only positive quad.form values\n",
    "A = np.array([ [2,4],[0,3] ])\n",
    "print('Eigenvalues: ')\n",
    "print(np.linalg.eig(A)[0])\n",
    "\n",
    "# print the quadratic form for some random vectors\n",
    "x,y = np.random.randn(2)\n",
    "print(f'\\nSome random quadratic form result:')\n",
    "A[0,0]*x**2 + (A[1,0]+A[0,1])*x*y + A[1,1]*y**2"
   ],
   "metadata": {
    "id": "v_kZW4_Hokdj"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# a matrix with both positive and negative quad.form values\n",
    "A = np.array([ [-9,4],[3,9] ])\n",
    "print('Eigenvalues: ')\n",
    "print(np.linalg.eig(A)[0])\n",
    "\n",
    "# print the quadratic form for some random vectors\n",
    "x,y = np.random.randn(2)\n",
    "print(f'\\nSome random quadratic form result:')\n",
    "A[0,0]*x**2 + (A[1,0]+A[0,1])*x*y + A[1,1]*y**2"
   ],
   "metadata": {
    "id": "v4UGlyP2oiGa"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "TBz7VbpmoiJn"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Generalized eigendecomposition"
   ],
   "metadata": {
    "id": "MPj51LQbp6Qe"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "n = 4\n",
    "\n",
    "# create symmetric matrices\n",
    "A = np.random.randn(n,n)\n",
    "A = A.T@A\n",
    "\n",
    "# impose a correlation between the two matrices (this improves numerical stability of the simultaneousl diagonalization)\n",
    "B = np.random.randn(n,n)\n",
    "B = B.T@B + A/10\n",
    "\n",
    "\n",
    "# using scipy\n",
    "from scipy.linalg import eigh\n",
    "evals,evecs = eigh(A,B)\n",
    "evals"
   ],
   "metadata": {
    "id": "omA9EZhWTHoK"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "4oxx8v10TH5c"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 1"
   ],
   "metadata": {
    "id": "dGnfrHRDTH70"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# create the matrix\n",
    "A = np.random.randn(5,5)\n",
    "A = A.T@A\n",
    "\n",
    "# compute its inverse\n",
    "Ai = np.linalg.inv(A)\n",
    "\n",
    "# eigenvalues of A and Ai\n",
    "eigvals_A  = np.linalg.eig(A)[0]\n",
    "eigvals_Ai = np.linalg.eig(Ai)[0]\n",
    "\n",
    "# compare them (hint: sorting helps!)\n",
    "print('Eigenvalues of A:')\n",
    "print(np.sort(eigvals_A))\n",
    "\n",
    "print(' ')\n",
    "print('Eigenvalues of inv(A):')\n",
    "print(np.sort(eigvals_Ai))\n",
    "\n",
    "print(' ')\n",
    "print('Reciprocal of evals of inv(A):')\n",
    "print(np.sort(1/eigvals_Ai))"
   ],
   "metadata": {
    "id": "ai_jjY8cHVxT"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "BOFlJNKUHV3s"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 2"
   ],
   "metadata": {
    "id": "dsHNYTEzHV73"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# the matrix\n",
    "M = np.array([ [-1,1],\n",
    "               [-1,2] ])\n",
    "\n",
    "# its eigenvalues and eigenvectors\n",
    "eigenvalues,eigenvectors = np.linalg.eig(M)\n",
    "\n",
    "# some random vectors\n",
    "notEigenvectors = np.random.randn(2,2)\n",
    "\n",
    "# multipy to create new vectors\n",
    "Mv = M @ eigenvectors\n",
    "Mw = M @ notEigenvectors\n",
    "\n",
    "\n",
    "\n",
    "## and now plot\n",
    "_,axs = plt.subplots(1,2,figsize=(10,6))\n",
    "\n",
    "# the two eigenvectors\n",
    "axs[0].plot([0,eigenvectors[0,0]],[0,eigenvectors[0,1]],'k',linewidth=2,label='$v_1$')\n",
    "axs[0].plot([0,Mv[0,0]],[0,Mv[0,1]],'k--',linewidth=2,label='$Mv_1$')\n",
    "\n",
    "axs[0].plot([0,eigenvectors[1,0]],[0,eigenvectors[1,1]],'r',linewidth=2,label='$v_2$')\n",
    "axs[0].plot([0,Mv[1,0]],[0,Mv[1,1]],'r--',linewidth=2,label='$Mv_2$')\n",
    "\n",
    "# the two non-eigenvectors\n",
    "axs[1].plot([0,notEigenvectors[0,0]],[0,notEigenvectors[0,1]],'k',linewidth=2,label='$w_1$')\n",
    "axs[1].plot([0,Mw[0,0]],[0,Mw[0,1]],'k--',linewidth=2,label='$Mw_1$')\n",
    "\n",
    "axs[1].plot([0,notEigenvectors[1,0]],[0,notEigenvectors[1,1]],'r',linewidth=2,label='$w_2$')\n",
    "axs[1].plot([0,Mw[1,0]],[0,Mw[1,1]],'r--',linewidth=2,label='$Mw_2$')\n",
    "\n",
    "\n",
    "# adjust the graphs a bit\n",
    "for i in range(2):\n",
    "  axs[i].axis('square')\n",
    "  axs[i].set_xlim([-1.5,1.5])\n",
    "  axs[i].set_ylim([-1.5,1.5])\n",
    "  axs[i].grid()\n",
    "  axs[i].legend()\n",
    "\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "xdLZT8KJIXEg"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "pB_minX89kTK"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 3"
   ],
   "metadata": {
    "id": "g39ycc-kIXKK"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# instructions don't specify matrix size; I'll use n=5\n",
    "N = 5\n",
    "\n",
    "# to store the reconstruction accuracies\n",
    "reconAcc = np.zeros(4)\n",
    "\n",
    "\n",
    "# Create a symmetric random-integers matrix\n",
    "A = np.random.randn(N,N)\n",
    "A = np.round( A.T+A )\n",
    "\n",
    "# diagonalize the matrix\n",
    "d,V  = np.linalg.eig(A)\n",
    "D    = np.diag(d)\n",
    "\n",
    "# demonstrate reconstruction accuracy\n",
    "# remember that inv(V)=V.T!\n",
    "Arecon = V @ D @ V.T\n",
    "print(np.round( A-Arecon ,4))\n",
    "\n",
    "reconAcc[0] = np.sqrt(np.sum( (A-Arecon)**2 ))\n",
    "print(f'\\nFrobenius distance: {reconAcc[0]}')"
   ],
   "metadata": {
    "id": "ZG5ftkYkIZGK"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# create D-tilde\n",
    "Dtild = np.diag( d[np.random.permutation(N)] )\n",
    "\n",
    "# test reconstruction accuracy\n",
    "Arecon = V @ Dtild @ V.T\n",
    "print(np.round( A-Arecon ,4))\n",
    "\n",
    "reconAcc[1] = np.sqrt(np.sum( (A-Arecon)**2 ))\n",
    "print(f'\\nFrobenius distance: {reconAcc[1]}')"
   ],
   "metadata": {
    "id": "FjJ2vtY5IXM2"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "### swap only the two largest eigenvalues\n",
    "evals_sort_idx = np.argsort(d) # note: default is to sort \n",
    "i = evals_sort_idx[np.r_[np.arange(N-2),N-1,N-2]][::-1]\n",
    "\n",
    "# create D-tilde\n",
    "Dtild = np.diag( d[i] )\n",
    "\n",
    "# test reconstruction accuracy\n",
    "Arecon = V @ Dtild @ V.T\n",
    "print(np.round( A-Arecon ,4))\n",
    "\n",
    "reconAcc[2] = np.sqrt(np.sum( (A-Arecon)**2 ))\n",
    "print(f'\\nFrobenius distance: {reconAcc[2]}')"
   ],
   "metadata": {
    "id": "KgUEYyUgHUyX"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "### swap only the two smallest eigenvalues\n",
    "evals_sort_idx = np.argsort(d) # note: default is to sort \n",
    "i = evals_sort_idx[np.r_[1,0,np.arange(2,N)]][::-1]\n",
    "\n",
    "# create D-tilde\n",
    "Dtild = np.diag( d[i] )\n",
    "\n",
    "# test reconstruction accuracy\n",
    "Arecon = V @ Dtild @ V.T\n",
    "print(np.round( A-Arecon ,4))\n",
    "\n",
    "reconAcc[3] = np.sqrt(np.sum( (A-Arecon)**2 ))\n",
    "print(f'\\nFrobenius distance: {reconAcc[3]}')"
   ],
   "metadata": {
    "id": "xtpJH7O7HU2Y"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# now for the plot\n",
    "\n",
    "plt.figure(figsize=(8,6))\n",
    "\n",
    "plt.bar(range(4),reconAcc)\n",
    "plt.xticks(range(4),labels=['None','All','Largest two','Smallest two'])\n",
    "plt.ylabel('Frobenius distance to original matrix')\n",
    "plt.xlabel('Type of eigenvalue swapping')\n",
    "plt.title('Reconstruction accuracy')\n",
    "\n",
    "plt.savefig('Figure_13_03.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "EvQzptn8HU5a"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "arnF5i6MHU8B"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 4"
   ],
   "metadata": {
    "id": "_Mv0Bml7HU-h"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "nIter = 123\n",
    "matsize = 42\n",
    "evals = np.zeros((nIter,matsize),dtype=complex)\n",
    "\n",
    "# create the matrices and get their scaled eigenvalues\n",
    "for i in range(nIter):\n",
    "  A = np.random.randn(matsize,matsize)\n",
    "  evals[i,:] = np.linalg.eig(A)[0] / np.sqrt(matsize)\n",
    "\n",
    "\n",
    "\n",
    "# and show in a plot\n",
    "plt.figure(figsize=(6,6))\n",
    "\n",
    "plt.plot(np.real(evals),np.imag(evals),'ko',markerfacecolor='white')\n",
    "plt.xlim([-1.5,1.5])\n",
    "plt.ylim([-1.5,1.5])\n",
    "plt.xlabel('Real')\n",
    "plt.ylabel('Imag')\n",
    "plt.savefig('Figure_13_04.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "KusdvR3mWS5M"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "MhPPCaM_WS8W"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 5"
   ],
   "metadata": {
    "id": "yOPG26-rWS_T"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# get the null_space function from scipy\n",
    "from scipy.linalg import null_space\n",
    "\n",
    "\n",
    "# Create a symmetric matrix\n",
    "N = 3\n",
    "A = np.random.randn(N,N)\n",
    "A = A@A.T\n",
    "\n",
    "# eigendecompose\n",
    "evals,evecs = np.linalg.eig(A)\n",
    "\n",
    "# compare the eigenvectors with N(A-lI)\n",
    "for i in range(N):\n",
    "\n",
    "  # get the null space vector of the shifted matrix\n",
    "  nullV = null_space( A-evals[i]*np.eye(N) )\n",
    "\n",
    "  # check for a match with the eigenvector via correlation (normalizes for magnitudes)\n",
    "  r = np.corrcoef(nullV.T,evecs[[i],:])[0,1]\n",
    "\n",
    "  # and print (abs(r))\n",
    "  print(f'Correlation between N(A-lI) and evec {i}: {np.abs(r):.2f}')"
   ],
   "metadata": {
    "id": "sr1SZfl8Y23O"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "e9YG3YKVY26G"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 6"
   ],
   "metadata": {
    "id": "IGtulTl8Y28u"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Create the Lambda matrix with positive values\n",
    "Lambda = np.diag( np.random.rand(4)*5 )\n",
    "\n",
    "# create Q\n",
    "Q,_ = np.linalg.qr( np.random.randn(4,4) )\n",
    "\n",
    "# reconstruct to a matrix\n",
    "A = Q @ Lambda @ Q.T\n",
    "\n",
    "# the matrix minus its transpose should be zeros (within precision error)\n",
    "np.round( A-A.T ,5)"
   ],
   "metadata": {
    "id": "RQa0NuvMgmGt"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# check eigenvalues against Lambda (sorting is helpful!)\n",
    "print(np.sort(np.diag(Lambda)))\n",
    "print(np.sort(np.linalg.eig(A)[0]))"
   ],
   "metadata": {
    "id": "584KCUoKgmJq"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "FHrMopF7gmMf"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 7"
   ],
   "metadata": {
    "id": "X9kkEUjxgmPW"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Refer back to the code for Chapter 12, exercise 4."
   ],
   "metadata": {
    "id": "haBkmQRbh_3K"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "RYYonQqHh_6D"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 8"
   ],
   "metadata": {
    "id": "IRZQEjaCgmSC"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# correlation matrix\n",
    "R = np.array([[ 1,.2,.9],\n",
    "              [.2, 1,.3],\n",
    "              [.9,.3, 1] ])\n",
    "\n",
    "# eigendecomposition\n",
    "d,V = np.linalg.eig(R)\n",
    "D = np.diag(d)\n",
    "\n",
    "# create new data with imposed correlation\n",
    "X = V @ np.sqrt(D) @ np.random.randn(3,10000)\n",
    "\n",
    "np.corrcoef(X)"
   ],
   "metadata": {
    "id": "Jwu_qHBNiBUK"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "cO_9R_c0rHIW"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 9"
   ],
   "metadata": {
    "id": "9rvnt6dWrG71"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# now whiten\n",
    "Y = X.T @ V @ np.linalg.inv(np.sqrt(D))\n",
    "\n",
    "# and check the correlations\n",
    "np.round( np.corrcoef(Y.T) ,3)"
   ],
   "metadata": {
    "id": "CmjHJE1giBWf"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "zqa9Xu7xY2_R"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 10"
   ],
   "metadata": {
    "id": "AeztM214Y3Bv"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# two symmetric matrices and GED\n",
    "n = 5\n",
    "A = np.random.randn(n,n)\n",
    "A = A.T@A\n",
    "B = np.random.randn(n,n)\n",
    "B = B.T@B\n",
    "\n",
    "evals,evecs = eigh( A,B )\n",
    "\n",
    "# eigenvectors times their transposes\n",
    "VV  = evecs.T @ evecs\n",
    "VBV = evecs.T @ B @ evecs\n",
    "\n",
    "\n",
    "# show in an image\n",
    "_,axs = plt.subplots(1,2,figsize=(10,6))\n",
    "\n",
    "axs[0].imshow(VV,cmap='gray')\n",
    "axs[0].set_title('$\\mathbf{V}^T\\mathbf{V}$')\n",
    "\n",
    "axs[1].imshow(VBV,cmap='gray')\n",
    "axs[1].set_title('$\\mathbf{V}^T\\mathbf{B}\\mathbf{V}$')\n",
    "\n",
    "plt.savefig('Figure_13_05.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "w_6xkc0DtS1o"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "YIpE7M0ItS7j"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 11"
   ],
   "metadata": {
    "id": "ckTJ1gcytS-q"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# create the matrix\n",
    "A = np.random.randint(-14,15,(4,4))\n",
    "\n",
    "\n",
    "# diagonalize\n",
    "d,V = np.linalg.eig(A)\n",
    "V   = V*np.pi\n",
    "D   = np.diag(d)\n",
    "Vi  = np.linalg.inv(V)\n",
    "\n",
    "\n",
    "# test for accurate reconstruction\n",
    "print('Reconstructed minus original:')\n",
    "print( np.round(V@D@Vi - A,3) )\n",
    "print(' ')\n",
    "\n",
    "# norms of the eigenvectors\n",
    "for i in range(A.shape[0]):\n",
    "  norm = np.sqrt(np.sum(V[:,i]*np.conj(V[:,i])))\n",
    "  print(f'Eigenvector {i} has norm {norm}')\n",
    "\n",
    "\n",
    "# Discussion: Scaling V doesn't matter because that scalar is normalized out in the matrix inverse."
   ],
   "metadata": {
    "id": "5ydkl4Vo99Ix"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "## repeat for a symmetric matrix using V' instead of inv(V)\n",
    "# create the matrix\n",
    "A = np.random.randint(-14,15,(4,4))\n",
    "A = A.T@A\n",
    "\n",
    "\n",
    "# diagonalize\n",
    "d,V = np.linalg.eig(A)\n",
    "V = V*np.pi\n",
    "D = np.diag(d)\n",
    "Vi = V.T\n",
    "\n",
    "\n",
    "# test for accurate reconstruction\n",
    "print('Reconstructed minus original:')\n",
    "print( np.round(V@D@Vi - A,3) )\n",
    "print(' ')\n",
    "\n",
    "# norms of the eigenvectors\n",
    "for i in range(A.shape[0]):\n",
    "  norm = np.sqrt(np.sum(V[:,i]*np.conj(V[:,i])))\n",
    "  print(f'Eigenvector {i} has norm {norm}')\n",
    "\n",
    "\n",
    "# Discussion: Scaling V *does* matter because V is not explicitly inverted!"
   ],
   "metadata": {
    "id": "lF7868llDGeq"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# "
   ],
   "metadata": {
    "id": "wJmjk6fQtTBu"
   },
   "execution_count": null,
   "outputs": []
  }
 ]
}