lcurve.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ea6d9c2b",
   "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\n",
    "\n",
    "# Check the L-Curve for Ridge Regression\n",
    "with a matrix $\\textbf{X}$ - full column rank, but rather high condition number.\n",
    "\n",
    "Using toy data, we calculate and plot the L-curve and find an optimum regularization parameter $\\lambda_{opt}$, i.e. where L-curve has maximum curvature. For this $\\lambda_{opt}$ we get a suitable trade-off between prediction error and magnitude of model parameter vector. L-curve is typically plotted log-log because of rather large range of numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9c24554",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy.linalg import inv, matrix_rank, cond\n",
    "from scipy.linalg import svd, diagsvd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e92f286e",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "001db254",
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = np.array([0, 0])\n",
    "tmp = 1 - 1e-4\n",
    "cov = [[1, tmp], [tmp, 1]]\n",
    "M, N = 2**7, mu.shape[0]\n",
    "\n",
    "X = rng.multivariate_normal(mu, cov, size=M, method='cholesky')\n",
    "[U, s, Vh] = svd(X)\n",
    "Uh = U.T.conj()\n",
    "V = Vh.T.conj()\n",
    "S = diagsvd(s, M, N)\n",
    "I = np.eye(N)\n",
    "\n",
    "n = rng.standard_normal(size=M)\n",
    "\n",
    "y = (2 * X[:, 0] + X[:, 1] + n)[:, None]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "199f2ee6",
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix_rank(X), s.shape[0], cond(X), s[0] / s[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee083563",
   "metadata": {},
   "outputs": [],
   "source": [
    "lmb = np.logspace(-2.2, 1.9, 2**7)\n",
    "lcurve = np.zeros([2, lmb.shape[0]])\n",
    "\n",
    "for cnt, val in enumerate(lmb):\n",
    "    X_reg_left_inv = V @ inv(S.T @ S + val * I) @ S.T  @ Uh\n",
    "    theta_hat = X_reg_left_inv @ y\n",
    "    tmp = y - X @ theta_hat\n",
    "    lcurve[0, cnt] = tmp.T.conj() @ tmp\n",
    "    lcurve[1, cnt] = theta_hat.T.conj() @  theta_hat\n",
    "\n",
    "curvature = np.gradient(np.gradient(lcurve[1, :]))\n",
    "lmb_opt_idx = np.argmin(np.abs(curvature))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7a92ebe",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.loglog(lcurve[0, :], lcurve[1, :], 'o:', ms=2)\n",
    "plt.loglog(lcurve[0, lmb_opt_idx], lcurve[1, lmb_opt_idx], 'C3o')\n",
    "plt.text(lcurve[0, 0], lcurve[1, 0], r'$\\lambda\\rightarrow 0$')\n",
    "plt.text(lcurve[0, -1], lcurve[1, -1], r'$\\lambda\\rightarrow \\infty$')\n",
    "plt.text(lcurve[0, lmb_opt_idx], lcurve[1, lmb_opt_idx],\n",
    "         r'$\\lambda_\\mathrm{opt}$')\n",
    "plt.xlabel(r'$||\\bf{y} - \\bf{X \\hat{\\theta}}||_2^2$')\n",
    "plt.ylabel(r'$||\\bf{\\hat{\\theta}}||_2^2$')\n",
    "plt.title(r'L-Curve for $\\hat{\\bf{\\theta}}(\\lambda)$')\n",
    "plt.xlim(100, 200)\n",
    "plt.ylim(1, 100)\n",
    "plt.grid(True, which='both', ls=':', lw=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29f17512",
   "metadata": {},
   "outputs": [],
   "source": [
    "# l-curve for tikzpicture in ddasp_exercise_slides.tex\n",
    "print('coordinates {', sep='', end='')\n",
    "for cnt, val in enumerate(lmb):\n",
    "    print('(', np.log10(lcurve[0, cnt]), ',', np.log10(\n",
    "        lcurve[1, cnt]), ')', sep='', end='')\n",
    "print('};', sep='', end='')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c3277b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# l-curve values of opt lambda for tikzpicture in ddasp_exercise_slides.tex\n",
    "print('coordinates {', sep='', end='')\n",
    "print('(', np.log10(lcurve[0, lmb_opt_idx]), ',', np.log10(\n",
    "    lcurve[1, lmb_opt_idx]), ')', sep='', end='')\n",
    "print('};', sep='', end='')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2470b8ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# l-curve values for small lambda for tikzpicture in ddasp_exercise_slides.tex\n",
    "print('coordinates {', sep='', end='')\n",
    "print('(', np.log10(lcurve[0, 0]), ',', np.log10(\n",
    "    lcurve[1, 0]), ')', sep='', end='')\n",
    "print('};', sep='', end='')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4c60168",
   "metadata": {},
   "outputs": [],
   "source": [
    "# l-curve values for large lambda for tikzpicture in ddasp_exercise_slides.tex\n",
    "print('coordinates {', sep='', end='')\n",
    "print('(', np.log10(lcurve[0, -1]), ',', np.log10(\n",
    "    lcurve[1, -1]), ')', sep='', end='')\n",
    "print('};', sep='', end='')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0897f31e",
   "metadata": {},
   "source": [
    "## Copyright\n",
    "\n",
    "- the notebooks are provided as [Open Educational Resources](https://en.wikipedia.org/wiki/Open_educational_resources)\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",
    "- feel free to use the notebooks for your own purposes\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
}