gradient_descent_momentum.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": "markdown",
   "id": "dcd5fabd",
   "metadata": {},
   "source": [
    "# Gradient Descent with Momentum"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbb4aa25",
   "metadata": {},
   "source": [
    "## Analytical Loss Function in 2D\n",
    "\n",
    "Suppose, that the (made up) loss function of a model with two parameters $\\beta_1$ and $\\beta_2$ is analytically given as\n",
    "\n",
    "$$\\mathcal{L}(\\beta_1, \\beta_2) = (\\beta_1 - 2)^2 + (\\beta_2 - 1)^4 - (\\beta_2 -1)^2$$\n",
    "\n",
    "In this toy example there is no data dependency involved, which is not how things work in practice, but it is  good to understand the essence of finding a minimum numerically.\n",
    "\n",
    "In order to find **potential minima**, and thereby the **optimum model parameters** $\\hat{\\beta_1}$ and $\\hat{\\beta_2}$, we need to solve gradient for zero\n",
    "\n",
    "$$\\nabla \\mathcal{L} = \n",
    "\\begin{bmatrix}\n",
    "\\frac{\\partial \\mathcal{L}}{\\partial \\beta_1}\\\\\n",
    "\\frac{\\partial \\mathcal{L}}{\\partial \\beta_2}\n",
    "\\end{bmatrix}=\n",
    "\\mathbf{0}$$\n",
    "\n",
    "The required partial derivatives of first order are \n",
    "\n",
    "$$\\frac{\\partial \\mathcal{L}}{\\partial \\beta_1} = 2 (\\beta_1 - 2)^1$$\n",
    "\n",
    "$$\\frac{\\partial \\mathcal{L}}{\\partial \\beta_2} = 4 (\\beta_2 - 1)^3 - 2(\\beta_2 -1)^1$$\n",
    "\n",
    "A check with the Hessian of $\\mathcal{L}(\\beta_1, \\beta_2)$ yields whether we deal with a minimum, maximum, saddle point or neither of them for each of the zero gradient conditions.\n",
    "\n",
    "We get **first minimum** at\n",
    "$$\\beta_{1,min1} = 2\\qquad \\beta_{2,min1} = 1+\\frac{1}{\\sqrt{2}}$$\n",
    "\n",
    "We get **second minimum** at\n",
    "$$\\beta_{1,min2} = 2\\qquad \\beta_{2,min2} = 1-\\frac{1}{\\sqrt{2}}$$\n",
    "\n",
    "Both minima yield the same function value\n",
    "$$\\mathcal{L}(\\beta_{1,min}, \\beta_{2,min}) = -\\frac{1}{4},$$\n",
    "so we deal actually with **two optimum models**, as there is **no global minimum** with only one lowest function value.\n",
    "\n",
    "We have **one saddle point**, that actually separates the two minima at\n",
    "$$\\beta_{1,saddle} = 2\\qquad \\beta_{2,saddle} = 1$$\n",
    "\n",
    "with function value\n",
    "$$\\mathcal{L}(\\beta_{1,saddle}, \\beta_{2,saddle}) = 0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5752a87",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits import mplot3d\n",
    "from scipy.signal import freqz, tf2zpk\n",
    "\n",
    "matplotlib_widget_flag = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8b0d840",
   "metadata": {},
   "outputs": [],
   "source": [
    "# analytical solutions from above\n",
    "minimum1, fminimum1 = np.array([2, 1 + 1 / np.sqrt(2)]), -1 / 4\n",
    "minimum2, fminimum2 = np.array([2, 1 - 1 / np.sqrt(2)]), -1 / 4\n",
    "saddle, fsaddle = np.array([2, 1]), 0\n",
    "\n",
    "\n",
    "def get_gradient(beta):\n",
    "    beta_gradient = np.array(\n",
    "        [2 * (beta[0] - 2), 4 * (beta[1] - 1) ** 3 - 2 * (beta[1] - 1)]\n",
    "    )\n",
    "    return beta_gradient"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67bc3812",
   "metadata": {},
   "source": [
    "## Plot the Loss Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f54b6d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x, y):\n",
    "    # our loss function from above\n",
    "    # only with nicer to read variables x,y\n",
    "    # instead of beta1, beta2\n",
    "    return (x - 2) ** 2 + (y - 1) ** 4 - (y - 1) ** 2\n",
    "\n",
    "\n",
    "x, y = np.linspace(0, 4, 2**7), np.linspace(-1 / 2, 5 / 2, 2**7)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "Z = f(X, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a13e8731",
   "metadata": {},
   "outputs": [],
   "source": [
    "col_min, col_max, no_col = -1 / 2, 5, 12\n",
    "col_tick = np.linspace(col_min, col_max, no_col, endpoint=True)\n",
    "cmap = plt.cm.magma_r\n",
    "norm = mpl.colors.BoundaryNorm(col_tick, cmap.N)\n",
    "\n",
    "if matplotlib_widget_flag:\n",
    "    %matplotlib widget\n",
    "fig = plt.figure()\n",
    "ax = plt.axes(projection=\"3d\")\n",
    "c = ax.plot_surface(\n",
    "    X, Y, Z, rstride=1, cstride=1, cmap=cmap, norm=norm, edgecolor=\"none\"\n",
    ")\n",
    "ax.plot(minimum1[0], minimum1[1], fminimum1, \"C0o\")\n",
    "ax.plot(minimum2[0], minimum2[1], fminimum2, \"C0o\")\n",
    "ax.plot(saddle[0], saddle[1], fsaddle, \"C0o\")\n",
    "cbar = fig.colorbar(\n",
    "    c, ax=ax, ticks=col_tick[:: no_col // 10], label=r\"$\\mathcal{L}$\"\n",
    ")\n",
    "ax.set_xlim(x[0], x[-1])\n",
    "ax.set_ylim(y[0], y[-1])\n",
    "ax.set_zlim(-1 / 4, 5)\n",
    "ax.set_xlabel(r\"$\\beta_1$\")\n",
    "ax.set_ylabel(r\"$\\beta_2$\")\n",
    "ax.set_zlabel(r\"$\\mathcal{L}$\")\n",
    "ax.view_init(elev=60, azim=-40)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cb567ee",
   "metadata": {},
   "source": [
    "## Gradient Descent Update Rule \n",
    "We could (and in ML problems we often need) to find the minima numerically.\n",
    "The most straightforward and simple numerical solver is the so called **gradient descent** (GD), a first order method.\n",
    "\n",
    "It uses the (analytically known) gradient $\\nabla\\mathcal{L}$, evaluates it for an actual $\\beta_{actual}$ and updates subsequently into direction of negative gradient, i.e. the (or rather a?!?) minimum that we want to find.\n",
    "\n",
    "This iterative procedure can be written as\n",
    "$$(1):\\quad \\beta_{new} = \\beta_{actual} - \\mathrm{step size} \\cdot \\nabla\\mathcal{L}\\bigg|_{\\beta_{actual}}\\quad(2):\\quad\\beta_{new} \\rightarrow \\beta_{actual}\\quad(3): \\mathrm{go to}\\,(1)$$\n",
    "repeated until we hopefully converged to the $\\beta$ that represents the minimum.\n",
    "\n",
    "In practice GD is not often used, as it is not very robust and many things can go wrong. Let us check this with some illustrative examples below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "813e447d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_gradient_descent(beta, step_size):\n",
    "    beta_gradient = get_gradient(beta)\n",
    "    beta -= beta_gradient * step_size\n",
    "    return beta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63fa751f",
   "metadata": {},
   "source": [
    "## Essence of Gradient Descent with Momentum\n",
    "\n",
    "One potential problem of straightforward GD is the potentially high fluctuations in the gradient update.\n",
    "The idea behind adding certain *momentum* onto the update rule is to increase the inertia of the GD path and thus decrease the fluctuations.\n",
    "We could do this with a simple lowpass filter applying onto the calculated gradient signal, such that fast direction changes are low pass filtered and thereby cannot change so fast (as we cut high frequencies in the spectrum)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8f706a9",
   "metadata": {},
   "source": [
    "A simple one-pole lowpass filter as known from a basic DSP course does the job for the upcoming toy examples.\n",
    "\n",
    "This one-pole filter needs one storage element, which we could call a intermediate variable $z$. The momentum $0\\leq m \\leq 1$ is controlled via the cut-off frequency of the lowpass.\n",
    "\n",
    "The iterative procedure can be written as\n",
    "\\begin{align}\n",
    "&(1):\\quad \\beta_{new} = \\beta_{actual} - \\mathrm{step size} \\cdot z_{actual}\\quad\\\\\n",
    "&(2):\\quad z_{new} = (1-m)\\cdot\\nabla\\mathcal{L}\\bigg|_{\\beta_{new}} + m \\cdot z_{actual}\\\\\n",
    "&(3):\\quad\\beta_{new} \\rightarrow \\beta_{actual}\\quad z_{new} \\rightarrow {z}_{actual}\\\\\n",
    "&(4):\\quad\\mathrm{go to}\\,(1)\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "repeated until we hopefully converged to the $\\beta$ that represents the minimum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5646a5a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_gradient_descent_momentum(beta, z, step_size, momentum_coeff):\n",
    "    beta = beta - step_size * z\n",
    "    beta_gradient = get_gradient(beta)\n",
    "    # this is a one pole filter in DSP\n",
    "    z = (1 - momentum_coeff) * beta_gradient + momentum_coeff * z\n",
    "    return beta, z"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36e1783a",
   "metadata": {},
   "source": [
    "### Lowpass Analysis\n",
    "\n",
    "The line (2) in the iterative procedure belongs to the code line \n",
    "`z = (1-momentum_coeff) * beta_gradient + momentum_coeff * z`,\n",
    "which is a recursive equation, i.e. a difference equation. From DSP we know how to analyze such an LTI system.\n",
    "With \n",
    "$$b_0 = (1-\\mathrm{momentum\\_coeff}) = 1-m$$ \n",
    "and\n",
    "$$a_1 = -\\mathrm{momentum\\_coeff} = - m$$\n",
    "we rewrite this difference equation in typical DSP fashion\n",
    "\n",
    "$$y[k] = b_0 x[k] - a_1 y[k-1]$$\n",
    "\n",
    "The z-transform\n",
    "\n",
    "$$Y(z) = b_0 X(z) - a_1 z^{-1} Y(z)$$\n",
    "\n",
    "$$Y(z) + a_1 z^{-1} Y(z) = b_0 X(z)$$\n",
    "\n",
    "$$(1 + a_1 z^{-1} ) Y(z) = b_0 X(z)$$\n",
    "\n",
    "yields the transfer function\n",
    "\n",
    "$$H(z) = \\frac{Y(z)}{X(z)} = \\frac{b_0}{1 + a_1 z^{-1}} = \\frac{b_0 z}{z + a_1}$$\n",
    "\n",
    "with a zero in origin $z_0 = 0$ and one real pole (on the $\\Re(z)$ axis) at $z_\\infty = -a_1=+\\mathrm{momentum\\_coeff}$ \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46d1fc28",
   "metadata": {},
   "outputs": [],
   "source": [
    "momentum_coeff = 0.25  # arbitrary number\n",
    "# the higher the lower the cutoff frequency\n",
    "# what does this mean in terms of the momentum for the GD?!\n",
    "\n",
    "b0 = 1 - momentum_coeff\n",
    "a0, a1 = 1, -momentum_coeff\n",
    "\n",
    "b, a = np.array([b0, 0]), np.array([a0, a1])\n",
    "w, h = freqz(b, a)\n",
    "plt.plot(w / np.pi, 20 * np.log10(np.abs(h)))\n",
    "plt.xlim(0, 1)\n",
    "plt.ylim(-10, 1)\n",
    "plt.xlabel(r\"$\\Omega / \\pi$\")\n",
    "plt.ylabel(\"level in dB\")\n",
    "plt.title(\"lowpass filter acting on gradient -> GD with higher inertia\")\n",
    "plt.grid(True)\n",
    "\n",
    "z, p, k = tf2zpk(b, a)\n",
    "print(\"zero and pole as analytically derived:\", z, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8390fbc",
   "metadata": {},
   "source": [
    "### Lowpass Filter Alternative\n",
    "\n",
    "We might also consider other lowpass filter characteristics. The first-order system from above has one more degree of freedom, namely the coefficient $b_1\\neq 0$ for $x[k-1]$.\n",
    "The typical DSP analysis shows\n",
    "\n",
    "$$y[k] = b_0 x[k] + b_1 x[k-1] - a_1 y[k-1]$$\n",
    "\n",
    "$$Y(z) = b_0 X(z) + b_1 X(z) z^{-1} - a_1 Y(z) z^{-1}$$\n",
    "\n",
    "$$Y(z) + a_1 Y(z) z^{-1} = b_0 X(z) + b_1 X(z) z^{-1}$$\n",
    "\n",
    "$$Y(z) (1 + a_1  z^{-1}) = X(z) (b_0 + b_1 z^{-1})$$\n",
    "\n",
    "$$H(z) = \\frac{b_0 + b_1 z^{-1}}{1 + a_1  z^{-1}} = \\frac{b_0 z + b_1 }{z + a_1}$$\n",
    "\n",
    "Ensuring 0 dB gain at DC needs the condition\n",
    "$$H(z=1) = \\frac{b_0 + b_1 }{1+a_1} = 1$$\n",
    "which can be used to scale the coefficients accordingly.\n",
    "\n",
    "See a toy example below, where we aim for level $-\\infty$ dB at half sampling frequency, i.e. a zero at $z_0=-1$. This is achieved for $b_0=b_1=1$ (and all scaled versions). The cutoff-frequency is controlled by the momentum coeff, which shifts the real pole $z_\\infty$ of the filter just as in the first lowpass filter option.\n",
    "\n",
    "A nice homework would be to implement this filter characteristics into `my_gradient_descent_momentum(beta, z, step_size, momentum_coeff)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec0473bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "momentum_coeff = 0.25  # arbitrary number\n",
    "# the higher the lower the cutoff frequency\n",
    "# what does this mean in terms of the momentum for the GD?!\n",
    "\n",
    "b0, b1 = 1.0, 1.0\n",
    "a0, a1 = 1.0, -momentum_coeff\n",
    "k = np.abs((b0 + b1) / (1 + a1))\n",
    "print(1 / k)\n",
    "b, a = 1 / k * np.array([b0, b1]), np.array([a0, a1])\n",
    "print(b)\n",
    "\n",
    "w, h = freqz(b, a)\n",
    "z, p, k = tf2zpk(b, a)\n",
    "print(\"zero, pole, k as analytically derived:\", z, p, k)\n",
    "\n",
    "plt.figure()\n",
    "plt.subplot(211)\n",
    "plt.plot(w / np.pi, 20 * np.log10(np.abs(h)))\n",
    "plt.xlim(0, 1)\n",
    "plt.ylim(-40, 1)\n",
    "plt.xlabel(r\"$\\Omega / \\pi$\")\n",
    "plt.ylabel(\"level in dB\")\n",
    "plt.title(\"lowpass alternative acting on gradient -> GD with higher inertia\")\n",
    "plt.grid(True)\n",
    "plt.subplot(212)\n",
    "plt.plot(w / np.pi, np.angle(h))\n",
    "plt.xlim(0, 1)\n",
    "plt.xlabel(r\"$\\Omega / \\pi$\")\n",
    "plt.ylabel(\"pahse in rad\")\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd0e4077",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_contour_of_my_gradient_descent():\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.contour(X, Y, Z, cmap=\"magma_r\")\n",
    "    ax.plot(beta_path_gd[0], beta_path_gd[1], \"C0d-\", lw=0.5, ms=4, label=\"GD\")\n",
    "    ax.plot(\n",
    "        beta_path_gd_momentum[0],\n",
    "        beta_path_gd_momentum[1],\n",
    "        \"C8o-\",\n",
    "        lw=0.5,\n",
    "        ms=4,\n",
    "        label=\"GD with momentum\",\n",
    "    )\n",
    "    ax.plot(minimum1[0], minimum1[1], \"kx\", ms=10)\n",
    "    ax.plot(minimum2[0], minimum2[1], \"kx\", ms=10)\n",
    "    ax.plot(saddle[0], saddle[1], \"kx\", ms=10)\n",
    "    ax.axis(\"equal\")\n",
    "    ax.set_xlim(0, 4)\n",
    "    ax.set_ylim(-0.25, 2.5)\n",
    "    ax.set_xlabel(r\"$\\beta_1$\")\n",
    "    ax.set_ylabel(r\"$\\beta_2$\")\n",
    "    ax.legend()\n",
    "    ax.grid(\"True\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16da4461",
   "metadata": {},
   "source": [
    "### Gradient Descent vs. Gradient Descent Momentum\n",
    "\n",
    "GD and GD with momentum both find minimum 1. Higher inertia leads another path for momentum GD, which therefore needs slightly more updates to converge. This problem is not critical for GD."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf675363",
   "metadata": {},
   "outputs": [],
   "source": [
    "steps = 2**5\n",
    "beta_gd = np.array([4, 1.075])\n",
    "step_size = 1 / 5\n",
    "\n",
    "z = np.array([0, 0])\n",
    "momentum_coeff = 0.33\n",
    "\n",
    "beta_gd_momentum = np.copy(beta_gd)\n",
    "beta_path_gd = np.zeros((2, steps + 1))\n",
    "beta_path_gd[:, 0] = beta_gd\n",
    "beta_path_gd_momentum = np.zeros((2, steps + 1))\n",
    "beta_path_gd_momentum[:, 0] = beta_gd_momentum\n",
    "\n",
    "for step in range(steps):\n",
    "    beta_gd = my_gradient_descent(beta_gd, step_size)\n",
    "    beta_path_gd[:, step + 1] = beta_gd\n",
    "\n",
    "    beta_gd_momentum, z = my_gradient_descent_momentum(\n",
    "        beta_gd_momentum, z, step_size, momentum_coeff\n",
    "    )\n",
    "    beta_path_gd_momentum[:, step + 1] = beta_gd_momentum\n",
    "\n",
    "# and plot\n",
    "plot_contour_of_my_gradient_descent()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57c29171",
   "metadata": {},
   "source": [
    "### Gradient Descent to Minimum 2 vs. Gradient Descent Momentum to Minimum 1\n",
    "\n",
    "GD finds its path to minimum 2 straightforward, as the GD parameters are well chosen.\n",
    "Momentum GD moves away from minimum 2 into the upper half and hence towards minimum 1. This is only possible due to the higher inertia imposed onto the gradient signal.\n",
    "This higher inertia then however leads to a 'crazy' GD path until minimum 1 is finally found."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab382881",
   "metadata": {},
   "outputs": [],
   "source": [
    "steps = 2**10\n",
    "beta_gd = np.array([0.0, -0.5])\n",
    "step_size = 1 / 50\n",
    "\n",
    "z = np.array([0, 0])\n",
    "momentum_coeff = 0.99\n",
    "\n",
    "beta_gd_momentum = np.copy(beta_gd)\n",
    "beta_path_gd = np.zeros((2, steps + 1))\n",
    "beta_path_gd[:, 0] = beta_gd\n",
    "beta_path_gd_momentum = np.zeros((2, steps + 1))\n",
    "beta_path_gd_momentum[:, 0] = beta_gd_momentum\n",
    "\n",
    "for step in range(steps):\n",
    "    beta_gd = my_gradient_descent(beta_gd, step_size)\n",
    "    beta_path_gd[:, step + 1] = beta_gd\n",
    "\n",
    "    beta_gd_momentum, z = my_gradient_descent_momentum(\n",
    "        beta_gd_momentum, z, step_size, momentum_coeff\n",
    "    )\n",
    "    beta_path_gd_momentum[:, step + 1] = beta_gd_momentum\n",
    "\n",
    "# and plot\n",
    "plot_contour_of_my_gradient_descent()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5132abe9",
   "metadata": {},
   "source": [
    "### Gradient Descent vs. Gradient Descent Momentum\n",
    "\n",
    "Both GDs do not converge when they start on the saddle point ridge.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a469708",
   "metadata": {},
   "outputs": [],
   "source": [
    "steps = 2**6\n",
    "beta_gd = np.array([4.0, 1.0])\n",
    "step_size = 1 / 10\n",
    "\n",
    "z = np.array([0, 0])\n",
    "momentum_coeff = 0.95\n",
    "\n",
    "beta_gd_momentum = np.copy(beta_gd)\n",
    "beta_path_gd = np.zeros((2, steps + 1))\n",
    "beta_path_gd[:, 0] = beta_gd\n",
    "beta_path_gd_momentum = np.zeros((2, steps + 1))\n",
    "beta_path_gd_momentum[:, 0] = beta_gd_momentum\n",
    "\n",
    "for step in range(steps):\n",
    "    beta_gd = my_gradient_descent(beta_gd, step_size)\n",
    "    beta_path_gd[:, step + 1] = beta_gd\n",
    "\n",
    "    beta_gd_momentum, z = my_gradient_descent_momentum(\n",
    "        beta_gd_momentum, z, step_size, momentum_coeff\n",
    "    )\n",
    "    beta_path_gd_momentum[:, step + 1] = beta_gd_momentum\n",
    "\n",
    "# and plot\n",
    "plot_contour_of_my_gradient_descent()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7522f5e2",
   "metadata": {},
   "source": [
    "### Exploding Gradient Descent vs. Stable Momentum Gradient Descent\n",
    "\n",
    "Here, the higher inertia imposed with momentum helps that the GD path does not explode.\n",
    "\n",
    "The straightforward GD does not converge, but rather explodes. The momentum GD finds its way to minimum 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97d002c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "steps = 2**4\n",
    "beta_gd = np.array([4.0, 2.5])\n",
    "step_size = 0.29\n",
    "z = np.array([0, 0])\n",
    "momentum_coeff = 0.08  # 0.03 vs. 0.09\n",
    "\n",
    "beta_gd_momentum = np.copy(beta_gd)\n",
    "beta_path_gd = np.zeros((2, steps + 1))\n",
    "beta_path_gd[:, 0] = beta_gd\n",
    "beta_path_gd_momentum = np.zeros((2, steps + 1))\n",
    "beta_path_gd_momentum[:, 0] = beta_gd_momentum\n",
    "\n",
    "for step in range(steps):\n",
    "    beta_gd = my_gradient_descent(beta_gd, step_size)\n",
    "    beta_path_gd[:, step + 1] = beta_gd\n",
    "\n",
    "    beta_gd_momentum, z = my_gradient_descent_momentum(\n",
    "        beta_gd_momentum, z, step_size, momentum_coeff\n",
    "    )\n",
    "    beta_path_gd_momentum[:, step + 1] = beta_gd_momentum\n",
    "\n",
    "# and plot\n",
    "plot_contour_of_my_gradient_descent()"
   ]
  },
  {
   "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
}