tutorials/linear_regression_gradient_descent.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fdb3af1b",
   "metadata": {},
   "source": [
    "## Gradient Descent and Linear Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b18e109",
   "metadata": {},
   "source": [
    "**Copyright (c) Meta Platforms, Inc. and affiliates.**\n",
    "\n",
    "This source code is licensed under the MIT license found in the\n",
    "LICENSE file in the root directory of this source tree."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2514d65",
   "metadata": {},
   "source": [
    "We are going to focus on linear regression here, using the least squares methodology and gradient descent. Validating and analyzing linear function models is beyond the scope of this notebook. Refer to the note below to learn more. \n",
    "\n",
    "A **linear regression** fits a linear function through data, often using a least squares method, to make predictions on new data. In this notebook we will learn how to build a linear regression in Kotlin. Along the way we will use gradient descent and DiffKT to perform automatic differentiation.\n",
    "\n",
    "However we are just going to focus on fitting a linear regression here, using the least squares methodology and gradient descent. Validating and analyzing linear function models is beyond the scope of this notebook, and refer to the callout box below to learn more. \n",
    "\n",
    "> Just fitting a line through some points is by no means a guarantee that it will produce good predictions. Concepts like statistical significance, correlation, and prediction intervals can provide tools for further validation from a statistics standpoint. In machine learning, a train/test split is often performed where the linear function is fitted to a majority portion of the data, known as the training data, and measured for prediction performance on the remaining portion, known as the test data.\n",
    "\n",
    "To visualize the linear regression and the least squares method, let's use [Lets-Plot](https://blog.jetbrains.com/kotlin/2020/12/lets-plot-in-kotlin/). It comes with the Kotlin Jupyter kernel which can be installed with the command below, which is also required to run this notebook. \n",
    "\n",
    "```shell\n",
    "conda install kotlin-jupyter-kernel -c jetbrains\n",
    "```\n",
    "\n",
    "Then lets enable Lets-Plot in this notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "025cc660",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "   <div id=\"7oeOnq\"></div>\n",
       "   <script type=\"text/javascript\" data-lets-plot-script=\"library\">\n",
       "       if(!window.letsPlotCallQueue) {\n",
       "           window.letsPlotCallQueue = [];\n",
       "       }; \n",
       "       window.letsPlotCall = function(f) {\n",
       "           window.letsPlotCallQueue.push(f);\n",
       "       };\n",
       "       (function() {\n",
       "           var script = document.createElement(\"script\");\n",
       "           script.type = \"text/javascript\";\n",
       "           script.src = \"https://cdn.jsdelivr.net/gh/JetBrains/lets-plot@v2.3.0/js-package/distr/lets-plot.min.js\";\n",
       "           script.onload = function() {\n",
       "               window.letsPlotCall = function(f) {f();};\n",
       "               window.letsPlotCallQueue.forEach(function(f) {f();});\n",
       "               window.letsPlotCallQueue = [];\n",
       "               \n",
       "               \n",
       "           };\n",
       "           script.onerror = function(event) {\n",
       "               window.letsPlotCall = function(f) {};\n",
       "               window.letsPlotCallQueue = [];\n",
       "               var div = document.createElement(\"div\");\n",
       "               div.style.color = 'darkred';\n",
       "               div.textContent = 'Error loading Lets-Plot JS';\n",
       "               document.getElementById(\"7oeOnq\").appendChild(div);\n",
       "           };\n",
       "           var e = document.getElementById(\"7oeOnq\");\n",
       "           e.appendChild(script);\n",
       "       })();\n",
       "   </script>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%use lets-plot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18bb8b55",
   "metadata": {},
   "source": [
    "Also let's bring in the DiffKT library to aid our gradient descent implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d1b23168",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "@file:DependsOn(\"../kotlin/api/build/libs/api.jar\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be712fe4",
   "metadata": {},
   "source": [
    "Lets first declare some data and plot it here using Kotlin code. We will also count the number of records $ n $ in our data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e531f66f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "   <div id=\"5YXUiu\"></div>\n",
       "   <script type=\"text/javascript\" data-lets-plot-script=\"plot\">\n",
       "       (function() {\n",
       "           var plotSpec={\n",
       "\"mapping\":{\n",
       "\"x\":\"x\",\n",
       "\"y\":\"y\"\n",
       "},\n",
       "\"ggsize\":{\n",
       "\"width\":600.0,\n",
       "\"height\":500.0\n",
       "},\n",
       "\"kind\":\"plot\",\n",
       "\"scales\":[{\n",
       "\"aesthetic\":\"x\",\n",
       "\"breaks\":[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0],\n",
       "\"name\":\"\"\n",
       "},{\n",
       "\"aesthetic\":\"y\",\n",
       "\"breaks\":[4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0,20.0,21.0,22.0,23.0,24.0,25.0,26.0],\n",
       "\"name\":\"\"\n",
       "}],\n",
       "\"layers\":[{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"shape\":1.0,\n",
       "\"size\":3.0,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"point\",\n",
       "\"data\":{\n",
       "}\n",
       "}],\n",
       "\"data\":{\n",
       "\"x\":[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0],\n",
       "\"y\":[5.0,10.0,10.0,15.0,14.0,15.0,19.0,18.0,25.0,23.0]\n",
       "}\n",
       "};\n",
       "           var plotContainer = document.getElementById(\"5YXUiu\");\n",
       "           window.letsPlotCall(function() {{\n",
       "               LetsPlot.buildPlotFromProcessedSpecs(plotSpec, -1, -1, plotContainer);\n",
       "           }});\n",
       "       })();    \n",
       "   </script>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "val xData = listOf(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0)\n",
    "val yData = listOf(5.0, 10.0, 10.0, 15.0, 14.0, 15.0, 19.0, 18.0, 25.0, 23.0)\n",
    "val n = xData.size\n",
    "\n",
    "val xScaleMin = 1\n",
    "val xScaleMax = 10\n",
    "val yScaleMin = 4\n",
    "val yScaleMax = 26\n",
    "\n",
    "// create the scatterplot \n",
    "val p = letsPlot { x=xData; y=yData; } + ggsize(600, 500) + geomPoint(shape=1, size=3) +\n",
    "    scaleXContinuous(\"\", (xScaleMin..xScaleMax).toList()) + \n",
    "    scaleYContinuous(\"\", (yScaleMin..yScaleMax).toList())\n",
    "\n",
    "p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47dfffde",
   "metadata": {},
   "source": [
    "Just looking at this plot we can see there is a visible linear correlation to this data. But how do we fit a linear function to it? \n",
    "\n",
    "First let's draw a line through the points. We will establish how we fit this line later but let's draw a line with a slope of $ 1.93939 $ and an intercept of $ 4.73333 $. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "801f6aa0",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "   <div id=\"QwlIm8\"></div>\n",
       "   <script type=\"text/javascript\" data-lets-plot-script=\"plot\">\n",
       "       (function() {\n",
       "           var plotSpec={\n",
       "\"mapping\":{\n",
       "\"x\":\"x\",\n",
       "\"y\":\"y\"\n",
       "},\n",
       "\"ggsize\":{\n",
       "\"width\":600.0,\n",
       "\"height\":500.0\n",
       "},\n",
       "\"kind\":\"plot\",\n",
       "\"scales\":[{\n",
       "\"aesthetic\":\"x\",\n",
       "\"breaks\":[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0],\n",
       "\"name\":\"\"\n",
       "},{\n",
       "\"aesthetic\":\"y\",\n",
       "\"breaks\":[4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0,20.0,21.0,22.0,23.0,24.0,25.0,26.0],\n",
       "\"name\":\"\"\n",
       "}],\n",
       "\"layers\":[{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"shape\":1.0,\n",
       "\"size\":3.0,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"point\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"color\":\"red\",\n",
       "\"intercept\":4.73333,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"abline\",\n",
       "\"slope\":1.93939,\n",
       "\"data\":{\n",
       "}\n",
       "}],\n",
       "\"data\":{\n",
       "\"x\":[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0],\n",
       "\"y\":[5.0,10.0,10.0,15.0,14.0,15.0,19.0,18.0,25.0,23.0]\n",
       "}\n",
       "};\n",
       "           var plotContainer = document.getElementById(\"QwlIm8\");\n",
       "           window.letsPlotCall(function() {{\n",
       "               LetsPlot.buildPlotFromProcessedSpecs(plotSpec, -1, -1, plotContainer);\n",
       "           }});\n",
       "       })();    \n",
       "   </script>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fun testLine(x: Double) = 1.93939 * x + 4.73333\n",
    "\n",
    "p + geomABLine(slope=1.93939, intercept=4.73333, color=\"red\") "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e66b9401",
   "metadata": {},
   "source": [
    "Note that the line is not fit perfectly through the points, simply because the points do not exist on a straight line. There is a linear pattern but noise has caused them to be above or below this line by some margin. Note that the differences between the line and the points are called **residuals**. Let's plot the residuals as shown below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "456b5367",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "   <div id=\"vTtuOe\"></div>\n",
       "   <script type=\"text/javascript\" data-lets-plot-script=\"plot\">\n",
       "       (function() {\n",
       "           var plotSpec={\n",
       "\"mapping\":{\n",
       "\"x\":\"x\",\n",
       "\"y\":\"y\"\n",
       "},\n",
       "\"ggsize\":{\n",
       "\"width\":600.0,\n",
       "\"height\":500.0\n",
       "},\n",
       "\"kind\":\"plot\",\n",
       "\"scales\":[{\n",
       "\"aesthetic\":\"x\",\n",
       "\"breaks\":[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0],\n",
       "\"name\":\"\"\n",
       "},{\n",
       "\"aesthetic\":\"y\",\n",
       "\"breaks\":[4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0,20.0,21.0,22.0,23.0,24.0,25.0,26.0],\n",
       "\"name\":\"\"\n",
       "}],\n",
       "\"layers\":[{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"shape\":1.0,\n",
       "\"size\":3.0,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"point\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"shape\":1.0,\n",
       "\"size\":3.0,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"point\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"color\":\"red\",\n",
       "\"intercept\":4.73333,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"abline\",\n",
       "\"slope\":1.93939,\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":1.0,\n",
       "\"y\":5.0,\n",
       "\"yend\":6.67272,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":1.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":2.0,\n",
       "\"y\":10.0,\n",
       "\"yend\":8.61211,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":2.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":3.0,\n",
       "\"y\":10.0,\n",
       "\"yend\":10.5515,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":3.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":4.0,\n",
       "\"y\":15.0,\n",
       "\"yend\":12.49089,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":4.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":5.0,\n",
       "\"y\":14.0,\n",
       "\"yend\":14.43028,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":5.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":6.0,\n",
       "\"y\":15.0,\n",
       "\"yend\":16.36967,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":6.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":7.0,\n",
       "\"y\":19.0,\n",
       "\"yend\":18.30906,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":7.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":8.0,\n",
       "\"y\":18.0,\n",
       "\"yend\":20.24845,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":8.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":9.0,\n",
       "\"y\":25.0,\n",
       "\"yend\":22.187839999999998,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":9.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":10.0,\n",
       "\"y\":23.0,\n",
       "\"yend\":24.127229999999997,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":10.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "}],\n",
       "\"data\":{\n",
       "\"x\":[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0],\n",
       "\"y\":[5.0,10.0,10.0,15.0,14.0,15.0,19.0,18.0,25.0,23.0]\n",
       "}\n",
       "};\n",
       "           var plotContainer = document.getElementById(\"vTtuOe\");\n",
       "           window.letsPlotCall(function() {{\n",
       "               LetsPlot.buildPlotFromProcessedSpecs(plotSpec, -1, -1, plotContainer);\n",
       "           }});\n",
       "       })();    \n",
       "   </script>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var residualGraph = p + geomPoint(shape=1, size=3) + \n",
    "    geomABLine(slope=1.93939, intercept=4.73333, color=\"red\")  \n",
    "\n",
    "// generate residuals \n",
    "for ((x, actualY) in xData.zip(yData)) { \n",
    "    residualGraph += geomSegment(x=x, y=actualY, xend=x, yend=testLine(x), size=1)\n",
    "}\n",
    "\n",
    "residualGraph"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc66202d",
   "metadata": {},
   "source": [
    "Numerically speaking, we can calculate the residuals by taking the difference between the actual $ y $ values and the predicted $ \\hat{y} $ values as shown below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e94af15e",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Actual Y: 5.0, Predicted Y: 6.67272, Residual: -1.67272\n",
      "Actual Y: 10.0, Predicted Y: 8.61211, Residual: 1.3878900000000005\n",
      "Actual Y: 10.0, Predicted Y: 10.5515, Residual: -0.5515000000000008\n",
      "Actual Y: 15.0, Predicted Y: 12.49089, Residual: 2.5091099999999997\n",
      "Actual Y: 14.0, Predicted Y: 14.43028, Residual: -0.4302799999999998\n",
      "Actual Y: 15.0, Predicted Y: 16.36967, Residual: -1.3696699999999993\n",
      "Actual Y: 19.0, Predicted Y: 18.30906, Residual: 0.6909400000000012\n",
      "Actual Y: 18.0, Predicted Y: 20.24845, Residual: -2.2484499999999983\n",
      "Actual Y: 25.0, Predicted Y: 22.187839999999998, Residual: 2.812160000000002\n",
      "Actual Y: 23.0, Predicted Y: 24.127229999999997, Residual: -1.1272299999999973\n"
     ]
    }
   ],
   "source": [
    "for ((x,actualY) in xData.zip(yData)) { \n",
    "    val predictedY = testLine(x.toDouble())\n",
    "    val residual = actualY - predictedY\n",
    "    println(\"Actual Y: $actualY, Predicted Y: $predictedY, Residual: $residual\")\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25b72068",
   "metadata": {},
   "source": [
    "How do we measure these residuals in total to see how well the line is fit to the points? We do not necessarily want to sum up the residuals, as the negatives will cancel out the positives. We do not want to sum the absolute values either, because absolute values are mathematically difficult to work with. Those sharp corners at the absolute value are not Calculus-friendly, an attribute needed for gradient descent.\n",
    "\n",
    "Instead we square each of the residuals and sum them, giving us the total error. In mathematical notation, each prediction $ \\hat{y}_i $ is subtracted from the actual value $ y_i $, and then squared. This is what we call our **loss function**.\n",
    "\n",
    "$ E = \\sum_{i=0}^{n}(y_i - \\hat{y}_i)^2 $\n",
    "\n",
    "We can visualize this literally by taking each residual and overlaying a square with a matching side length. The sum of all those square areas will be our **loss**, or total error measuring how far off our points are from the line in aggregate. Notice that by squaring the residuals, we penalize larger residuals by making them even larger ($ e.g. 2^2 = 4 $ but $9^2 = 81 $). The larger the number, the more amplified it will be. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fb665680",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "   <div id=\"v8YMav\"></div>\n",
       "   <script type=\"text/javascript\" data-lets-plot-script=\"plot\">\n",
       "       (function() {\n",
       "           var plotSpec={\n",
       "\"mapping\":{\n",
       "\"x\":\"x\",\n",
       "\"y\":\"y\"\n",
       "},\n",
       "\"ggsize\":{\n",
       "\"width\":600.0,\n",
       "\"height\":500.0\n",
       "},\n",
       "\"kind\":\"plot\",\n",
       "\"scales\":[{\n",
       "\"aesthetic\":\"x\",\n",
       "\"breaks\":[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0],\n",
       "\"name\":\"\"\n",
       "},{\n",
       "\"aesthetic\":\"y\",\n",
       "\"breaks\":[4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0,20.0,21.0,22.0,23.0,24.0,25.0,26.0],\n",
       "\"name\":\"\"\n",
       "}],\n",
       "\"layers\":[{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"shape\":1.0,\n",
       "\"size\":3.0,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"point\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"shape\":1.0,\n",
       "\"size\":3.0,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"point\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"color\":\"red\",\n",
       "\"intercept\":4.73333,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"abline\",\n",
       "\"slope\":1.93939,\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":1.0,\n",
       "\"y\":5.0,\n",
       "\"yend\":6.67272,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":1.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":2.0,\n",
       "\"y\":10.0,\n",
       "\"yend\":8.61211,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":2.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":3.0,\n",
       "\"y\":10.0,\n",
       "\"yend\":10.5515,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":3.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":4.0,\n",
       "\"y\":15.0,\n",
       "\"yend\":12.49089,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":4.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":5.0,\n",
       "\"y\":14.0,\n",
       "\"yend\":14.43028,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":5.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":6.0,\n",
       "\"y\":15.0,\n",
       "\"yend\":16.36967,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":6.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":7.0,\n",
       "\"y\":19.0,\n",
       "\"yend\":18.30906,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":7.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":8.0,\n",
       "\"y\":18.0,\n",
       "\"yend\":20.24845,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":8.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":9.0,\n",
       "\"y\":25.0,\n",
       "\"yend\":22.187839999999998,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":9.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"size\":1.0,\n",
       "\"x\":10.0,\n",
       "\"y\":23.0,\n",
       "\"yend\":24.127229999999997,\n",
       "\"position\":\"identity\",\n",
       "\"xend\":10.0,\n",
       "\"geom\":\"segment\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"shape\":1.0,\n",
       "\"size\":3.0,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"point\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"color\":\"red\",\n",
       "\"intercept\":4.73333,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"abline\",\n",
       "\"slope\":1.93939,\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"ymin\":5.0,\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"xmin\":1.0,\n",
       "\"ymax\":6.67272,\n",
       "\"xmax\":1.6433538461538462,\n",
       "\"color\":\"orange\",\n",
       "\"alpha\":0.3,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"rect\",\n",
       "\"fill\":\"orange\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"ymin\":8.61211,\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"xmin\":2.0,\n",
       "\"ymax\":10.0,\n",
       "\"xmax\":2.5338038461538464,\n",
       "\"color\":\"orange\",\n",
       "\"alpha\":0.3,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"rect\",\n",
       "\"fill\":\"orange\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"ymin\":10.0,\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"xmin\":3.0,\n",
       "\"ymax\":10.5515,\n",
       "\"xmax\":3.2121153846153847,\n",
       "\"color\":\"orange\",\n",
       "\"alpha\":0.3,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"rect\",\n",
       "\"fill\":\"orange\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"ymin\":12.49089,\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"xmin\":4.0,\n",
       "\"ymax\":15.0,\n",
       "\"xmax\":4.965042307692308,\n",
       "\"color\":\"orange\",\n",
       "\"alpha\":0.3,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"rect\",\n",
       "\"fill\":\"orange\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"ymin\":14.0,\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"xmin\":5.0,\n",
       "\"ymax\":14.43028,\n",
       "\"xmax\":5.165492307692308,\n",
       "\"color\":\"orange\",\n",
       "\"alpha\":0.3,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"rect\",\n",
       "\"fill\":\"orange\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"ymin\":15.0,\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"xmin\":6.0,\n",
       "\"ymax\":16.36967,\n",
       "\"xmax\":6.526796153846154,\n",
       "\"color\":\"orange\",\n",
       "\"alpha\":0.3,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"rect\",\n",
       "\"fill\":\"orange\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"ymin\":18.30906,\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"xmin\":7.0,\n",
       "\"ymax\":19.0,\n",
       "\"xmax\":7.2657461538461545,\n",
       "\"color\":\"orange\",\n",
       "\"alpha\":0.3,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"rect\",\n",
       "\"fill\":\"orange\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"ymin\":18.0,\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"xmin\":8.0,\n",
       "\"ymax\":20.24845,\n",
       "\"xmax\":8.864788461538462,\n",
       "\"color\":\"orange\",\n",
       "\"alpha\":0.3,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"rect\",\n",
       "\"fill\":\"orange\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"ymin\":22.187839999999998,\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"xmin\":9.0,\n",
       "\"ymax\":25.0,\n",
       "\"xmax\":10.081600000000002,\n",
       "\"color\":\"orange\",\n",
       "\"alpha\":0.3,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"rect\",\n",
       "\"fill\":\"orange\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"ymin\":23.0,\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"xmin\":10.0,\n",
       "\"ymax\":24.127229999999997,\n",
       "\"xmax\":10.433549999999999,\n",
       "\"color\":\"orange\",\n",
       "\"alpha\":0.3,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"rect\",\n",
       "\"fill\":\"orange\",\n",
       "\"data\":{\n",
       "}\n",
       "}],\n",
       "\"data\":{\n",
       "\"x\":[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0],\n",
       "\"y\":[5.0,10.0,10.0,15.0,14.0,15.0,19.0,18.0,25.0,23.0]\n",
       "}\n",
       "};\n",
       "           var plotContainer = document.getElementById(\"v8YMav\");\n",
       "           window.letsPlotCall(function() {{\n",
       "               LetsPlot.buildPlotFromProcessedSpecs(plotSpec, -1, -1, plotContainer);\n",
       "           }});\n",
       "       })();    \n",
       "   </script>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "// Overlay squares onto residuals to visualize sum of squares\n",
    "var squaresGraph = residualGraph + geomPoint(shape = 1, size = 3) + \n",
    "    geomABLine(slope=1.93939, intercept = 4.73333, color=\"red\")\n",
    "\n",
    "for ((x,actualY) in xData.zip(yData)) {\n",
    "    \n",
    "    val predictedY = testLine(x) \n",
    "    val residual = actualY - predictedY\n",
    "    squaresGraph += geomRect(xmin=x, \n",
    "                             ymin=if (residual < 0) actualY else predictedY, \n",
    "                             xmax=x+abs(residual) * (xScaleMax.toDouble() / yScaleMax), // account for axis scaling\n",
    "                             ymax=if (residual < 0) predictedY else actualY, \n",
    "                             color=\"orange\", \n",
    "                             fill=\"orange\",\n",
    "                             alpha=.3, \n",
    "                            )\n",
    "}\n",
    "\n",
    "squaresGraph"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "744e96ba",
   "metadata": {},
   "source": [
    "Now that we have an objective function, how do we minimize it? What is the line that produces the smallest total area of these squares as possible? One way we can find this is by means of **gradient descent**, which optimizes the parameters by following their derivatives to the minimum loss.  \n",
    "\n",
    "Given we have the $ x $ and $ y $ values provided in the data, we are left to solve the slope $ m $ and intercept $ b $ in the function $ y = mx + b $. \n",
    "\n",
    "$$\n",
    "y = mx + b\n",
    "$$\n",
    "\n",
    "By taking the gradients of our loss function (the sum of squared errors), we can find the gradients for the slope and intercept. We then adjust the parameters by subtracting a fraction of the gradient (known as the **learning rate**) until a gradient of 0 is reached for all parameters. When all gradients are 0, that means the loss has been minimized. When we plot the loss for a simple linear regression with one slope and intercept, we can see it has a convex surface. We use gradient descent to find the slope of 0 for all gradients, putting us at the lowest point of the plane.  \n",
    "\n",
    "![](resources/rWkCTQwALn.png)\n",
    "\n",
    "\n",
    "We need to get to the lowest point on this surface to minimize loss. Thankfully we can use DiffKT to perform this automatic differentiation for us and guide us down the gradients to find the optimal linear regression parameters. Let's first start by declaring two scalar tensors holding our slope and intercept coefficients for our line. We will initialize them at $ 0 $, and gradient descent will optimize them. Let's also re-declare our `xData` and `yData` as tensors and express the sum of squares operation below in the declared `f()` function. The function accepts the slope and intercepts as parameters and returns the sum of squares. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6f466e6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import org.diffkt.*\n",
    "\n",
    "val xData = tensorOf(1f, 2f, 3f, 4f, 5f, 6f, 7f, 8f, 9f, 10f)\n",
    "val yData = tensorOf(5f, 10f, 10f, 15f, 14f, 15f, 19f, 18f, 25f, 23f)\n",
    "\n",
    "// Start m and b at 0\n",
    "var slope: DScalar = FloatScalar.ZERO\n",
    "var intercept: DScalar = FloatScalar.ZERO\n",
    "\n",
    "// calculate sum of squares of the error with given slope and intercept for a line\n",
    "fun f(slope: DScalar, intercept: DScalar): DScalar =\n",
    "    (((xData * slope) + intercept) - yData).pow(2).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b91dd01",
   "metadata": {},
   "source": [
    "To apply gradient descent to this squared error `f()` function, we will adjust the `slope` and `intercept` iteratively. For 1000 iterations, we will pass the `f()` function to the `reverseDerivative()` function to get the gradients for the slope and intercept respectively. Since we have more inputs than outputs (2 inputs, 1 output) for this function, we will get more efficiency using reverse differentiation. In a loop, we will then take the gradients and multiply them by a learning rate `lr` of $ .0025 $. We will then subtract that from the respective coefficients. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "886412f4",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "slope=1.9394095 intercept=4.7332234"
     ]
    }
   ],
   "source": [
    "// The learning rate\n",
    "val lr = .0025F\n",
    "\n",
    "// The number of iterations to perform gradient descent\n",
    "val iterations = 1000\n",
    "\n",
    "// Perform gradient descent\n",
    "for (i in 0..iterations) {\n",
    "\n",
    "    // get gradients for line slope and intercept\n",
    "    val (slopeGradient, interceptGradient) = reverseDerivative(slope, intercept, ::f)\n",
    "\n",
    "    // update m and b by subtracting the (learning rate) * (slope)\n",
    "    slope -= slopeGradient * lr\n",
    "    intercept -= interceptGradient * lr\n",
    "}\n",
    "print(\"slope=$slope intercept=$intercept\") "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6f76467",
   "metadata": {},
   "source": [
    "After the iterations complete, we should get a slope of $ 1.9394095 $ and an intercept of $ 4.7332234 $ which provides the solved slope and intercept respectively. Let's graph this line against our scatterplot. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f5cb2b6c",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "   <div id=\"iQqVbd\"></div>\n",
       "   <script type=\"text/javascript\" data-lets-plot-script=\"plot\">\n",
       "       (function() {\n",
       "           var plotSpec={\n",
       "\"mapping\":{\n",
       "\"x\":\"x\",\n",
       "\"y\":\"y\"\n",
       "},\n",
       "\"ggsize\":{\n",
       "\"width\":600.0,\n",
       "\"height\":500.0\n",
       "},\n",
       "\"kind\":\"plot\",\n",
       "\"scales\":[{\n",
       "\"aesthetic\":\"x\",\n",
       "\"breaks\":[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0],\n",
       "\"name\":\"\"\n",
       "},{\n",
       "\"aesthetic\":\"y\",\n",
       "\"breaks\":[4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0,20.0,21.0,22.0,23.0,24.0,25.0,26.0],\n",
       "\"name\":\"\"\n",
       "}],\n",
       "\"layers\":[{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"shape\":1.0,\n",
       "\"size\":3.0,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"point\",\n",
       "\"data\":{\n",
       "}\n",
       "},{\n",
       "\"mapping\":{\n",
       "},\n",
       "\"stat\":\"identity\",\n",
       "\"color\":\"red\",\n",
       "\"intercept\":4.7332234382629395,\n",
       "\"position\":\"identity\",\n",
       "\"geom\":\"abline\",\n",
       "\"slope\":1.9394094944000244,\n",
       "\"data\":{\n",
       "}\n",
       "}],\n",
       "\"data\":{\n",
       "\"x\":[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0],\n",
       "\"y\":[5.0,10.0,10.0,15.0,14.0,15.0,19.0,18.0,25.0,23.0]\n",
       "}\n",
       "};\n",
       "           var plotContainer = document.getElementById(\"iQqVbd\");\n",
       "           window.letsPlotCall(function() {{\n",
       "               LetsPlot.buildPlotFromProcessedSpecs(plotSpec, -1, -1, plotContainer);\n",
       "           }});\n",
       "       })();    \n",
       "   </script>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p + geomABLine(slope=(slope as FloatScalar).value.toDouble(), \n",
    "               intercept = (intercept as FloatScalar).value.toDouble(), \n",
    "               color=\"red\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "187d503c",
   "metadata": {},
   "source": [
    "Looks pretty good! That line is comfortably fitted along those points. Now we have done linear regression by applying gradient descent using DiffKT."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Kotlin",
   "language": "kotlin",
   "name": "kotlin"
  },
  "language_info": {
   "codemirror_mode": "text/x-kotlin",
   "file_extension": ".kt",
   "mimetype": "text/x-kotlin",
   "name": "kotlin",
   "nbconvert_exporter": "",
   "pygments_lexer": "kotlin",
   "version": "1.7.0-dev-3303"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}