exercise07_linear_regression_LS_vs_SVD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "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",
    "# Tutorial Selected Topics in Audio Signal Processing\n",
    "\n",
    "Master Course #24512\n",
    "\n",
    "- lecture: https://github.com/spatialaudio/selected-topics-in-audio-signal-processing-lecture\n",
    "- tutorial: https://github.com/spatialaudio/selected-topics-in-audio-signal-processing-exercises\n",
    "\n",
    "Feel free to contact lecturer frank.schultz@uni-rostock.de"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Least Squares Error Regression vs. Orthogonal Regression (via SVD)\n",
    "\n",
    "using 2D toy example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "np.random.seed(10)\n",
    "N = 2**10\n",
    "\n",
    "mean = [0, 0]\n",
    "\n",
    "cov = [[1, 0.75], [0.75, 1]]\n",
    "# cov = [[1, 0.2], [0.2, 0.25]]\n",
    "# cov = [[1, 0.2], [0.2, 0.75]]\n",
    "# cov = [[1, 0.9999999], [0.9999999, 1]]\n",
    "# cov = [[1, 0.0001], [0.0001, 1]]\n",
    "\n",
    "data = np.random.multivariate_normal(mean, cov, N)\n",
    "# make mean free for fair comparison of both approaches\n",
    "data = data - np.mean(data)\n",
    "print(\"dim of data\", data.shape)\n",
    "\n",
    "# Least Squares Regression X beta = meas, we want to find beta\n",
    "X = np.array([np.ones(N), data[:, 0]]).T\n",
    "print(\"dim of matrix A for LS:\", X.shape)\n",
    "meas = data[:, 1]\n",
    "# analytical LS solution using left inverse\n",
    "beta = (np.linalg.inv((X.T @ X)) @ X.T) @ meas\n",
    "print(\"intercept, slope for LS\", beta)\n",
    "# unit vector along the slope of LS reg line\n",
    "v_ls = [1, beta[1]]\n",
    "v_ls = v_ls / np.linalg.norm(v_ls, 2)\n",
    "\n",
    "# SVD Regression\n",
    "u, s, vh = np.linalg.svd(data)\n",
    "# right sing vector with largest sing value indicates the slope of regression line\n",
    "v_max_sing = vh[:, 0].T\n",
    "m_svd = v_max_sing[1] / v_max_sing[0]\n",
    "n_svd = 0  # by definition of spanning U,V spaces\n",
    "print(\"intercept, slope for svd\", [n_svd, m_svd])\n",
    "\n",
    "# plot\n",
    "x_predict = np.linspace(-10, 10, 2)\n",
    "plt.figure(figsize=(6, 6))\n",
    "plt.plot(data[:, 0], data[:, 1], \"C0o\", ms=2, label=\"data\")\n",
    "plt.plot(\n",
    "    x_predict,\n",
    "    beta[0] + beta[1] * x_predict,\n",
    "    \"C1-\",\n",
    "    label=\"LS regression line fit\",\n",
    ")\n",
    "# plt.plot(v_max_sing[0] * x_predict, v_max_sing[1] * x_predict,\n",
    "#         'C3-', label='SVD regression line fit')  # ==\n",
    "plt.plot(\n",
    "    x_predict, n_svd + m_svd * x_predict, \"C3-\", label=\"SVD regression line fit\"\n",
    ")\n",
    "plt.axis(\"square\")\n",
    "plt.xlim(-4, 4)\n",
    "plt.ylim(-4, 4)\n",
    "plt.xticks(np.arange(-4, 5))\n",
    "plt.yticks(np.arange(-4, 5))\n",
    "plt.xlabel(\"data[:,0]\")\n",
    "plt.ylabel(\"data[:,1]\")\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "\n",
    "# we hope to find at least one data point:\n",
    "tmp1 = data[data[:, 0] > 3, :]\n",
    "tmp1 = tmp1[0]  # get first data point of it\n",
    "plt.plot(tmp1[0], tmp1[1], \"C0o\", ms=7)  # draw it large\n",
    "\n",
    "# connection line between this single data point and the\n",
    "# regression lines\n",
    "# SVD line regression: data point projection -> orthogonal\n",
    "tmp2 = np.inner(tmp1, v_max_sing) * v_max_sing\n",
    "plt.plot([tmp1[0], tmp2[0]], [tmp1[1], tmp2[1]], \"C3-.\", lw=2)\n",
    "# LS line regression: x of data point == x of regression line\n",
    "tmp2 = beta[0] + beta[1] * tmp1[0]\n",
    "plt.plot([tmp1[0], tmp1[0]], [tmp1[1], tmp2], \"C1-.\", lw=2)\n",
    "\n",
    "# check the minimization criteria:\n",
    "# for squared sums\n",
    "ss_ls, ss_svd = 0, 0\n",
    "for n in range(data.shape[0]):\n",
    "    x_data, y_data = data[n, 0], data[n, 1]\n",
    "    # get squared sums for both approaches\n",
    "    ss_ls += (y_data - beta[0] + beta[1] * x_data) ** 2\n",
    "    ss_svd += (y_data - n_svd + m_svd * x_data) ** 2\n",
    "    # compare data points vs. prediction\n",
    "    # print(y_data, x[0] + x[1]*x_data)  # LS reg\n",
    "    # print(y_data, n_svd + m_svd*x_data+)  # SVD reg\n",
    "ss_ls *= 1 / data.shape[0]\n",
    "ss_svd *= 1 / data.shape[0]\n",
    "print(\"\\nsquared sums:\")\n",
    "print(\"ss_ls =\", ss_ls, \"< ss_svd =\", ss_svd)\n",
    "\n",
    "# for squared orth dist\n",
    "sod_ls, sod_svd = 0, 0\n",
    "for n in range(data.shape[0]):\n",
    "    # get actual data point\n",
    "    tmp1 = data[n, :]\n",
    "\n",
    "    # for SVD\n",
    "    # project down to v => length, create weighted v\n",
    "    tmp2 = np.inner(tmp1, v_max_sing) * v_max_sing\n",
    "    # squared distance between tmp1 and tmp2\n",
    "    sod_svd += np.linalg.norm(tmp2 - tmp1, 2) ** 2\n",
    "\n",
    "    # for LS\n",
    "    # offset data by intercept of LS regression line\n",
    "    tmp1[1] += beta[0]\n",
    "    # then we can use straightforward projection\n",
    "    tmp2 = np.inner(tmp1, v_ls) * v_ls\n",
    "    # squared distance between tmp1 and tmp2\n",
    "    sod_ls += np.linalg.norm(tmp2 - tmp1, 2) ** 2\n",
    "\n",
    "sod_ls *= 1 / data.shape[0]\n",
    "sod_svd *= 1 / data.shape[0]\n",
    "print(\"\\nsquared orthogonal distances:\")\n",
    "print(\"sod_ls =\", sod_ls, \"> sod_svd =\", sod_svd)"
   ]
  },
  {
   "cell_type": "markdown",
   "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": {
  "interpreter": {
   "hash": "1743232f157dd0954c61aae30535e75a2972519a625c7e796bafe0cd9a07bf7e"
  },
  "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": 4
}