03_linear_gradient_descent.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "gpuClass": "standard"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "#### Recap"
      ],
      "metadata": {
        "id": "vu9d5t3-QuPc"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We have three data points and we want to fit a `linear model` to explain them\n",
        "\n",
        "$$(x_1, y_1)=(3,5),\\, (x_2, y_2)=(4,6),\\, (x_3, y_3)=(1,4)$$\n",
        "\n",
        "We write the model as\n",
        "\n",
        "$$y=A\\theta$$\n",
        "\n",
        "where\n",
        "\n",
        "$$y=\\begin{bmatrix} 5 \\\\ 6 \\\\ 4 \\end{bmatrix},\\, A=\\begin{bmatrix} 3 & 1 \\\\ 4 & 1 \\\\ 1 & 1\\end{bmatrix},\\, \\theta=\\begin{bmatrix} a \\\\ b \\end{bmatrix}$$\n",
        "\n",
        "and we can use the following to obtain $\\theta$\n",
        "\n",
        "$$\\theta=(A^TA)^{-1}A^Ty$$\n",
        "\n",
        "where $A^T$ is `transpose` of $A$\n",
        "\n",
        "However, the problem with such method is that computing inverse is not (always) easy or desirable, especially when the size of matrix is large\n",
        "\n",
        "As a result, in practice, `iterative numerical` approach is often used"
      ],
      "metadata": {
        "id": "E4EKYzZdSRf3"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### `Numerical` approach"
      ],
      "metadata": {
        "id": "MSzgkb2BRAtP"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Numerical method is commonly done `iteratively`\n",
        "\n",
        "* We `guess` $a$ and $b$ initially\n",
        "* We plug in guessed values and compute $A\\theta$ to get `predicted` output $y_{pred}$\n",
        "* We compared the output $y_{pred}$ with the `actual` data $y_{act}$, and compute the difference $\\|A\\theta-y_{act}\\|^2$\n",
        "* We call this difference the `cost function` or `loss function` $J$\n",
        "* Given data, $J$ is a `function of model parameters` $a$ and $b$\n",
        "* We adjust $a$ and $b$ in the next iteration in a `direction` such that $J$ will \"likely\" be reduced\n",
        "* This direction is the direction of `negative gradient` of $J$ w.r.t. $a$ and $b$, as the gradient of function always points to the direction where its value increases\n",
        "\n",
        "This algorithm is refered to as `gradient descent`"
      ],
      "metadata": {
        "id": "Wezg0HPWUAXB"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Derivation of `gradient` w.r.t. parameters"
      ],
      "metadata": {
        "id": "teq4PTiiRgpu"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "To compute the loss function, we can plug in each data sample and compute the prediction error corresponding to each, before adding them together\n",
        "\n",
        "$$J=\\frac{1}{m}\\sum_{i=1}^{m}\\,(y_{\\text{pred},i}-y_{\\text{act},i})^2$$\n",
        "\n",
        "where $m$ is the number of data samples (here, $m=3$) used to find the model parameters, and $y_{\\text{pred,i}}=a\\cdot x_i+b$\n",
        "\n",
        "There are also other options to define loss function, see `notes`"
      ],
      "metadata": {
        "id": "kQcBtwixf23h"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "The `gradient` w.r.t. $a$ and $b$ is then\n",
        "\n",
        "$$\\begin{align*}\n",
        "\\frac{\\partial J}{\\partial a}&=\\frac{1}{m}\\sum_{i=1}^{m}\\,(a\\cdot x_i+b-y_{\\text{act},i})\\cdot x_i\\\\\n",
        "\\frac{\\partial J}{\\partial b}&=\\frac{1}{m}\\sum_{i=1}^{m}\\,(a\\cdot x_i+b-y_{\\text{act},i})\\cdot 1\n",
        "\\end{align*}$$"
      ],
      "metadata": {
        "id": "K0UJJwjsgv8P"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Matrix notation"
      ],
      "metadata": {
        "id": "N4WHljXKR8bY"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "These two equations can be written in `matrix` form\n",
        "\n",
        "For example, `dJa=np.mean((y_pred)-(y_act))*x` where `y_pred` is a vector\n",
        "\n",
        "$$\\text{y_pred}=\\begin{bmatrix} y_{\\text{pred},1} \\\\ y_{\\text{pred},2} \\\\ y_{\\text{pred},3} \\end{bmatrix}$$\n",
        "\n",
        "`y_act` is also a vector\n",
        "\n",
        "$$\\text{y_act}=\\begin{bmatrix} y_{\\text{act},1} \\\\ y_{\\text{act},2} \\\\ y_{\\text{act},3} \\end{bmatrix}$$\n",
        "\n",
        "so is `x`\n",
        "\n",
        "$$\\text{x}=\\begin{bmatrix} x_1 \\\\ x_2 \\\\ x_3 \\end{bmatrix}$$\n",
        "\n",
        "When we do `dJa=np.mean((y_pred)-(y_act))*x`, numpy first carries out element-wise operation\n",
        "\n",
        "$$\\begin{bmatrix} (y_{\\text{pred},1}-y_{\\text{act},1})\\cdot x_1 \\\\ (y_{\\text{pred},2}-y_{\\text{act},2})\\cdot x_2 \\\\ (y_{\\text{pred},3}-y_{\\text{act},3})\\cdot x_3 \\end{bmatrix}$$\n",
        "\n",
        "Then, it takes the `average` of all elements, making it equivalent to $\\frac{\\partial J}{\\partial a}$ as described above\n",
        "\n",
        "Same goes for $\\frac{\\partial J}{\\partial b}$"
      ],
      "metadata": {
        "id": "lJz-erGGhY7X"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Adjustment of model parameters"
      ],
      "metadata": {
        "id": "-jsoNk1vcWZD"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "With the gradient, we can adjust $a$ and $b$ following\n",
        "\n",
        "$$a \\leftarrow a-\\eta\\frac{\\partial J}{\\partial a}, b \\leftarrow b-\\eta\\frac{\\partial J}{\\partial b}$$\n",
        "\n",
        "where $\\eta$ is the step size of adjustment that is often called `learning rate`"
      ],
      "metadata": {
        "id": "R4lhHatPavMB"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Put everything together"
      ],
      "metadata": {
        "id": "t9cqwBoilIod"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "plt.style.use('dark_background')\n",
        "np.set_printoptions(formatter={'float': '{: 0.4f}'.format})"
      ],
      "metadata": {
        "id": "nJgOaNMgTmV0"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "data_points = np.array([[3, 5], [4, 6], [1, 4]])\n",
        "x = data_points[:, 0] # a vector of three elements (3,)\n",
        "y = data_points[:, 1] # also a vector of three elements (3,)"
      ],
      "metadata": {
        "id": "G60A9dFAWwUf"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "learning_rate = 0.05\n",
        "iterations = 2000\n",
        "\n",
        "# Initial guess\n",
        "a = 0.0\n",
        "b = 0.0"
      ],
      "metadata": {
        "id": "oh4IkRA-lGQP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "for i in range(iterations):\n",
        "    # Make predictions using current a and b\n",
        "    y_pred = a * x + b\n",
        "\n",
        "    # Compute gradients using vectors directly\n",
        "    da = np.mean((y_pred - y) * x)\n",
        "    db = np.mean(y_pred - y)\n",
        "\n",
        "    # Update a and b in negative gradient direction\n",
        "    a -= learning_rate * da\n",
        "    b -= learning_rate * db\n",
        "\n",
        "# Final parameters\n",
        "print(f\"Final parameters: a = {a:.4f}, b = {b:.4f}\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Gj7-k7KplSqT",
        "outputId": "4dee925d-fdb8-4a55-e21c-e7d59c8596e3"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Final parameters: a = 0.6429, b = 3.2857\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "To see whether our `numerical` approach is doing well, we can use NumPy command `np.linalg.inv()` to compute the results analytically"
      ],
      "metadata": {
        "id": "enxVROePXMx-"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "A = np.hstack([x.reshape(-1, 1), np.ones((len(x), 1))]) # this is a matrix of size (3,2)\n",
        "theta = np.linalg.inv(A.T@A)@A.T@y\n",
        "theta"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "uAleBa5tW86g",
        "outputId": "66c9cb4c-f27a-4a23-8727-09525f4e9c1f"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([ 0.6429,  3.2857])"
            ]
          },
          "metadata": {},
          "execution_count": 5
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "We see the results are the same"
      ],
      "metadata": {
        "id": "vsLZDFD0ma_X"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Finally, we plot the results"
      ],
      "metadata": {
        "id": "NiRHaz5UXYT8"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def prediction(x, theta):\n",
        "    return theta[0]*x + theta[1]\n",
        "\n",
        "fig, ax=plt.subplots(figsize=(6,4))\n",
        "\n",
        "plt.scatter(A[:,0], y, color='cyan', label='Data', edgecolors='w')\n",
        "ax.plot([0,6],[prediction(0, theta), prediction(6, theta)],'r-', label='Model', zorder=-5)\n",
        "\n",
        "ax.set_xlabel('x')\n",
        "ax.set_ylabel('y')\n",
        "plt.tight_layout()\n",
        "plt.legend()\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 407
        },
        "id": "x9RauGzfXWir",
        "outputId": "be39fb21-7bf8-4367-ee90-ac522d12ec6e"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 600x400 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAGGCAYAAACNCg6xAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABAv0lEQVR4nO3deXRUVbr+8W9KCEoowCuakEDCJDMyBGyihqBAC+1ElAuKtiDtANqKY2vo9sa2r9AOLSIi2qCI4g+19SIqKggCRgg4MENA5mAMARSSICEJyf79cYpIgEAlpLJreD5rvauzT52qeqtWq8/aZ9c+YYBBRERERE7LZbsBERERkUCh4CQiIiLiJQUnERERES8pOImIiIh4ScFJRERExEsKTiIiIiJeUnASERER8ZKCk4iIiIiXatluwIbo6Gjy8/NttyEiIiJ+wu1289NPP532vJALTtHR0WRlZdluQ0RERPxMTEzMacNTyAWnozNNMTExmnUSERER3G43WVlZXuWCkAtOR+Xn5ys4iYiISKVocbiIiIiIlxScRERERLyk4CQiIiLipZBd43Q6devWpVGjRoSFhdluJegYY9i3bx+HDh2y3YqIiEilKDgdJywsjNtuu43evXvbbiXoLVq0iGnTpmGMsd2KiIiIVxScjnPbbbeRlJTEu+++y8aNGzly5IjtloJOrVq1aNu2LYMHDwbg9ddft9yRiIiIdxScjhEREUHv3r159913mTNnju12gtrWrVsBGDJkCO+8844u24mISEDQ4vBjnHfeeQBs3LjRcieh4ej33KhRI8udiIiIeEfB6RhHF4Lr8lzNOPo9awG+iIgECgUnERERES8pOImIiIhfcwG34B+hxR96kGpw9Gf9xhiKiorYvXs38+bN47bbbqvUpbBhw4axf/9+H3YqIiLivS7AMuAt4B67rQAKTj7jcrlISkrixhtvJCkpCZfL91/1Z599RlRUFM2aNWPAgAEsXLiQCRMm8Mknn3DWWWf5/P1FRESqS13gWeBboAdwAMiz2ZCHgpMPJCcns2XHDhYtWsTMmTNZtGgRW3bsIDk52afvW1hYSE5ODj/99BMrV65k3LhxXHfddfzhD39g+PDhADzwwAOsWbOGgwcPkpmZyaRJk4iIiAAgKSmJN954g4YNG5bNXqWmpgJwyy238O2335KXl0d2djZvv/02559/vk8/j4iIhKYBwHrgYZx9k94F2gHTbTbloeBUzZKTk3n//fdZGxNDT6Ae0BNYGxPD+++/7/PwdLyFCxeyatUqrr/+egBKS0u577776NChA8OGDeOKK67gmWeeAWDp0qWMHj2a3NxcoqKiiIqK4rnnngOgdu3aPP7443Tu3JmBAwfSrFkz3njjjRr9LCIiEtwigZnAp0AzYCdwFXAjsNteWycwoVRut9sYY4zb7T7hsbi4OPPmm2+auLi4Kr22y+Uy2zIzzeySEhN2dMrGU2HGmNklJWbrzp3G5XJV++eaNm2amTVr1kkfmzlzplm/fv1JH7vhhhvM3r17y8bDhg0z+/fvP+37xcfHG2OMiYiIqHLPZ/p9q1QqlSo4KgzMHWB+AWPAHAHzHJiIGnr/U2WD40szTtUoMTGR5k2bMtblwhz3mAHGuVy0iI0lMTGxRvsKCwsrux9cnz59mD9/Pj/++CN5eXm89dZbNGrUiHPOOeeUr9GtWzc++ugjdu7cSV5eHosXLwYgNjbW5/2LiEjwagssBv4NnAt8h7Om6WHgV4t9VUTBqRo1btwYgHUVPL7uuPNqSrt27di+fTtxcXF88sknrFmzhhtuuIH4+Hjuucf5jUJ4eHiFz69bty5z584lLy+Pm2++mR49epRdcjzV80RERCpSB/g7sBpIBA4CD+Asb1lpsa/TsRqctm/fznFXtDDG8NJLL1X4nEGDBpGRkUFBQQFr1qxhwIABNdjxqWVnZwPQsYLHOx53Xk24/PLLueiii/jggw+Ij4/H5XLx0EMPsXz5cjZv3kx0dHS584uKik74BV7btm1p1KgRjz32GF9//TWbNm3iggsuqLHPICIiwSUJJzD9DxAOfAy0B14ASuy15RWrwalHjx5li5CjoqLo27cvAP/5z39Oen5CQgIzZ87ktddeo2vXrnz44Yd8+OGHdOjQoSbbrlBaWhrbd+1iTGkpx++cFAaklJayLTOTtLQ0n7x/nTp1iIyMJDo6mq5du5KSksLs2bP5+OOPefPNN9myZQvh4eHce++9NG/enFtuuYWRI0eWe40dO3bgdru54oorOO+88zjnnHPIzMyksLCw7HnXXHMNjz/+uE8+g4iIBK//Al4DFgFtgJ+AQcC1wC57bVWa9UVhR2v8+PFm8+bNFT7+zjvvmI8//rjcsfT0dDN58uRqWQBWHYuVk5OTTUlJiZldUmJ6GmPqGWN6ehaGl5SUmOTkZJ98d9OmTTNHFRUVmZycHDNv3jwzfPhwExYWVnbe/fffb7Kyssyvv/5qPvvsM3PLLbcYY4xp0KBB2Tkvv/yy2bt3rzHGmNTUVAOYG2+80Wzbts0UFBSYJUuWmKuvvtoYY0znzp2r3LMWh6tUKlXo1M1g9uAs/jZgJoGp7wd9QeUWh2O72aNVu3Zts3fvXpOSklLhOTt37jSjR48ud+yJJ54wq1atqpYvp7r+Q56cnGy2ZWaaY23dudNnoSlQS8FJpVKpgr9agJnLb4FpLZgEP+jr2KpMcKqFnxg4cCANGzY85d5AUVFR5OTklDuWk5NDVFRUhc8JDw+nTp06ZWO3233GvZ7OrFmzmD17NomJiTRu3Jjs7GzS0tIoLS31+XuLiIj4g1rAQ0AqcA5wGPgHzm7gxRb7OlN+E5z+9Kc/8dlnn1X7wumUlBSeeOKJan1Nb5SWlpb9ZF9ERCSU/A5ne4GLPOMFwEhgi7WOqo9fbEcQGxtL3759mTp16inP2717N5GRkeWORUZGsnt3xfuJjhs3jvr165dVTExMtfQsIiIi5bmBicBSnNC0D7gV6EtwhCbwk+B02223sWfPHubMmXPK89LT0+nTp0+5Y/369SM9Pb3C5xQVFZGfn1+uREREpHoNBDKAP+OEi+k4m1u+ZbEnX7AenMLCwrjtttuYPn06JSXld2+YPn06Y8eOLRtPmDCB/v378+CDD9KmTRtSU1Pp3r37Kfd9EhEREd9pAszyVAywGegDDAd+tteWz1gPTn379iUuLo7XX3/9hMdiY2PL7bKdnp7O0KFDufPOO1m9ejWDBg1i4MCBrF+/viZbFhERCXku4F5gA85sUzHwvziX6L6015bPWV8c/sUXXxAWdvx2kY7LL7/8hGPvv/8+77//vq/bEhERkQpcBEwBLvaMlwJ3AqEwjWF9xklEREQCQ13gaeB7nNCUi/NrucsIjdAECk7ipaSkJIwxNGjQwOvnbN++ndGjR/uwKxERqSlX4tys/i84l6v+A7QDXsXZGTJUKDgFiWnTpmGMYfLkySc89tJLL2GMYdq0aRY6ExGRQHYB8DbwOdAcyASuBgYDNXfLev+h4BREMjMzufHGGzn77LPLjtWpU4ehQ4eyc+dOi52JiEigCQP+BGwEhgIlwPNAe+DUmwcFNwWnILJixQp27drF9ddfX3bs+uuvJzMzk5UrV5YdCw8PZ8KECeTk5FBQUEBaWhrdu3cv91oDBgxg06ZNHDp0iC+//JJmzZqd8H6XXnopX331FYcOHSIzM5MJEyZQt25dn30+ERGpGW2BRcBU4FxgBc6apoeAX+215RcUnILM66+/zm233VY2HjFixAmX6J555hluuOEGhg0bRrdu3diyZQtz587l3HPPBaBJkyb83//9Hx9//DFdunRh6tSp/POf/yz3Gi1atODzzz/ngw8+4KKLLmLIkCFcdtll2lNLRCSAhePcW2410AsnJD2IE5pWWOzL31i/K3FN1qnugBwXF2fefPNNExcXV+54XUtVmc81bdo0M2vWLNOoUSNTUFBgYmNjTWxsrDl06JA577zzzKxZs8y0adNM3bp1TWFhobnpppt+u9NzrVrmxx9/NA8//LABzFNPPWXWrVtX7vXHjRtnjDGmQYMGBjBTpkwxr7zySrlzLr30UnPkyBFTp04dA5jt27eb0aNHV9hzRd+3SqVSqWq+eoHJAGM89QmYOD/oqybqVNng+LK+j5O/q4u9ackI4FAln7Nv3z7mzJnD8OHDCQsLY86cOfz88297t7Zs2ZLw8HCWLFlSduzIkSN88803tGvXDoB27dqxfPnycq97/G1tOnfuzEUXXcTNN99cdiwsLIyzzjqL5s2bs3Hjxkp2LiIiNpwLPIuznglgN3Afzq/m5EQKTkHo9ddfL7tkds899/jkPerVq8err77Kiy++eMJjmZmZPnlPERGpXjcB44FIz/gV4DGc/Znk5BScTuMQzsyPrfeuis8//5zw8HCMMcydO7fcY1u3bqWwsJBLL720LODUqlWLHj168MILLwCQkZHBtddeW+55PXv2LDdesWIF7du3Z+vWrVXsUkREbGkOvAz094zXA3cBSyp8hhyl4OSFqgYYW0pLS8suu5WWlpZ77NChQ0yePJlnn32WX375hczMTP7yl79Qt25dXnvtNQBeeeUVHnroIZ555hmmTp1KfHw8w4cPL/c6Tz/9NMuWLWPixIlMnTqVX3/9lfbt29OvXz/uvffeGvmcIiJSObWAB4AncJaiHMa5v9wzOPeak9PTr+qCVH5+Pvn5+Sd97LHHHuODDz7grbfeYsWKFbRq1Yorr7ySAwcOALBr1y5uuOEGBg4cyOrVqxk5ciRjxowp9xpr164lKSmJ1q1bk5aWxsqVK3nyySf56aeffP3RRESkCnoA3+GEpLo4N+K9CHgKhabKsr6avSarKr+qU/mm9H2rVCqV78sNZgKYEpxfy+0DM8wP+vKn0q/qREREhOuAl4AmnvGbOJtY7rPWUeBTcBIREQkyMcBEINkz3gKMBBZY6yh4aI2TiIhIkHAB9wAbcEJTMTAW6IRCU3XRjJOIiEgQ6AT8Gzi6eUw6cCewzlpHwUkzTiIiIgHsHGAczr3keuJsXnk3cCkKTb6gGadjGGMAZ0NI8b2j3/PR711ERCqnH85u3y084/dxbpeSba2j4KcZp2Mcvadb27ZtLXcSGo5+z/v26fcdIiKVcT4wA5iHE5p2AdcC/41Ck69pauUYv/76K4sWLWLw4MEAbNy4kSNHjljuKvjUqlWLtm3bMnjwYBYtWsShQ4G2N7uIiD23Ac8B/wWU4Px67nHgoM2mQoiC03GmTZsGwJAhQyx3EvwWLVpU9n2LiMiptQZeBXp7xiuBO4DvbTUUosJwdsIMGW63m7y8POrXr1/hLUkA6tatS6NGjQgLC6vB7kKDMYZ9+/ZppklExAvhwGPAGKAO8CuQCryAM+MkZ87bbACacarQoUOHyMzMtN2GiIiEsMtwthho5xl/ivOLuZ3WOhItDhcREfEzDXECUxpOaMoBhgBXodBkm4KTiIiIHxkCZOCsXwInQLUF3rPWkRxLl+pERET8QDPgZWCAZ7wBuAv42lZDclKacRIREbHoLOBhYD1OaCoE/gfoikKTP7IenKKjo3nrrbfKfmW1Zs0a4uPjKzw/KSkJY8wJFRkZWYNdi4iInLnuwLfAs0BdYBFwEfAPoMheW3IKVi/VNWzYkCVLlrBw4UIGDBjA3r17ufDCC9m/f/9pn9u6dWvy8vLKxnv27PFlqyIiItWmHvC/wJ9xZpx+AR4C3rDYk3jHanB69NFH2bVrFyNGjCg7tmPHDq+eu2fPHnJzc33UmYiIiG9cA0wCmnrGM4AHgb3WOpLKsHqp7tprr+W7777jvffeIycnhxUrVnD77bd79dxVq1bx008/MW/ePC655JIKzwsPD8ftdpcrERGRmhaNcxPej3BC0zbg98AfUWgKNMZWFRQUmIKCAvPUU0+ZLl26mDvuuMMcOnTI3HrrrRU+p3Xr1ubOO+803bp1MwkJCea1114zRUVFpmvXric9PzU11ZyM2+229rlVKpVKFTrlAnM3mFwwBkwxmHFgzvGD3lROud3uymQDe40WFhaaJUuWlDs2YcIEs3Tp0kq9zqJFi8ybb7550sfCw8ON2+0uq+joaAUnlUqlUtVIdQSzFCcwGTDpYDr5QV+q8lWZ4GT1Ul12djYbNmwodywjI4PY2NhKvc4333xDq1atTvpYUVER+fn55UpERMSXzgbGAiuABCAPZyH4pcBai33JmbO6OHzJkiW0adOm3LHWrVuzc2flNpTv0qUL2dnZ1dmaiIhIlfQFXgFaesb/B9wL/GStI6lu1qbGunfvboqKikxKSopp2bKluemmm8zBgwfN0KFDy84ZO3asmT59etl49OjR5tprrzUtW7Y0HTp0MOPHjzdHjhwxV1xxRbVPx6lUKpVK5W01AvMmv12W2wXmWj/oS3X6Cpg1ToC56qqrzJo1a0xBQYHZsGGDuf3228s9Pm3aNLNw4cKy8SOPPGI2b95sDh06ZPbt22e+/PJL07t3b199OSqVSqVSnbaGg9mHE5hKwLwApp4f9KXyriqTDcI8f4QMt9tNXl4e9evX13onERE5IxcCrwKXe8argDtxdgOXwFGZbGD9lisiIiKBpjbwN2ANTmg6BPwF6IFCU7CzujhcREQk0FwK/Bto7xl/DowCdthqSGqUZpxERES80BDn13Jf44SmHOAmYAAKTaFEM04iIiKnMRiYAER5xlNxLs3tt9aR2KLgJCIiUoE4nBvyXuUZb8RZ/J1mrSOxTZfqREREjnMW8CCwHic0FQKpQGcUmkKdZpxERESOEY+z+LubZ7wYuAvYZK0j8SeacRIREQEigOeB5Tih6RfgTzjbDSg0yVGacRIRkZB3Nc5apqO3mP9/wAPAHmsdib9ScBIRkZDVGOfXcv/tGW8D7gbmWutI/J0u1YmISMgJA0YCGTih6QjwNNARhSY5Nc04iYhISOmAs/j7Es/4G5wtBlZb60gCiWacREQkJJwN/C+wEic05QP3AgkoNIn3NOMkIiJB7wqc26Vc6BnPwglNWdY6kkClGScREQlajYDpwAKc0JQFJAPXo9AkVaPgJCIiQelWnMXftwKlwESgHfChxZ4k8OlSnYiIBJVWOJfl+njGa3AWfy+31pEEE804iYhIUKgNjAHW4oSmAuBRnFuoKDRJddGMk4iIBLxLgFdx9mECmAeMwtnQUqQ6acZJREQCVgPgZWAJTmjaA9wMXIlCk/iGZpxERCQgDQJexLltCsBrwF9wbs4r4isKTiIiElBicW7Ie7VnvAm4C1hsrSMJJbpUJyIiAeEs4AFgPU5oKgL+DnRGoUlqjmacRETE73UFpuD8Qg7gK5xZpo3WOpJQpRknERHxWxHAc8C3OKFpP3AH0BuFJrFDM04iIuKX/oDzi7k4z3gmzqW6HGsdiSg4iYiIn4kCJgCDPeMdOHsyfW6rIZFjWL9UFx0dzVtvvcW+ffs4dOgQa9asIT4+/pTPSUpK4vvvv+fw4cNs3ryZYcOG1VC3IiKhzeVykZSUxI033khSUhIuV/X9ZyQMZ91SBk5oOgI8C3RAoUn8i7FVDRs2NNu3bzevv/666dGjh2nWrJnp16+fadGiRYXPadasmTl48KB57rnnTNu2bc0999xjiouLze9//3uv3tPtdhtjjHG73dY+t0qlUgViJScnm22ZmeZY2zIzTXJy8hm/dnswX4MxnvoGTGc/+Myq0KhKZgN7jY4bN8589dVXlXrOP//5T7N27dpyx2bOnGk+++wzX3w5KpVKpcIJTSUlJWZ2SYn5nTEmwhjzO2PM7JISU1JSUuXwVAfMP8AU4gSmPDD3gXH5wWdWhU4FTHBav369ef755817771ncnJyzIoVK8ztt99+yucsXrzYjB8/vtyx4cOHmwMHDvjiy1GpVKqQL5fLZbZlZprZJSUmzBjDMRXmCU9bd+40LperUq97OZgf+G2W6UMwTfzg86pCryqTDayucWrRogWjRo1i8+bNXHnllUyePJkXX3yRW2+9tcLnREVFkZNT/jcVOTk5NGjQgLPPPvuE88PDw3G73eVKRES8l5iYSPOmTRnrcmGOe8wA41wuWsTGkpiY6NXrnQdMA74ELgR+Aq4HBgI/VlvXIr5hNTi5XC5WrFjBX//6V1atWsWUKVOYMmUKI0eOrLb3SElJIS8vr6yysrKq7bVFREJB48bO3eDWVfD4uuPOO5VbcBZ/DwdKcW6d0g6YdYY9itQUq8EpOzubDRs2lDuWkZFBbGxshc/ZvXs3kZGR5Y5FRkaSm5vL4cOHTzh/3Lhx1K9fv6xiYmKqp3kRkRCRnZ0NQMcKHu943Hkn0xL4AngLOB9YC1wK/BnIq6Y+RWqC1eC0ZMkS2rRpU+5Y69at2blzZ4XPSU9Pp0+fPuWO9evXj/T09JOeX1RURH5+frkSERHvpaWlsX3XLsaUlhJ23GNhQEppKdsyM0lLSzvhubWBFJyg1Bco8Iy7Act827aIz1hbjNW9e3dTVFRkUlJSTMuWLc1NN91kDh48aIYOHVp2ztixY8306dPLxke3I3j66adNmzZtzKhRo7QdgUqlUvm4jv1VXU9jTD1jTM/T/KquJ5g1/Lb4ex6YFn7wWVSq4ytgflUHmKuuusqsWbPGFBQUmA0bNpzwq7pp06aZhQsXljuWlJRkVqxYYQ4fPmy2bNlihg0b5qsvR6VSqVSeOtk+Tlt37jwhNNUHMwlMCU5g2gPmZj/oX6WqqCqTDcI8f4QMt9tNXl4e9evX12U7EZFKcrlcJCYm0rhxY7Kzs0lLS6O0tLTs8euBiUC0ZzwNeBj4peZbFfFaZbKB7lUnIiJeKy0tZfHixSccbwq8BFzrGf+Ac/uURTXWmUjNsH6vOhERCVwuYDSwASc0FQFPAheh0CTBSTNOIiJSJV2AfwM9POOvgTtx9mkSCVaacRIRkUqpCzwLfIsTmg7gXJbrhUKTBD/NOImIiNcGAC8DzTzjd4H7gd2W+hGpaQpOIiJyWpHAC8CNnvEO4G7gM0v9iNiiS3UiIlKhMOAOnEtwNwIlwHNABxSaJDRpxklERE6qHfAqkOgZf4ez+HultY5E7NOMk4iIlFMH+DuwCic0HcRZx9QThSYRzTiJiEiZJJxZpqO3X/8YuAfYZa0jEf+iGScREeG/gNdwNq1sA/wEDMLZ1FKhSeQ3Ck4iIiHuZmAjMMIzfhlnfdMH1joS8V+6VCciEqJaAJOB33vG63AWf6db60jE/2nGSUQkxNQCHsUJSr8HDgNjgG4oNImcjmacRERCyO9w7i93kWe8ABgJbLHWkUhg0YyTiEgIcAMTgaU4oWkfcCvQF4UmkcrQjJOISJBLxglNMZ7xdOAh4GdrHYkELgUnEZEg1QR4CbjOM96Mc1nuS2sdiQQ+XaoTEQkyLuBeYANOaCoG/hfnEp1Ck8iZ0YyTiEgQ6Yyz+Ptiz3gJcBew3lpHIsFFM04iIkGgLvA0zo14LwZycS7LJaLQJFKdNOMkIhLgrsTZyLK5Z/weMBrYba0jkeClGScRkQB1AfD/gM9xQlMmcDUwBIUmEV9RcBIRCTBhwJ9w7i93E1ACPA+0B+ZY7EskFOhSnYhIAGkLvAr08oxXAHd4/ldEfE8zTiIiASAcSAVW44SmX4EHcRaCKzSJ1BzNOImI+LleOLNMbT3jOcA9wE5rHYmELs04iYj4qXOBqcBinNC0GxiMswBcoUnEDqvBKTU1FWNMucrIyKjw/GHDhp1wfkFBQQ12LCJSM24CMnAWgQO8ghOe/mOtIxEBP7hUt27dOvr27Vs2PnLkyCnPz83NpU2bNmVjY4zPehMRqWnNcfZkutIzXo+z8/cSax2JyLGsB6cjR46Qk5Pj9fnGmEqdLyISCGrhLPZOxdkF/DDO/eWewbnXnIj4B+trnC688EKysrLYunUrM2bMoGnTpqc8v169euzYsYPMzEw+/PBD2rdvf8rzw8PDcbvd5UpExJ/0wLlVytM4oelLnBvyPoVCk4i/sRqcli9fzvDhw+nfvz+jRo2iefPmpKWlUa9evZOev2nTJkaMGMF1113HLbfcgsvlYunSpcTExFT4HikpKeTl5ZVVVlaWrz6OiEiluIEJwDKcm/P+DAwH+gCb7bUlIqdh/KUaNGhgDhw4YEaMGOHV+bVq1TKbN282Tz75ZIXnhIeHG7fbXVbR0dHGGGPcbrf1z6tSqUK3rgOzC4zx1HQwjfygL5UqFMvtdnudDayvcTpWbm4uP/zwA61atfLq/CNHjrBy5cpTnl9UVERRUVF1tSgickZigIlAsme8BRgJLLDWkYhUhvU1TseKiIigZcuWZGdne3W+y+WiU6dOXp8vImKLC/gzsAEnNBUDY4FOKDSJBBKrwenZZ5+lV69exMXFkZCQwKxZsygpKWHmzJkATJ8+nbFjx5ad//jjj9OvXz+aN29O165dmTFjBnFxcUydOtXWRxAROa1OwFKcmab6QDrQDfgrzq/nRCRwWL1U16RJE2bOnMl5553H3r17+frrr+nZsyf79u0DIDY2ltLS0rLzzz33XKZMmUJUVBT79+/n+++/55JLLjnlppkiIracg7O9wEM4/7LNBVJwNrM0FvsSkaoLI8T++XW73eTl5VG/fn3y8/NttyMiQer3OBtZtvCM3wfuA7SwQMT/VCYb+NUaJxGRQHc+MAOYixOaMoFrgP9GoUkkGCg4iYhUkxHARuBmoAR4AegAfGKxJxGpXn61HYGISCBqDbwK9PaMVwJ3AN/bakhEfEYzTiIiVRQO/A+wBic0/Qo8jHMLFYUmkeCkGScRkSq4DPg30M4z/hS4G9hprSMRqQmacRIRqYSGOIEpDSc07QaGAFeh0CQSChScRES8NARn8fcdnvHRGaf3rHUkIjVNl+pERE6jGfAyMMAz3gDcBXxtqyERsUYzTiIiFTgLZ7H3epzQVAg8DnRFoUkkVGnGSUTkJLrjXIrr6hkvwpll+sFWQyLiFzTjJCJyjHo4G1cuwwlNvwC3AZej0CQimnESESlzDTAJaOoZzwAeBPZa60hE/E2lZ5zeeOMNEhMTfdGLiIgV0Tg34f0IJzRtxblJ7x9RaBKR8iodnBo0aMD8+fP54YcfSElJITo62hd9iYj4nAtn08oM4AbgCDAO6AR8YbEvEfFflQ5OycnJxMTEMHnyZIYMGcKOHTv49NNPueGGG6hVS1f+RCQwdAKW4Fyaq4+zpqkbMAYosNiXiPg/cybVtWtX8+KLL5pDhw6ZPXv2mOeff960atXqjF7Tl+V2u40xxrjdbuu9qFSqmq+zwYwFUwTGgMkFczcYlx/0plKp7FRlssEZ/aouKiqKfv360a9fP0pKSvj000/p1KkTGzZs4P777z+TlxYRqXZ9gXVAClAb+D+cnb9fBkot9iUigaVSqaxWrVrm+uuvNx9//LEpLCw03377rbnrrrvKpbSBAweaX375xXqCPFlpxkmlCr1qBOZNnBkmA2YXmGv9oC+VSuUfVZlsUOlFSdnZ2bhcLmbOnMnFF1/M6tWrTzhn4cKFHDhwoLIvLSJS7YYDzwHn4cwqTQT+Bhy02JOIBK5KB6cHHniA//znPxQWFlZ4Tm5uLi1atDijxkREzsSFwKs4G1cCrALuBL611ZCIBIVKr3GaMWPGKUOTiIhNtXFmlNbghKZDwCNADxSaROTMaf8AEQkal+LcX669Z/w5MArYYashEQk6ulediAS8hsArwNc4oSkHuAkYgEKTiFQvzTiJSEAbDEwAojzjKcCjwH5rHYlIMFNwEpGAFIez/9IfPOONOIu/06x1JCKhQJfqRCSgnAU8CKzHCU2FQCrQGYUmEfE9zTiJSMCIx1n83c0zXgzcBWyy1pGIhBqrM06pqakYY8pVRkbGKZ8zaNAgMjIyKCgoYM2aNQwYMKCGuhURWyKA54HlOKHpF+BPONsNVDU0uVwukpKSuPHGG0lKSsLl0gS8iHjH2hbnqampZu3atSYyMrKszjvvvArPT0hIMMXFxebhhx82bdu2NU8++aQpLCw0HTp08Mm26iqVyn5dDWYnv90uZQaYC87wNZOTk822zExzrG2ZmSY5Odn651WpVDVflcwG9hpNTU01K1eu9Pr8d955x3z88cfljqWnp5vJkyf76stRqVSWqjGY//BbYNoK5vfV8LrJycmmpKTEzC4pMb8zxkQYY35njJldUmJKSkoUnlSqEKzKZAPrc9MXXnghWVlZbN26lRkzZtC0adMKz01ISGD+/Pnljs2dO5eEhARftykiNSQMGAlkAIOAI8DTQEdg3hm+tsvl4l8TJvAJMNDlYjnwK84lwIEuF58Az73wgi7biUiFrP7bYfny5QwfPpz+/fszatQomjdvTlpaGvXq1Tvp+VFRUeTk5JQ7lpOTQ1RU1EnPBwgPD8ftdpcrEfFPHXA2sZwMNMAJNPHAY0BBNbx+YmIizZs2ZazLhTnuMQOMc7loERtLYmJiNbybiAQjq8Hp888/5/3332ft2rXMmzePP/zhDzRs2JDBgwdX23ukpKSQl5dXVllZWdX22iJSPc4GngJWApcA+cC9nr/XVOP7NG7cGIB1FTy+7rjzRESO51fz0bm5ufzwww+0atXqpI/v3r2byMjIcsciIyPZvXt3ha85btw46tevX1YxMTHV2rOInJkrgLXAGJwb9M4C2gEvAaXV/F7Z2dmAc9nvZDoed56IyPH8KjhFRETQsmXLCv+llZ6eTp8+fcod69evH+np6RW+ZlFREfn5+eVKROxrBEwHFgCtgCwgGbje87cvpKWlsX3XLsaUlhJ23GNhQEppKdsyM0lL01aaIlIxa6vYn332WdOrVy8TFxdnEhISzLx588yePXtMo0aNDGCmT59uxo4dW3Z+QkKCKSoqMg8++KBp06aNSU1N1XYEKlUA1q1g9uL8Wq4EzItg3DX03sf+qq6nMaaeMaanflWnUoV0Bcx2BDNnzjRZWVnm8OHDZteuXWbmzJmmRYsWZY8vXLjQTJs2rdxzBg0aZDZu3GgOHz5s1q5dawYMGODLL0elUlVjtQIzn9+2GFgN5mILfZxsH6etO3cqNKlUIVqVyQZhnj9ChtvtJi8vj/r16+uynUgNqQ38BfgbzkLwQ8DfcXYDP2KpJ5fLRWJiIo0bNyY7O5u0tDRKS6t7VZWIBILKZAPdq05EfOoSnPvLdfCM5wKjgO3WOnKUlpayePFiy12ISKDxq8XhIhI8GuDsx7QEJzTtAW4G+mM/NImIVJVmnESk2g0CXgSO7ob0GvAIsN9aRyIi1UPBSUSqTSwwCbjaM94E3AXogpiIBAtdqhORM3YW8ACwHic0FeEs/u6MQpOIBBfNOInIGekKTMG5pxzAVzizTButdSQi4juacRKRKokA/gV8ixOa9gN3AL1RaBKR4KUZJxGptD8ALwNxnvFMnEt1OdY6EhGpGQpOIuK1KGACMNgz3g7cDXxurSMRkZqlS3UiclphOOuWMnBC0xHgWaAjCk0iElo04yQip9QeZ+fvSz3jb3HWMq221pGIiD2acRKRk6oD/ANYiROa8oHRQE8UmkQkdGnGSUROcDnwKnChZzwb+DPwo7WORET8g2acRKTMecA04Euc0JQFXA8MRKFJRAQUnETE4xac/ZeGA6U4t05pD8yy2JOIiL/RpTqRENcSeAXo6xmvBe4EllnrSETEf2nGSSRE1QZScIJSX6AAeAzohkKTiEhFNOMkEoIScBZ/d/KMvwBGAtusdSQiEhg04yQSQurjrF36Gic07cVZ2/R7FJpERLyhGSeREHE9MBGI9oynAQ8Dv1jrSEQk8Cg4iQS5psBLwLWe8Q84t09ZZKshEZEApkt1IkHKhbPT9wac0FQEPAlchEKTiEhVacZJJAh1AaYA3T3jr3G2GMiw1ZCISJDQjJNIEKkLPAt8hxOaDuAEpl4oNImIVAfNOIkEiQHAy0Azz/hd4H5gt6V+RESCkYKTSICLBCYAQzzjHcDdwGe2GhIRCWK6VCcSoMKAO3AuwQ0BSoDngA4oNImI+IpmnEQCUDucnb8TPePvcNYyrbTWkYhIaPCbGadHH30UYwzjx4+v8Jxhw4ZhjClXBQUFNdiliF11gL8Dq3BC00GcdUw9UWgSEakJfjHj1L17d+666y5Wr1592nNzc3Np06ZN2dgY48vWRPxGb+AV4Oj/+z8G7gF22WpIRCQEWZ9xioiI4O233+aOO+5g//79pz3fGENOTk5Z7dmzpwa6FLHnv4DXgIU4oeknYBDOppYKTSIiNct6cJo0aRJz5sxhwYIFXp1fr149duzYQWZmJh9++CHt27c/5fnh4eG43e5yJRIobgY2AiOAUpztBtoBH9hsSkQkhFkNTkOGDKFbt26kpKR4df6mTZsYMWIE1113Hbfccgsul4ulS5cSExNT4XNSUlLIy8srq6ysrOpqX8RnWgBzgRnA+cA64DKcS3N5FvsSEREwNqpJkyZm9+7dplOnTmXHFi5caMaPH+/1a9SqVcts3rzZPPnkkxWeEx4ebtxud1lFR0cbY4xxu91WPrdKdaqqBeZRMIfAGDAFYFLA1PaD3lQqlSpYy+12e50NrC0Oj4+PJzIykhUrVpQdq1WrFr169eLPf/4zderUobS09JSvceTIEVauXEmrVq0qPKeoqIiioqJq61vEV34H/BvnJrwAC4CRwBZrHYmIyPGsBacFCxbQsWPHcsemTZvGxo0befrpp08bmgBcLhedOnXi008/9VWbIj7nBsbi7PbtAvYBDwJv2WxKREROylpwOnjwIOvXry937Ndff+Xnn38uOz59+nSysrIYM2YMAI8//jjLli1jy5YtNGzYkEceeYS4uDimTp1a4/2LVIdkYCJwdJXeG8DDwM+2GhIRkVPyi32cKhIbG1tu5uncc89lypQpREVFsX//fr7//nsuueQSMjJ033cJLE2Al4DrPOPNOJflvrTWkYiIeCMMZ7FTyHC73eTl5VG/fn3y8/NttyMhxoXzy7incC7RFQNPe8aHLfYlIhLKKpMN/HrGSSSYdMZZ/H2xZ7wE5/5yG6x1JCIilWV9A0yRYFcXeAbnRrwXAwdwLsslotAkIhJoNOMk4kNXApOB5p7xe8BoYLe1jkRE5ExoxknEBy4A/h/wOU5oygSuBoag0CQiEsgUnESqURjwJ5z7y90ElADPA+2BORb7EhGR6qFLdSLVpC3wKtDLM14B3OH5XxERCQ6acRI5Q3WAJ4DVOKHpV5ydvy9GoUlEJNhoxknkDPTCmWVq6xl/grNPU6a1jkRExJc04yRSBecCU4HFOKEpGxgMXINCk4hIMFNwEqmkm3AWf//JM34FaAf8x1pHIiJSU3SpTsRLzXH2ZLrSM14P3IWzA7iIiIQGzTiJnEYt4C/AOpzQdBj4G9AVhSYRkVCjGSeRU+gBTMG5zxzAlzi3S9lsrSMREbFJM04iJ+EGXgSW4YSmn4HhQB8UmkREQplmnESOcx3wEtDEM34TeAjYZ60jERHxFwpOIh4xwEQg2TPegnNZboG1jkRExN/oUp2EPBfwZ2ADTmgqBp4COqHQJCIi5WnGSULaRcC/gd95xktxthhYZ60jERHxZ5pxkpB0DvBP4Huc0JQLjAIuQ6FJREQqphknCTm/x9nIsoVn/D5wH85tU0RERE5FM04SMs4HZgBzcUJTJs695f4bhSYREfGOgpOEhBE495e7GSgBXgA6AJ9Y7ElERAKPLtVJUGsNvAr09oxXAnfgrG0SERGpLM04SVAKB/4HWIMTmn7F2cSyBwpNIiJSdZpxkqCTiDPL1M4z/hS4G9hprSMREQkWmnGSoNEQZ0+mr3BC025gCHAVCk0iIlI9FJwkKAzBWfx9h2f8b5zw9J61jkREJBj5TXB69NFHMcYwfvz4U543aNAgMjIyKCgoYM2aNQwYMKCGOhRfcblcJCUlceONN5KUlITL5f3/LZvhXIp7B4jEuW1KIs7u3weqv1URERGM7erevbvZtm2bWbVqlRk/fnyF5yUkJJji4mLz8MMPm7Zt25onn3zSFBYWmg4dOnj9Xm632xhjjNvttv65VZjk5GSzLTPTHGtbZqZJTk4+5fPOAvMwmF/BGDCHwfwNTLgffCaVSqVSBVZVMhvYbTYiIsJs2rTJ9OnTxyxcuPCUwemdd94xH3/8cblj6enpZvLkyb76clQ+rOTkZFNSUmJml5SY3xljIowxvzPGzC4pMSUlJRWGp+5gVuIEJgPmSzCt/eDzqFQqlSowqzLZwPqlukmTJjFnzhwWLDj9fegTEhKYP39+uWNz584lISHBV+2Jj7hcLv41YQKfAANdLpbjbBmw3DP+BHjuhRfKXbarh7Nx5XKgC/AzcBtwBfBDjXYvIiKhyup2BEOGDKFbt2706NHDq/OjoqLIyckpdywnJ4eoqKgKnxMeHk6dOnXKxm63u2rNSrVKTEykedOm3IQT4Y9lgHEuF+mxsSQmJrJ48WKuBV4CmnrOmQE8COytuZZFRETsLQ5v0qQJEyZM4Oabb6awsNBn75OSkkJeXl5ZZWVl+ey9xHuNGzcGYF0Fjx893r5xYz4AZuOEpq1AP+CPKDSJiEjNsxac4uPjiYyMZMWKFRQXF1NcXEzv3r257777KC4uPukvq3bv3k1kZGS5Y5GRkezevbvC9xk3bhz169cvq5iYmGr/LFJ52dnObXU7VvD40ePPZmdzPVAMjAM6AfMreI6IiIivWQtOCxYsoGPHjnTp0qWsvv32W95++226dOlCaWnpCc9JT0+nT58+5Y7169eP9PT0Ct+nqKiI/Pz8ciX2paWlsX3XLsaUlhJ23GNhwJjSUoozM4lIS2MZEA+MAQpqvFMREZHyrK9mP1rH/6pu+vTpZuzYsWXjhIQEU1RUZB588EHTpk0bk5qaqu0IAriO/VVdT2NMPWNMT2PMRyUlprSkxPyanGzuBuPyg15VKpVKFbwVUL+qO5XY2NiytTDgzDgNHTqUO++8k9WrVzNo0CAGDhzI+vXrLXYpVTVr1iwGDRpEp6ws0oF8IB3o/+OPLB80iAtnzeJl4MS5RxERETvCcBJUyHC73eTl5VG/fn1dtvMDjYDxLhe3JCZC48bszc7mzrQ0PjzJpVoRERFfqEw2sLodgYS24cBzwHmlpZQuXsxE4G/AQatdiYiIVEzBSWrchcCrwOWe8SrgTuBbWw2JiIh4ya/XOElwCQceB9bihKZDwCNADxSaREQkMGjGSWrEpcC/gfae8WfA3cAOWw2JiIhUgWacxKcaAq8AX+OEphzgJuAPKDSJiEjg0YyT+MxgYAJw9E6CU4BHgf3WOhIRETkzCk5S7eKAl3FmlQA24iz+TrPWkYiISPXQpTqpNmcBDwLrcUJTIZAKdEahSUREgoNmnKRaxOMs/u7mGS8G7gI2WetIRESk+mnGSc5IPWA8sBwnNP0C/AlnuwGFJhERCTaacZIquxqYBMR6xm8DDwB7rXUkIiLiWwpOUmmNgReBQZ7xNmAUMM9aRyIiIjVDl+rEa2HASCADJzQdAZ4GOqLQJCIioUEzTuKVDjiLvy/xjJfjbDGwxlpHIiIiNU8zTnJKZwNPAStxQlM+cK/nb4UmEREJNZpxkgpdAbwKtPKMZ+GEpixrHYmIiNilGSc5QSNgOrAAJzT9CCQD16PQJCIioU3BScq5FWfx961AKTAR5+a8H1rsSURExF/oUp0AzszSK0Afz3g1zuLvb6x1JCIi4n804xTiagN/BdbihKZDwF+A7ig0iYiIHE8zTiHsEpwtBjp4xnNxNrLcbq0jERER/6YZpxDUAJgMLMEJTXuAm4H+KDSJiIicimacQswgnNulNPaMXwMeAfZb60hERCRwKDiFiFicG/Je7RlvBO4CvrLWkYiISODRpbogdxbwALABJzQVAU8AnVFoEhERqSzNOAWxbjiLv+M9469wZpk2WutIREQksGnGKQhFAP/C2U4gHmf90u1AbxSaREREzoRmnILMVThrmeI845k4l+pyrHUkIiISPKzOOI0cOZLVq1eTm5tLbm4uS5cupX///hWeP2zYMIwx5aqgoKAGO/ZfUcC7wCc4oWk7MAAYikKTiIhIdbE64/Tjjz/y2GOPsXnzZsLCwhg2bBizZ8+ma9eubNiw4aTPyc3NpU2bNmVjY0xNteuXwnBujfJPoCFwBBiPswD8kLWuREREgpfxp/r555/NiBEjTvrYsGHDzP79+8/o9d1utzHGGLfbbf2znmm1B/M1GOOpb8B09oO+VCqVSqUKpKpMNvCbxeEul4shQ4YQERFBenp6hefVq1ePHTt2kJmZyYcffkj79u1P+brh4eG43e5yFejqAP8AVgKXAvnAfUBPnJvzioiIiO9YTXkdO3Y0+fn5pri42Ozfv98MGDCgwnN79uxp/vjHP5rOnTubXr16mY8++sgcOHDAxMTEVPic1NRUczKBOuN0OZgf+G2W6UMwTfygL5VKpVKpArUqeTXKbrO1a9c2LVu2NN26dTNjx441e/bsMe3atfPqubVq1TKbN282Tz75ZIXnhIeHG7fbXVbR0dEBGZzOAzON3wLTj2CS/aAvlUqlUqkCvSoTnKxvR1BcXMzWrVsBWLFiBT169GD06NGMHDnytM89cuQIK1eupFWrVhWeU1RURFFRUbX1a8MfgeeBRkAp8DLwVyDPZlMiIiIhyG/WOB3lcrmoU6eO1+d26tSJ7OxsH3dlR0vgC+BNnNC0FrgEuBeFJhERERuszjiNHTuWzz77jMzMTNxuN0OHDqV3795ceeWVAEyfPp2srCzGjBkDwOOPP86yZcvYsmULDRs25JFHHiEuLo6pU6fa/BjVrjbwMPA4cA5QAPwdZzfwIxb7EhERCXVWg9MFF1zAm2++SePGjcnNzWXNmjVceeWVzJ8/H4DY2FhKS0vLzj/33HOZMmUKUVFR7N+/n++//55LLrmEjIwMWx+h2iUArwKdPOMvgJHANmsdiYiIyFFhOIudQobb7SYvL4/69euTn59vu50y9YFxOCHJBezFuVXK2zabEhERCQGVyQbWF4cL3AC8CER7xq8DjwC/WOtIRERETkbByaKmODfkvcYz/gG4C1hkqyERERE5Jb/7VV0ocAGjgQ04oakIeBK4CIUmERERf6YZpxrWBZgCdPeMv8a5SW/wLG8XEREJXppxqiF1gWeB73BC0wGcwNQLhSYREZFAoRmnGjAAZ7fvZp7xu8D9wG5L/YiIiEjVKDj5UCQwARjiGe8A7gY+s9WQiIiInBFdqvOBMOAOnEtwQ4AS4DmgAwpNIiIigUwzTtWsHfBv4DLP+DuctUwrrXUkIiIi1UUzTtXoWmAVTmg6iLOOqScKTSIiIsFCM07VKA3n13LLgXuAXVa7ERERkeqm4FSN9gNdgZ9sNyIiIiI+oUt11UyhSUREJHgpOImIiIh4ScFJRERExEsKTiIiIiJeUnASERER8ZKCk4iIiIiXFJxEREREvKTgJCIiIuIlBScRERERLyk4iYiIiHhJwUlERETESwpOIiIiIl4K2Zv8ut1u2y2IiIiIH6hMJgi54HT0y8nKyrLciYiIiPgTt9tNfn7+Kc8JA0zNtOM/oqOjT/vFVJXb7SYrK4uYmBifvUcw0vdWdfruqk7fXdXoe6s6fXdVUxPfm9vt5qeffjrteSE34wR49cWcqfz8fP1DUQX63qpO313V6burGn1vVafvrmp8+b15+7paHC4iIiLiJQUnERERES8pOFWzwsJCnnjiCQoLC223ElD0vVWdvruq03dXNfreqk7fXdX40/cWkovDRURERKpCM04iIiIiXlJwEhEREfGSgpOIiIiIlxScqtHdd9/N9u3bKSgoYNmyZfTo0cN2SwEhMTGRjz76iKysLIwxXHfddbZbCgiPPfYY33zzDXl5eeTk5DBr1ixat25tuy2/N3LkSFavXk1ubi65ubksXbqU/v37224rID366KMYYxg/frztVvxeamoqxphylZGRYbutgBAdHc1bb73Fvn37OHToEGvWrCE+Pt5aPwpO1WTw4ME8//zz/P3vf6dbt26sXr2auXPncv7559tuze9FRESwevVq7rnnHtutBJSkpCQmTZpEz5496devH7Vr12bevHnUrVvXdmt+7ccff+Sxxx4jPj6e7t278+WXXzJ79mzat29vu7WA0r17d+666y5Wr15tu5WAsW7dOqKiosrqsssus92S32vYsCFLliyhuLiYAQMG0L59ex566CH2799vtS+jOvNatmyZmThxYtk4LCzM/Pjjj+bRRx+13lsglTHGXHfdddb7CMRq1KiRMcaYxMRE670EWv38889mxIgR1vsIlIqIiDCbNm0yffr0MQsXLjTjx4+33pO/V2pqqlm5cqX1PgKtxo0bZ7766ivrfRxbmnGqBrVr1yY+Pp758+eXHTPGMH/+fBISEix2JqGkQYMGAPzyyy+WOwkcLpeLIUOGEBERQXp6uu12AsakSZOYM2cOCxYssN1KQLnwwgvJyspi69atzJgxg6ZNm9puye9de+21fPfdd7z33nvk5OSwYsUKbr/9dqs9KThVg0aNGlGrVi1ycnLKHc/JySEqKspSVxJKwsLCeOGFF/j6669Zv3697Xb8XseOHcnPz6ewsJBXXnmF5ORkrTfx0pAhQ+jWrRspKSm2Wwkoy5cvZ/jw4fTv359Ro0bRvHlz0tLSqFevnu3W/FqLFi0YNWoUmzdv5sorr2Ty5Mm8+OKL3HrrrVb7sj7tFejVuHFjY4wxPXv2LHf86aefNsuWLbPeXyCVLtVVrV5++WWzfft2ExMTY72XQKjatWubli1bmm7dupmxY8eaPXv2mHbt2lnvy9+rSZMmZvfu3aZTp05lx3SprmrVoEEDc+DAAV0iPk0VFhaaJUuWlDs2YcIEs3TpUms9acapGuzbt48jR44QGRlZ7nhkZCS7d++21JWEiokTJ3L11Vdz+eWXk5WVZbudgFBcXMzWrVtZsWIFY8aMYfXq1YwePdp2W34vPj6eyMhIVqxYQXFxMcXFxfTu3Zv77ruP4uJiXC79J8Vbubm5/PDDD7Rq1cp2K34tOzubDRs2lDuWkZFBbGyspY50qa5aFBcX8/3339OnT5+yY2FhYfTp00frJsSnJk6cSHJyMldccQU7duyw3U7Acrlc1KlTx3Ybfm/BggV07NiRLl26lNW3337L22+/TZcuXSgtLbXdYsCIiIigZcuWZGdn227Fry1ZsoQ2bdqUO9a6dWt27txpqSOH9am4YKjBgwebgoICc+utt5q2bduaV155xfzyyy/mggsusN6bv1dERITp3Lmz6dy5szHGmPvvv9907tzZNG3a1Hpv/lyTJk0y+/fvN7169TKRkZFldfbZZ1vvzZ9r7NixJjEx0cTFxZmOHTuasWPHmpKSEtO3b1/rvQVi6VKdd/Xss8+aXr16mbi4OJOQkGDmzZtn9uzZYxo1amS9N3+u7t27m6KiIpOSkmJatmxpbrrpJnPw4EEzdOhQm33Z/2KCpe655x6zY8cOc/jwYbNs2TJz8cUXW+8pECopKcmczLRp06z35s9VkWHDhlnvzZ9r6tSpZvv27ebw4cMmJyfHfPHFFwpNZ1AKTt7VzJkzTVZWljl8+LDZtWuXmTlzpmnRooX1vgKhrrrqKrNmzRpTUFBgNmzYYG6//Xar/YR5/hARERGR09AaJxEREREvKTiJiIiIeEnBSURERMRLCk4iIiIiXlJwEhEREfGSgpOIiIiIlxScRERERLyk4CQiIiLiJQUnERERES8pOImIiIh4ScFJRERExEsKTiISUho1akR2djYpKSllxxISEigsLOSKK66w2JmIBArrdz5WqVSqmqwBAwaYwsJCEx8fb+rVq2e2bNli/vWvf1nvS6VS+X+Fef4QEQkpL730En379uW7776jU6dO9OjRg6KiItttiYifU3ASkZB09tlns27dOpo2bUp8fDzr1q2z3ZKIBACtcRKRkNSyZUuio6NxuVw0a9bMdjsiEiA04yQiIad27dp88803rFq1ik2bNnH//ffTqVMn9u7da7s1EfFzCk4iEnKeeeYZBg0aROfOnTl48CCLFy8mNzeXa665xnZrIhIArK9QV6lUqpqqpKQkU1RUZC699NKyY3FxcebAgQNm5MiR1vtTqVT+XZpxEhEREfGSFoeLiIiIeEnBSURERMRLCk4iIiIiXlJwEhEREfGSgpOIiIiIlxScRERERLyk4CQiIiLiJQUnERERES8pOImIiIh4ScFJRERExEsKTiIiIiJeUnASERER8dL/B9EHSH7nzsM9AAAAAElFTkSuQmCC\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### A slightly less empirical view of gradient descent..."
      ],
      "metadata": {
        "id": "FhI4UK9nceVe"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "At `kth` iteration, gradient descent as described above can be considered as the operations of\n",
        "\n",
        "* find the 1st order `approximation of loss function` at current $a^k, b^k$\n",
        "* find new $a^{k+1}, b^{k+1}$ `in the vicinity of` $a^k, b^k$ that `minimizes` this 1st order approximation"
      ],
      "metadata": {
        "id": "JP1wv2EHfVI6"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "To see this, suppose at `kth` iteration we have our parameters $\\theta^k=\\begin{bmatrix}a^k \\\\ b^k\\end{bmatrix}$\n",
        "\n",
        "We take the 1st order `Taylor approximation` of our cost function $J(\\theta)$ (recall $J$ is a function of model parameters) in the vicinity of $\\theta^k$\n",
        "\n",
        "$$\\begin{align*}J(\\theta)\\approx J(\\theta^k)+\\nabla J(\\theta^k)^T(\\theta-\\theta^k)\n",
        "\\end{align*}$$\n",
        "\n",
        "where $\\nabla J(\\theta^k)^T$ is the trasposed vector containing the gradients\n",
        "\n",
        "$$\\nabla J(\\theta^k)^T =\\begin{bmatrix} \\left.\\frac{\\partial J}{\\partial a}\\right|_{a=a^k} & \\left.\\frac{\\partial J}{\\partial b}\\right|_{b=b^k}\\end{bmatrix}$$"
      ],
      "metadata": {
        "id": "hSK6fNZmcmXE"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Since it is just an approximation, the idea is to minimize $J$ while not moving too far away from $\\theta^k$ (otherwise the approximation can break down)\n",
        "\n",
        "Therefore, we want to find $\\theta^{k+1}=\\begin{bmatrix}a^{k+1} \\\\ b^{k+1}\\end{bmatrix}$ for the `(k+1)th` iteration as a solution to minimize the following expression\n",
        "\n",
        "$$J(\\theta^{k})+\\nabla J(\\theta^{k})^T(\\theta-\\theta^{k})+\\frac{1}{2\\eta}\\|\\theta-\\theta^{k}\\|^2$$\n",
        "\n",
        "By having the following term in the expression\n",
        "\n",
        "$$\\frac{1}{2\\eta}\\|\\theta-\\theta^{k}\\|^2$$\n",
        "\n",
        "we basically ask to `not move too far away` from $\\theta^{k}$\n",
        "\n",
        "We take derivative of this expression w.r.t. $\\theta$ and set it to zero\n",
        "\n",
        "$$0+\\nabla J(\\theta^{k})+\\frac{1}{\\eta}(\\theta-\\theta^{k})=0$$\n",
        "\n",
        "and we get our `gradient descent equation`\n",
        "\n",
        "$$\\theta^{k+1}=\\theta^{k}-\\eta \\nabla J(\\theta^{k})$$\n",
        "\n",
        "or\n",
        "\n",
        "$$\\begin{bmatrix}a^{k+1} \\\\ b^{k+1}\\end{bmatrix}=\\begin{bmatrix}a^k \\\\ b^k\\end{bmatrix}-\\eta \\begin{bmatrix} \\left.\\frac{\\partial J}{\\partial a}\\right|_{a=a^k} \\\\ \\left.\\frac{\\partial J}{\\partial b}\\right|_{b=b^k}\\end{bmatrix}$$\n",
        "\n",
        "which is the same as described earlier\n",
        "\n",
        "$$a \\leftarrow a-\\eta\\frac{\\partial J}{\\partial a}, b \\leftarrow b-\\eta\\frac{\\partial J}{\\partial b}$$"
      ],
      "metadata": {
        "id": "c52E-s-vcoax"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "i4hwsSsonF89"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}