LA4DS_ch11.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "LA4DS_ch11.ipynb",
      "provenance": []
    },
    "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 11"
      ],
      "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",
        "# null space\n",
        "from scipy.linalg import null_space\n",
        "\n",
        "import sympy as sym\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": "code",
      "source": [],
      "metadata": {
        "id": "F8nKu5cQjyeq"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "## this code creates figure 2\n",
        "\n",
        "# data\n",
        "x = [ 1,2,3,4,5 ]\n",
        "y = [ 0,3,2,5,5 ]\n",
        "\n",
        "# model\n",
        "X = np.hstack((np.ones((5,1)),np.array(x,ndmin=2).T))\n",
        "yHat = X @ np.linalg.inv(X.T@X) @ X.T @ y\n",
        "\n",
        "# plot the data and predicted values\n",
        "plt.figure(figsize=(6,6))\n",
        "plt.plot(x,y,'ks',markersize=15,label='Observed data')\n",
        "plt.plot(x,yHat,'o-',color=[.6,.6,.6],linewidth=3,markersize=8,label='Predicted data')\n",
        "\n",
        "# plot the residuals (errors)\n",
        "for n,y,yHat in zip(x,y,yHat):\n",
        "  plt.plot([n,n],[y,yHat],'--',color=[.8,.8,.8],zorder=-10)\n",
        "\n",
        "plt.legend()\n",
        "plt.savefig('Figure_11_02.png',dpi=300)\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "NsQQmXVz2lDS"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "mWvZfl5qjyhT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Example in fake data"
      ],
      "metadata": {
        "id": "0c2ueLP3uz5b"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "numcourses = [13,4,12,3,14,13,12,9,11,7,13,11,9,2,5,7,10,0,9,7]\n",
        "happiness  = [70,25,54,21,80,68,84,62,57,40,60,64,45,38,51,52,58,21,75,70]\n",
        "\n",
        "plt.figure(figsize=(6,6))\n",
        "\n",
        "plt.plot(numcourses,happiness,'ks',markersize=15)\n",
        "plt.xlabel('Number of courses taken')\n",
        "plt.ylabel('General life happiness')\n",
        "plt.xlim([-1,15])\n",
        "plt.ylim([0,100])\n",
        "plt.grid()\n",
        "plt.xticks(range(0,15,2))\n",
        "plt.savefig('Figure_11_03.png',dpi=300)\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "QbDG-cjnjc6Q"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Build a statistical model\n",
        "\n",
        "# design matrix as a column vector\n",
        "X = np.array(numcourses,ndmin=2).T\n",
        "print(X.shape)\n",
        "\n",
        "# fit the model using the left-inverse\n",
        "X_leftinv = np.linalg.inv(X.T@X) @ X.T\n",
        "\n",
        "# solve for the coefficients\n",
        "beta = X_leftinv @ happiness\n",
        "beta"
      ],
      "metadata": {
        "id": "7oG9rOjRjagr"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# let's plot it!\n",
        "\n",
        "# predicted data\n",
        "pred_happiness = X@beta\n",
        "\n",
        "\n",
        "plt.figure(figsize=(6,6))\n",
        "\n",
        "# plot the data and predicted values\n",
        "plt.plot(numcourses,happiness,'ks',markersize=15)\n",
        "plt.plot(numcourses,pred_happiness,'o-',color=[.6,.6,.6],linewidth=3,markersize=8)\n",
        "\n",
        "# plot the residuals (errors)\n",
        "for n,y,yHat in zip(numcourses,happiness,pred_happiness):\n",
        "  plt.plot([n,n],[y,yHat],'--',color=[.8,.8,.8],zorder=-10)\n",
        "\n",
        "plt.xlabel('Number of courses taken')\n",
        "plt.ylabel('General life happiness')\n",
        "plt.xlim([-1,15])\n",
        "plt.ylim([0,100])\n",
        "plt.xticks(range(0,15,2))\n",
        "plt.legend(['Real data','Predicted data','Residual'])\n",
        "plt.title(f'SSE = {np.sum((pred_happiness-happiness)**2):.2f}')\n",
        "plt.savefig('Figure_11_04.png',dpi=300)\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "jxh-OgHUjaj1"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Build a statistical model with an intercept\n",
        "\n",
        "# design matrix as a column vector\n",
        "X = np.hstack((np.ones((20,1)),np.array(numcourses,ndmin=2).T))\n",
        "print(X.shape)\n",
        "\n",
        "# fit the model using the left-inverse\n",
        "X_leftinv = np.linalg.inv(X.T@X) @ X.T\n",
        "\n",
        "# solve for the coefficients\n",
        "beta = X_leftinv @ happiness\n",
        "beta"
      ],
      "metadata": {
        "id": "74HCXEvvjam4"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# let's plot it!\n",
        "\n",
        "# predicted data\n",
        "pred_happiness = X@beta\n",
        "\n",
        "\n",
        "plt.figure(figsize=(6,6))\n",
        "\n",
        "# plot the data and predicted values\n",
        "plt.plot(numcourses,happiness,'ks',markersize=15)\n",
        "plt.plot(numcourses,pred_happiness,'o-',color=[.6,.6,.6],linewidth=3,markersize=8)\n",
        "\n",
        "# plot the residuals (errors)\n",
        "for n,y,yHat in zip(numcourses,happiness,pred_happiness):\n",
        "  plt.plot([n,n],[y,yHat],'--',color=[.8,.8,.8],zorder=-10)\n",
        "\n",
        "plt.xlabel('Number of courses taken')\n",
        "plt.ylabel('General life happiness')\n",
        "plt.xlim([-1,15])\n",
        "plt.ylim([0,100])\n",
        "plt.xticks(range(0,15,2))\n",
        "plt.legend(['Real data','Predicted data','Residual'])\n",
        "plt.title(f'SSE = {np.sum((pred_happiness-happiness)**2):.2f}')\n",
        "plt.savefig('Figure_11_05.png',dpi=300)\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "jAjOzukzjapt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "7Cp04-Qw2k6T"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Exercise 1"
      ],
      "metadata": {
        "id": "kwU4W1d62lL5"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# compute residual\n",
        "res = happiness-pred_happiness\n",
        "\n",
        "\n",
        "# should be zero + some error\n",
        "print('Dot product: ' + str(np.dot(pred_happiness,res)) )\n",
        "print('Correlation: ' + str(np.corrcoef(pred_happiness,res)[0,1]))\n",
        "print(' ')\n",
        "\n",
        "\n",
        "# show in a plot\n",
        "plt.figure(figsize=(6,6))\n",
        "plt.plot(res,pred_happiness,'ko',markersize=12)\n",
        "plt.xlabel('Residual error')\n",
        "plt.ylabel('Model-predicted values')\n",
        "plt.title(f'r = {np.corrcoef(pred_happiness,res)[0,1]:.20f}')\n",
        "plt.savefig('Figure_11_06.png',dpi=300)\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "t9QF60gi2m38"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# correlation is smaller because we're dividing by the vector norms, e.g.,\n",
        "np.linalg.norm(res)"
      ],
      "metadata": {
        "id": "jCXHdxB3m30Y"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "s3O1fA8zCUXr"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Exercise 2"
      ],
      "metadata": {
        "id": "o3svhXLXm4bs"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# the residual is orthogonal to the entire column space of the design matrix.\n",
        "\n",
        "# I demonstrated this by showing that the residuals vector is in the left-null space of the design matrix.\n",
        "# I did that by using scipy.linalg.null_space to find the left-null space, augmenting that null-space basis\n",
        "# matrix by the residuals vector, and showing that the null space and augmented null space have the same rank.\n",
        "\n",
        "\n",
        "# compute the null space (via scipy.linalg)\n",
        "nullspace = null_space(X.T)\n",
        "\n",
        "\n",
        "# augment the residuals\n",
        "nullspaceAugment = np.hstack( (nullspace,res.reshape(-1,1)) )\n",
        "\n",
        "\n",
        "# print their ranks\n",
        "print(f'dim(  N(X)    ) = {np.linalg.matrix_rank(nullspace)}')\n",
        "print(f'dim( [N(X)|r] ) = {np.linalg.matrix_rank(nullspaceAugment)}')"
      ],
      "metadata": {
        "id": "TS5yPkwj2m62"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "fI3dJcb42m-B"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Exercise 3"
      ],
      "metadata": {
        "id": "0bPcsugrjast"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "### Uncomment these lines to use random data\n",
        "# random design matrix and data vector\n",
        "# M,N = 20,3\n",
        "# X = np.random.randn(M,N)\n",
        "# happiness = np.random.randn(M,1)\n",
        "\n",
        "\n",
        "# recreate the design matrix and solution via left-inverse\n",
        "X = np.hstack((np.ones((20,1)),np.array(numcourses,ndmin=2).T))\n",
        "beta1 = np.linalg.inv(X.T@X) @ X.T @ happiness\n",
        "\n",
        "\n",
        "# QR decomp\n",
        "Q,R = np.linalg.qr(X)\n",
        "\n",
        "# beta coefficients implemented as translation of the math\n",
        "beta2 = np.linalg.inv(R) @ (Q.T@happiness)\n",
        "\n",
        "# and using back-substitution via RREF\n",
        "# Q'y, but needs to be reshaped into a column vector\n",
        "tmp = (Q.T@happiness).reshape(-1,1)\n",
        "Raug = np.hstack( (R,tmp) ) # augment the matrix\n",
        "Raug_r = sym.Matrix(Raug).rref()[0] # this gets the matrix\n",
        "beta3 = np.array(Raug_r[:,-1]) # convert back to numpy\n",
        "\n",
        "\n",
        "print('Betas from left-inverse: ')\n",
        "print(np.round(beta1,3)), print(' ')\n",
        "\n",
        "print('Betas from QR with inv(R): ')\n",
        "print(np.round(beta2,3)), print(' ')\n",
        "\n",
        "print('Betas from QR with back-substitution: ')\n",
        "print(np.round(np.array(beta3.T).astype(float),3)) # transposed to facilitate visual inspection"
      ],
      "metadata": {
        "id": "LnO4XR0b2nAg"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# show the matrices\n",
        "print('Matrix R:')\n",
        "print(np.round(R,3)) # note that it's upper-triangular (as you know!)\n",
        "\n",
        "print(' ')\n",
        "print(\"Matrix R|Q'y:\")\n",
        "print(np.round(Raug,3))\n",
        "\n",
        "print(' ')\n",
        "print(\"Matrix RREF(R|Q'y):\")\n",
        "print(np.round(np.array(Raug_r).astype(float),3)) # convert to numpy floats"
      ],
      "metadata": {
        "id": "Y_Do_j3d2nDP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "a1D5ATv1Qu9s"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Exercise 4"
      ],
      "metadata": {
        "id": "uThOWkG7JMLz"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# happiness with outliers due to typos (oops!)\n",
        "happiness_oops1 = [170,25,54,21,80,68,84,62,57,40,60,64,45,38,51,52,58,21,75,70]\n",
        "happiness_oops2 = [70,25,54,21,80,68,84,62,57,40,60,64,45,38,51,52,58,21,75,170]\n",
        "\n",
        "\n",
        "# design matrix and its left-inverse (doesn't change with the data)\n",
        "X = np.hstack((np.ones((20,1)),np.array(numcourses,ndmin=2).T))\n",
        "X_leftinv = np.linalg.inv(X.T@X) @ X.T\n",
        "\n",
        "\n",
        "\n",
        "_,axs = plt.subplots(1,3,figsize=(16,5))\n",
        "\n",
        "for axi,y in zip(axs,[happiness,happiness_oops1,happiness_oops2]):\n",
        "\n",
        "  # compute the best-fit parameters\n",
        "  beta = X_leftinv @ y\n",
        "\n",
        "  # predicted data\n",
        "  pred_happiness = X@beta\n",
        "\n",
        "\n",
        "  # plot the data and predicted values\n",
        "  axi.plot(numcourses,y,'ks',markersize=15)\n",
        "  axi.plot(numcourses,pred_happiness,'o-',color=[.6,.6,.6],linewidth=3,markersize=8)\n",
        "\n",
        "  # plot the residuals (errors)\n",
        "  for n,yy,yHat in zip(numcourses,y,pred_happiness):\n",
        "    axi.plot([n,n],[yy,yHat],'--',color=[.8,.8,.8],zorder=-10)\n",
        "\n",
        "  # make the plot look nicer\n",
        "  axi.set(xlabel='Number of courses taken',ylabel='General life happiness',\n",
        "          xlim=[-1,15],ylim=[0,100],xticks=range(0,15,2))\n",
        "  axi.legend(['Real data','Predicted data','Residual'])\n",
        "  axi.set_title(f'SSE = {np.sum((pred_happiness-y)**2):.2f}')\n",
        "\n",
        "\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.savefig('Figure_11_07.png',dpi=300)\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "oh4rHkG6JNqk"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "ztcpxLKZJMQt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Exercise 5"
      ],
      "metadata": {
        "id": "-FGBU3kz2k9P"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# matrix size\n",
        "n = 6\n",
        "\n",
        "# some random \"design matrix\"\n",
        "X = np.random.randn(n,n)\n",
        "\n",
        "# the target matrix (identity)\n",
        "Y = np.eye(n)\n",
        "\n",
        "\n",
        "# find the best-fitting model one column at a time\n",
        "Xinv1 = np.zeros_like(X)\n",
        "\n",
        "for coli in range(n):\n",
        "  Xinv1[:,coli] = np.linalg.inv(X.T@X) @ X.T @ Y[:,coli]\n",
        "\n",
        "\n",
        "\n",
        "# repeat but without a loop\n",
        "Xinv2 = np.linalg.inv(X.T@X) @ X.T @ Y\n",
        "\n",
        "\n",
        "# and the inverse using inv()\n",
        "Xinv3 = np.linalg.inv(X)\n",
        "\n",
        "\n",
        "# visualize\n",
        "_,axs = plt.subplots(1,3,figsize=(10,6))\n",
        "\n",
        "# column-wise least-squares\n",
        "axs[0].imshow( Xinv1@X ,cmap='gray')\n",
        "axs[0].set_title('Via column-wise LS')\n",
        "\n",
        "# matrix-wise least-squares\n",
        "axs[1].imshow( Xinv2@X ,cmap='gray' )\n",
        "axs[1].set_title('Via matrix-wise LS')\n",
        "\n",
        "# inv()\n",
        "axs[2].imshow( Xinv3@X ,cmap='gray' )\n",
        "axs[2].set_title('Via inv()')\n",
        "\n",
        "\n",
        "# don't need the tick marks\n",
        "for a in axs: a.set(xticks=[],yticks=[])\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.savefig('Figure_11_08.png',dpi=300)\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "ecx7hGZe2lAT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# show they are equivalent\n",
        "# Note the relatively large rounding errors when comparing to inv() -- the left-inverse\n",
        "#   least-squares method is not a numerically stable method!\n",
        "\n",
        "\n",
        "print(Xinv1-Xinv2)\n",
        "print(' ')\n",
        "\n",
        "print(Xinv1-Xinv3)\n",
        "print(' ')\n",
        "\n",
        "print(Xinv2-Xinv3)"
      ],
      "metadata": {
        "id": "_XhAQVqtHroi"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "1TCcpX4TPgkN"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}