gradient_descent_on_least_squares.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4ff5888c",
   "metadata": {},
   "source": [
    "Sascha Spors,\n",
    "Professorship Signal Theory and Digital Signal Processing,\n",
    "Institute of Communications Engineering (INT),\n",
    "Faculty of Computer Science and Electrical Engineering (IEF),\n",
    "University of Rostock,\n",
    "Germany\n",
    "\n",
    "# Data Driven Audio Signal Processing - A Tutorial with Computational Examples\n",
    "\n",
    "Winter Semester 2023/24 (Master Course #24512)\n",
    "\n",
    "- lecture: https://github.com/spatialaudio/data-driven-audio-signal-processing-lecture\n",
    "- tutorial: https://github.com/spatialaudio/data-driven-audio-signal-processing-exercise\n",
    "\n",
    "Feel free to contact lecturer frank.schultz@uni-rostock.de"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5752a87",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ce8f602",
   "metadata": {},
   "source": [
    "# Least Squares Error with Gradient Descent\n",
    "\n",
    "##  Straightforward Gradient Descent (GD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36335ed5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for reproducible outputs\n",
    "rng = np.random.default_rng(12345)  # used for data creation and shuffling\n",
    "\n",
    "shuffled = True\n",
    "\n",
    "M = 2**10\n",
    "x = np.linspace(0, 2 * np.pi, M)  # lin increase\n",
    "\n",
    "if shuffled:\n",
    "    rng.shuffle(x)\n",
    "\n",
    "# design/feature matrix of the true model\n",
    "X = np.column_stack(\n",
    "    (np.cos(1 * x), np.sin(2 * x), np.cos(5 * x), np.cos(6 * x))\n",
    ")\n",
    "\n",
    "# add bias/intercept column to the design/feature matrix\n",
    "X = np.hstack((np.ones((M, 1)), X))\n",
    "N = X.shape[1]\n",
    "\n",
    "# some nice numbers for the true model parameters beta\n",
    "beta_true = np.array([3, 2, 1, 1 / 2, 1 / 4])\n",
    "# outcome of true model, i.e. linear combination\n",
    "y_true = (X @ beta_true)[:, None]\n",
    "# add measurement noise\n",
    "noise_scale = 0.5\n",
    "noise = rng.normal(loc=0, scale=noise_scale, size=(M, 1))\n",
    "y = y_true + noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a1f9c60",
   "metadata": {},
   "outputs": [],
   "source": [
    "beta = rng.normal(loc=0, scale=1, size=(N, 1))  # random init of beta coeff\n",
    "step_size = 1e-4\n",
    "steps = 2**6\n",
    "\n",
    "for step in range(steps):\n",
    "    # from the normal equations\n",
    "    beta_gradient = -2 * X.T @ y + 2 * X.T @ X @ beta\n",
    "    # we can rewrite to see:\n",
    "    # transposed feature matrix X times error term 'y - X @ beta'\n",
    "    # and weighted by -2\n",
    "    beta_gradient = -2 * X.T @ (y - X @ beta)\n",
    "    # the usual update rule\n",
    "    beta = beta - step_size * beta_gradient\n",
    "# check ratio between coeffs from GD and true coeff\n",
    "beta / beta_true[:, None]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ed55357",
   "metadata": {},
   "source": [
    "## Stochastic Gradient Descent (SGD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "859775a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# straighforward gradient descent parameters\n",
    "beta = rng.normal(loc=0, scale=1, size=(N, 1))  # random init of beta coeff\n",
    "step_size = 1e-4\n",
    "steps = 2**11  # more steps than above\n",
    "\n",
    "# additional stochastic gradient descent parameters / needed arrays\n",
    "idx = np.arange(M)\n",
    "Nrandom = 2**4  # but less data for each gradient calculation\n",
    "\n",
    "for step in range(steps):\n",
    "    random_idx = np.random.choice(idx, size=Nrandom, replace=False)\n",
    "    Xr, yr = X[random_idx, :], y[random_idx, :]\n",
    "    # from the normal equations\n",
    "    beta_gradient = -2 * Xr.T @ yr + 2 * Xr.T @ Xr @ beta\n",
    "    beta_gradient = -2 * Xr.T @ (yr - Xr @ beta)\n",
    "    # the usual update rule\n",
    "    beta = beta - step_size * beta_gradient\n",
    "# check ratio between coeffs from SGD and true coeff\n",
    "beta / beta_true[:, None]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d2c2d4e",
   "metadata": {},
   "source": [
    "## Stochastic Gradient Descent (SGD) with Momentum\n",
    "\n",
    "Try out at home: implement Stochastic Gradient Descent with momentum, cf. [gradient_descent_momentum.ipynb](gradient_descent_momentum.ipynb)\n",
    "\n",
    "## Further Reading\n",
    "\n",
    "We could also invent some adaptive rule for the `step_size`, i.e. large steps if far away from a suitable minimum and small steps if we converge towards this minimum.\n",
    "\n",
    "If we combine all theses ideas\n",
    "- momentum\n",
    "- random choice on which data to calculate gradient (stochastic, batch)\n",
    "- adaptive step size\n",
    "we get into the world of the currently used robust GD algorithms, such as ADAM."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "651b1eff",
   "metadata": {},
   "source": [
    "## Copyright\n",
    "\n",
    "- the notebooks are provided as [Open Educational Resources](https://en.wikipedia.org/wiki/Open_educational_resources)\n",
    "- feel free to use the notebooks for your own purposes\n",
    "- the text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/)\n",
    "- the code of the IPython examples is licensed under the [MIT license](https://opensource.org/licenses/MIT)\n",
    "- please attribute the work as follows: *Frank Schultz, Data Driven Audio Signal Processing - A Tutorial Featuring Computational Examples, University of Rostock* ideally with relevant file(s), github URL https://github.com/spatialaudio/data-driven-audio-signal-processing-exercise, commit number and/or version tag, year."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "myddasp",
   "language": "python",
   "name": "myddasp"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}