{
"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
}