LA4DS_ch10.ipynb

{
 "nbformat": 4,
 "nbformat_minor": 0,
 "metadata": {
  "colab": {
   "name": "LA4DS_ch10.ipynb",
   "provenance": [],
   "collapsed_sections": [],
   "authorship_tag": "ABX9TyPMViqAGFTfl4rudQdwV256"
  },
  "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 10"
   ],
   "metadata": {
    "id": "SbGFWGzkd44U"
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "pq-JEewsnKHm"
   },
   "outputs": [],
   "source": [
    ""
   ]
  },
  {
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# sympy library for RREF\n",
    "import sympy as sym\n",
    "\n",
    "# scipy for LU\n",
    "import scipy.linalg\n",
    "\n",
    "\n",
    "# used to create non-regular subplots\n",
    "import matplotlib.gridspec as gridspec\n",
    "\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": "GVWZlfT-nThD"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "MvATXzy6gxLF"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Working with matrix equations"
   ],
   "metadata": {
    "id": "gpPMdUOe6uOl"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# generate some matrices\n",
    "A = np.random.randn(4,4)\n",
    "B = np.random.randn(4,4)\n",
    "\n",
    "# solve for X\n",
    "# 1) inv(A)@A@X = inv(A)@B\n",
    "# 2) inv(A)@A@X = B@inv(A)\n",
    "\n",
    "X1 = np.linalg.inv(A) @ B\n",
    "X2 = B @ np.linalg.inv(A)\n",
    "\n",
    "# residual (should be zeros matrix)\n",
    "res1 = A@X1 - B\n",
    "res2 = A@X2 - B\n",
    "\n",
    "# which is correct?\n",
    "print('res1:'), print(' ')\n",
    "print( np.round(res1,10) ), print(' ')\n",
    "\n",
    "print('res2:'), print(' ')\n",
    "print( np.round(res2,10) )\n"
   ],
   "metadata": {
    "id": "aY1G7pax6wuc"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "9MQHoH3L6w18"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# RREF"
   ],
   "metadata": {
    "id": "bfu3STmJ6w5h"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# the augmented matrix\n",
    "M = np.array([ [1,1,4],[-1/2,1,2] ])\n",
    "\n",
    "# converted into a sympy matrix\n",
    "symMat = sym.Matrix(M)\n",
    "print(symMat)\n",
    "\n",
    "# RREF\n",
    "symMat.rref()[0] # just the first output to get the RREF matrix (the second output is the indices of the pivots per row)"
   ],
   "metadata": {
    "id": "A6g4eOq81ejx"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "Stg-pxfE1emn"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# LU"
   ],
   "metadata": {
    "id": "NUvmTdmM1epv"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# simple example with integers\n",
    "\n",
    "# a matrix\n",
    "A = np.array([ [2,2,4], [1,0,3], [2,1,2] ])\n",
    "\n",
    "# its LU decomposition via scipy (please ignore the first output for now)\n",
    "_,L,U = scipy.linalg.lu(A)\n",
    "\n",
    "# print them out\n",
    "print('L: ')\n",
    "print(L), print(' ')\n",
    "\n",
    "print('U: ')\n",
    "print(U), print(' ')\n",
    "\n",
    "print('A - LU: ')\n",
    "print(A - L@U) # should be zeros"
   ],
   "metadata": {
    "id": "0ZotHlkZKrfn"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "# matrix sizes\n",
    "m = 4\n",
    "n = 6\n",
    "\n",
    "A = np.random.randn(m,n)\n",
    "\n",
    "P,L,U = scipy.linalg.lu(A)\n",
    "\n",
    "# show the matrices\n",
    "fig,axs = plt.subplots(1,5,figsize=(13,4))\n",
    "\n",
    "axs[0].imshow(A,vmin=-1,vmax=1)\n",
    "axs[0].set_title('A')\n",
    "\n",
    "axs[1].imshow(np.ones((m,n)),cmap='gray',vmin=-1,vmax=1)\n",
    "axs[1].text(n/2,m/2,'=',ha='center',fontsize=30,fontweight='bold')\n",
    "# axs[1].axis('off')\n",
    "\n",
    "axs[2].imshow(P.T,vmin=-1,vmax=1)\n",
    "axs[2].set_title(r'P$^T$')\n",
    "\n",
    "axs[3].imshow(L,vmin=-1,vmax=1)\n",
    "axs[3].set_title('L')\n",
    "\n",
    "h = axs[4].imshow(U,vmin=-1,vmax=1)\n",
    "axs[4].set_title('U')\n",
    "\n",
    "for a in axs:\n",
    "  a.axis('off')\n",
    "  a.set_xlim([-.5,n-.5])\n",
    "  a.set_ylim([m-.5,-.5])\n",
    "\n",
    "\n",
    "fig.colorbar(h,ax=axs[-1],fraction=.05)\n",
    "plt.tight_layout()\n",
    "plt.savefig('Figure_10_01.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "6gKC5wE1UnkY"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "vo7h7VHCog2p"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 1"
   ],
   "metadata": {
    "id": "74yjyg7_-ZUV"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# Time-test!\n",
    "\n",
    "import time\n",
    "\n",
    "# start the timer\n",
    "tic = time.time()\n",
    "\n",
    "# run the test\n",
    "for i in range(1000):\n",
    "  A = np.random.randn(100,100)\n",
    "  P,L,U = scipy.linalg.lu(A)\n",
    "\n",
    "# stop the timer\n",
    "toc = time.time() - tic\n",
    "toc # print the result in seconds"
   ],
   "metadata": {
    "id": "Un3ZoiQ76w8X"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "_zZV-8T16w_j"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 2"
   ],
   "metadata": {
    "id": "4yx-wGm0nZLl"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# make a reduced-rank random matrix\n",
    "\n",
    "# sizes and rank\n",
    "M = 6\n",
    "N = 8\n",
    "r = 3\n",
    "\n",
    "# create the matrix\n",
    "A = np.random.randn(M,r) @ np.random.randn(r,N)\n",
    "\n",
    "# LU\n",
    "P,L,U = scipy.linalg.lu(A)\n",
    "\n",
    "# and plot\n",
    "_,axs = plt.subplots(1,3,figsize=(12,7))\n",
    "\n",
    "axs[0].imshow(A,vmin=-1,vmax=1,cmap='gray')\n",
    "axs[0].set_title(f'A, rank={np.linalg.matrix_rank(A)}')\n",
    "\n",
    "axs[1].imshow(L,vmin=-1,vmax=1,cmap='gray')\n",
    "axs[1].set_title(f'L, rank={np.linalg.matrix_rank(L)}')\n",
    "\n",
    "axs[2].imshow(U,vmin=-1,vmax=1,cmap='gray')\n",
    "axs[2].set_title(f'U, rank={np.linalg.matrix_rank(U)}')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('Figure_10_02.png',dpi=300)\n",
    "plt.show()"
   ],
   "metadata": {
    "id": "2_g2-nQfuypP"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    "np.round(L,2)"
   ],
   "metadata": {
    "id": "iXHURqD5uzxM"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "swpj20CVPR6N"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 3"
   ],
   "metadata": {
    "id": "noadx007yIVi"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# a matrix and its det\n",
    "M = 6\n",
    "A = np.random.randn(M,M)\n",
    "\n",
    "# LU\n",
    "P,L,U = scipy.linalg.lu(A)\n",
    "\n",
    "# determinant as the product of the diagonals of U\n",
    "detLU = np.prod( np.diag(U) ) * np.linalg.det(P)\n",
    "\n",
    "# check against the det function\n",
    "detNP = np.linalg.det(A)\n",
    "\n",
    "# compare\n",
    "print(detLU,detNP)\n",
    "print(detLU-detNP)"
   ],
   "metadata": {
    "id": "Vu6vRVJkyJdV"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "0poRjiIxUnnT"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 4"
   ],
   "metadata": {
    "id": "fWhk-PGPUnqK"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# matrix sizes\n",
    "m = 4\n",
    "A = np.random.randn(m,m)\n",
    "\n",
    "# LU decomposition\n",
    "P,L,U = scipy.linalg.lu(A)\n",
    "\n",
    "# inverse\n",
    "invViaLU = np.linalg.inv(U) @ np.linalg.inv(L) @ P.T\n",
    "\n",
    "# \"regular\" inverse\n",
    "invViaInv = np.linalg.inv(A)\n",
    "\n",
    "np.round( A@invViaLU ,10)"
   ],
   "metadata": {
    "id": "_mU1rgjEog0B"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "tHUkNFBkuzz-"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Exercise 5"
   ],
   "metadata": {
    "id": "8M49SES1uz2n"
   }
  },
  {
   "cell_type": "code",
   "source": [
    "# The reason is that writing out the equation leads to PtP in the middle, which is the identity matrix. \n",
    "# Conceptually, it means that any row swaps are undone when multiplying by the transpose.\n",
    "\n",
    "# create a matrix\n",
    "A = np.random.randn(4,4)\n",
    "\n",
    "# LUP\n",
    "P,L,U = scipy.linalg.lu(A)\n",
    "\n",
    "# compute AtA via LU\n",
    "AtA_lu = U.T @ L.T @ L @ U\n",
    "\n",
    "# direct computation\n",
    "AtA_direct = A.T @ A\n",
    "\n",
    "# compare to direct computation\n",
    "np.round( AtA_lu - AtA_direct ,10)"
   ],
   "metadata": {
    "id": "2dfDk9yhjycI"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "F8nKu5cQjyeq"
   },
   "execution_count": null,
   "outputs": []
  },
  {
   "cell_type": "code",
   "source": [
    ""
   ],
   "metadata": {
    "id": "0c2ueLP3uz5b"
   },
   "execution_count": null,
   "outputs": []
  }
 ]
}