exercise07_left_inverse_SVD_QR.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "489c7bcd-d5b9-4a10-beb9-e58120f4758b",
   "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": "201797c2-3cc5-4cce-b1b9-684286f3b759",
   "metadata": {},
   "source": [
    "# Exercise 7: Least Squares Error Solution / Left Inverse in SVD and QR Domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59a782b9-b4d8-482e-b834-31522998e4c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.linalg import svd, diagsvd, qr, inv, pinv, norm\n",
    "from numpy.linalg import matrix_rank\n",
    "\n",
    "np.set_printoptions(precision=2, floatmode=\"fixed\", suppress=True)\n",
    "\n",
    "matplotlib_widget_flag = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5bb19c8",
   "metadata": {},
   "source": [
    "For the given full-column rank feature matrix $\\mathbf{X}$ (tall/thin shape with independent columns with $M \\times N$ dimension, $M$ samples, $N$ features) and outcome vector ${\\beta}$ the linear set of equations\n",
    "\n",
    "$$\\mathbf{X} {\\beta} = \\mathbf{y}$$\n",
    "\n",
    "is to be solved for unknowns ${\\beta}$.\n",
    "We obviously cannot invert $\\mathbf{X}$, so we must find an optimum row space estimate $\\hat{\\beta}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c72192e3",
   "metadata": {},
   "source": [
    "### Least Squares Error Solution\n",
    "\n",
    "great material, strong recommendation:\n",
    "- Gilbert Strang (2020): \"Linear Algebra for Everyone\", Wellesley-Cambridge Press, Ch. 4.3\n",
    "- Gilbert Strang (2019): \"Linear Algebra and Learning from Data\", Wellesley-Cambridge Press, Ch. II.2\n",
    "\n",
    "We know for sure, that **pure row space** $\\hat{\\beta}$ maps to **pure column space** $\\hat{\\mathbf{y}}$ for the linear combination\n",
    "\n",
    "$$\\mathbf{X} \\hat{\\beta} = \\hat{\\mathbf{y}}$$\n",
    "\n",
    "The given (measured) outcome vector $\\mathbf{y}$ might not necessarily (in practical problems probably never!) live in pure column space of $\\mathbf{X}$.\n",
    "\n",
    "We therefore need an offset (error) vector to get there\n",
    "\n",
    "$$\\mathbf{X} \\hat{\\beta} + \\mathbf{e} = \\mathbf{y}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3c569bd",
   "metadata": {},
   "source": [
    "We want to find the $\\hat{\\beta}$ which yields the **smallest** $||\\mathbf{e}||_2^2$ or equivalently $||\\mathbf{e}||_2$.\n",
    "This is our optimization criterion, known as **least squares error**, because the sum of squares $||\\mathbf{e}||_2^2 = e_1^2 + e_2^2 + ... + e_M^2$ should be minimized.\n",
    "\n",
    "So, thinking in vector addition\n",
    "\n",
    "$$\\mathbf{y} = \\hat{\\mathbf{y}} + \\mathbf{e} \\rightarrow \\mathbf{e} = \\mathbf{y} - \\hat{\\mathbf{y}}$$\n",
    "\n",
    "we can geometrically figure (we should imagine and/or draw this is in 2D / 3D) that the smallest $||\\mathbf{e}||_2^2$ is achieved when we span a **right-angled triangle** using $\\mathbf{y}$ as hypotenuse."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43222278",
   "metadata": {},
   "source": [
    "Therefore, $\\hat{\\mathbf{y}} \\perp \\mathbf{e}$.\n",
    "Recall, that $\\hat{\\mathbf{y}}$ lives in column space of $\\mathbf{X}$, that's where we started from.\n",
    "Column space is orthogonal to left null space, so $\\hat{\\mathbf{y}} \\perp \\mathbf{e}$ tells us, that $\\mathbf{e}$ must live in left null space of $\\mathbf{X}$.\n",
    "This requirement can be formally written as $\\mathbf{X}^\\mathrm{H} \\mathbf{e} = \\mathbf{0}$ and then can be utilized as\n",
    "\n",
    "$$\\mathbf{X} \\hat{\\beta} + \\mathbf{e} = \\mathbf{y} \\qquad\\rightarrow\\qquad \\mathbf{X}^\\mathrm{H}\\mathbf{X} \\hat{\\beta} + \\mathbf{X}^\\mathrm{H}\\mathbf{e} = \\mathbf{X}^\\mathrm{H}\\mathbf{y} \\qquad\\rightarrow\\qquad \\mathbf{X}^\\mathrm{H}\\mathbf{X} \\hat{\\beta} = \\mathbf{X}^\\mathrm{H}\\mathbf{y}$$\n",
    "\n",
    "The last equation in the above line is known as normal equation(s)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c8b5af2",
   "metadata": {},
   "source": [
    "This can be solved using the left inverse of\n",
    "$$\\mathbf{X}^\\mathrm{H} \\mathbf{X}$$\n",
    "This matrix is full rank and therefore invertible.\n",
    "Hence, multiplying both sides of normal equation(s) from left\n",
    "\n",
    "$$(\\mathbf{X}^\\mathrm{H} \\mathbf{X})^{-1} (\\mathbf{X}^\\mathrm{H} \\mathbf{X}) \\hat{\\beta} = (\\mathbf{X}^\\mathrm{H} \\mathbf{X})^{-1} \\mathbf{X}^\\mathrm{H} \\mathbf{y}$$\n",
    "\n",
    "Since for left inverse $(\\mathbf{X}^\\mathrm{H} \\mathbf{X})^{-1} (\\mathbf{X}^\\mathrm{H} \\mathbf{X}) = \\mathbf{I}$ holds, we get the least-squares sense solution for $\\beta$ in the row space of $\\mathbf{X}$\n",
    "\n",
    "$$\\hat{\\beta} = (\\mathbf{X}^\\mathrm{H} \\mathbf{X})^{-1} \\mathbf{X}^\\mathrm{H} \\mathbf{y}$$\n",
    "\n",
    "using the **left inverse** of $\\mathbf{X}$ as\n",
    "\n",
    "$$\\mathbf{X}^{+L} = (\\mathbf{X}^\\mathrm{H} \\mathbf{X})^{-1} \\mathbf{X}^\\mathrm{H}$$\n",
    "\n",
    "So, optimum estimator for our unknowns (i.e. the model parameters)\n",
    "\n",
    "$$\\hat{\\beta} = \\mathbf{X}^{+L} \\mathbf{y}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86d4aa79",
   "metadata": {},
   "source": [
    "If we let matrix $\\mathbf{X}$ act on $\\hat{\\beta}$,\n",
    "\n",
    "$$\\mathbf{X} \\hat{\\beta} = \\mathbf{X}\\mathbf{X}^{+L} \\mathbf{y} = \\hat{\\mathbf{y}}$$\n",
    "\n",
    "we land in pure column space as initially started.\n",
    "Thus,\n",
    "\n",
    "$$\\mathbf{P}_\\mathrm{col space} = \\mathbf{X} \\mathbf{X}^{+L}$$\n",
    "\n",
    "is a **projection matrix**, that projects any $\\mathbf{y}$ into the column space of $\\mathbf{X}$ creating $\\hat{\\mathbf{y}}$, thus often referred to as the **hat matrix**.\n",
    "\n",
    "We can define another projection matrix\n",
    "\n",
    "$$\\mathbf{P}_\\mathrm{row space} = \\mathbf{X}^{+L} \\cdot \\mathbf{X},$$\n",
    "which maps any $\\beta$ to row space of $\\mathbf{X}$.\n",
    "From above we already know\n",
    "\n",
    "$$\\mathbf{P}_\\mathrm{row space} = (\\mathbf{X}^\\mathrm{H} \\mathbf{X})^{-1} \\mathbf{X}^\\mathrm{H} \\cdot \\mathbf{X} = \\mathbf{I}$$\n",
    "\n",
    "which makes sense, as the tall/thin, full column rank matrix does not span a null space, other than the $\\mathbf{0}$ vector.\n",
    "So, all possible $\\beta$ already live in the row space and must be projected exactly to itself. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e35aa0cd",
   "metadata": {},
   "source": [
    "### Least Squares Error Solution in SVD Domain\n",
    "\n",
    "great material, strong recommendation:\n",
    "- Gilbert Strang (2019): \"Linear Algebra and Learning from Data\", Wellesley-Cambridge Press, p.125ff\n",
    "\n",
    "The left inverse \n",
    "\n",
    "$$\\mathbf{X}^{+L} = (\\mathbf{X}^\\mathrm{H} \\mathbf{X})^{-1} \\mathbf{X}^\\mathrm{H}$$\n",
    "\n",
    "in terms of SVD\n",
    "\n",
    "$$\\mathbf{X}^{+L} = ((\\mathbf{U}\\mathbf{S}\\mathbf{V}^\\mathrm{H})^\\mathrm{H} \\mathbf{U}\\mathbf{S}\\mathbf{V}^\\mathrm{H})^{-1} (\\mathbf{U}\\mathbf{S}\\mathbf{V}^\\mathrm{H})^\\mathrm{H}$$\n",
    "\n",
    "$$\\mathbf{X}^{+L} = (\\mathbf{V}\\mathbf{S}^\\mathrm{H}\\mathbf{S}\\mathbf{V}^\\mathrm{H})^{-1} (\\mathbf{V}\\mathbf{S}^\\mathrm{H}\\mathbf{U}^\\mathrm{H})$$\n",
    "\n",
    "$$\\mathbf{X}^{+L} = \\mathbf{V} (\\mathbf{S}^\\mathrm{H}\\mathbf{S})^{-1} \\mathbf{V}^\\mathrm{H} \\mathbf{V}\\mathbf{S}^\\mathrm{H}\\mathbf{U}^\\mathrm{H}$$\n",
    "\n",
    "$$\\mathbf{X}^{+L} = \\mathbf{V} (\\mathbf{S}^\\mathrm{H}\\mathbf{S})^{-1} \\mathbf{S}^\\mathrm{H}\\mathbf{U}^\\mathrm{H}$$\n",
    "\n",
    "$$\\mathbf{X}^{+L} = \\mathbf{V} \\mathbf{S}^\\mathrm{+L} \\mathbf{U}^\\mathrm{H}$$\n",
    "\n",
    "allows for a convenient discussion, how singular values act when **mapping column space back to row space**.\n",
    "\n",
    "Considering only one singular value $\\sigma_i$ and its corresponding left/right singular vectors, the left inverse $\\mathbf{S}^\\mathrm{+L} = (\\mathbf{S}^\\mathrm{H}\\mathbf{S})^{-1} \\mathbf{S}^\\mathrm{H}$ reduces to\n",
    "\n",
    "$$\\frac{\\sigma_i}{\\sigma_i^2} = \\frac{1}{\\sigma_i}$$\n",
    "\n",
    "For **very, very small** $\\sigma_i$, the **inversion** thus leads to **huge values**, which might be not meaningful as this (these) weighted $\\mathbf{v}_i$ vector(s) then dominate(s) the row space solution. Small changes in $\\sigma_i$ \n",
    "(i.e. small changes in a badly conditioned feature matrix), then lead to comparably large changes in the row space solution $\\hat{\\beta}$. So-called ridge regression (aka Tikhonov regularization) is a straightforward workaround for ill-conditioned matrices. See stuff below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba47cf9f",
   "metadata": {},
   "source": [
    "### Least Squares Error Solution in QR Domain\n",
    "\n",
    "great material, strong recommendation:\n",
    "- Gilbert Strang (2020): \"Linear Algebra for Everyone\", Wellesley-Cambridge Press, p.170ff\n",
    "- Gilbert Strang (2019): \"Linear Algebra and Learning from Data\", Wellesley-Cambridge Press, p.128ff\n",
    "\n",
    "The normal equation(s)\n",
    "\n",
    "$$\\mathbf{X}^\\mathrm{H}\\mathbf{X} \\hat{\\beta} = \\mathbf{X}^\\mathrm{H}\\mathbf{y}$$\n",
    "\n",
    "can be conveniently given as QR decomposition (recall $\\mathbf{Q}^\\mathrm{H} \\mathbf{Q}=\\mathbf{I}$ due to Gram-Schmidt orthonormalization)\n",
    "\n",
    "$$(\\mathbf{Q R})^\\mathrm{H}\\mathbf{Q R} \\hat{\\beta} = (\\mathbf{Q R})^\\mathrm{H}\\mathbf{y}$$\n",
    "\n",
    "$$\\mathbf{R}^\\mathrm{H} \\mathbf{Q}^\\mathrm{H} \\mathbf{Q R} \\hat{\\mathbf{x}} = (\\mathbf{Q R})^\\mathrm{H}\\mathbf{y}$$\n",
    "\n",
    "$$\\mathbf{R}^\\mathrm{H} \\mathbf{R} \\hat{\\beta} = \\mathbf{R}^\\mathrm{H} \\mathbf{Q}^\\mathrm{H} \\mathbf{y}$$\n",
    "\n",
    "$$\\mathbf{R} \\hat{\\beta} = \\mathbf{Q}^\\mathrm{H} \\mathbf{y}$$\n",
    "\n",
    "We should not expect that algorithms solve\n",
    "\n",
    "$$\\hat{\\mathbf{x}} = \\mathbf{R}^{+L} \\mathbf{Q}^\\mathrm{H} \\mathbf{b}$$\n",
    "\n",
    "with the left inverse $\\mathbf{R}^{+L}$ of upper triangle $\\mathbf{R}$, we should not do this for non-toy-examples as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b9fc068",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(1)\n",
    "mean, stdev = 0, 0.01\n",
    "M = 100\n",
    "N = 3\n",
    "X = rng.normal(mean, stdev, [M, N])\n",
    "print(\"rank =\", matrix_rank(X), \"== number of cols =\", N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e083507b",
   "metadata": {},
   "outputs": [],
   "source": [
    "if matplotlib_widget_flag:\n",
    "    %matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f947dcf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "[Q, R] = qr(X)\n",
    "[U, s, Vh] = svd(X)\n",
    "print(\"sing vals\", s)\n",
    "V = Vh.conj().T\n",
    "\n",
    "# scipy function\n",
    "Xli_pinv = pinv(X)\n",
    "\n",
    "# manual normal equation solver\n",
    "Xli_man = inv(X.conj().T @ X) @ X.conj().T\n",
    "\n",
    "# SVD\n",
    "Si = diagsvd(1 / s, N, M)  # works if array s has only non-zero entries\n",
    "Xli_svd = V @ Si @ U.conj().T\n",
    "\n",
    "# QR\n",
    "Xli_qr = pinv(R) @ Q.conj().T\n",
    "\n",
    "print(\"pinv == inverse via normal eq?\", np.allclose(Xli_pinv, Xli_man))\n",
    "print(\"pinv == inverse via SVD?\", np.allclose(Xli_pinv, Xli_svd))\n",
    "print(\"pinv == inverse via QR?\", np.allclose(Xli_pinv, Xli_qr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb277f66",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create y from one column space entry and one left null space entry\n",
    "# note that we use unit length vectors for convenience: ||e||_2^2 = 1\n",
    "yh = U[:, 0]  # choose one of col space\n",
    "e = U[:, N]  # assuming rank N -> we choose last vector of left null space\n",
    "y = yh + e\n",
    "\n",
    "# find betah in the row space\n",
    "betah = Xli_pinv @ y  # only yh gets mapped back to row space\n",
    "# this is our LS solution betah\n",
    "\n",
    "print(Xli_pinv @ e)  # e is mapped to zero vec\n",
    "print(Xli_pinv @ y)\n",
    "print(Xli_pinv @ yh)\n",
    "print(betah)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4cfafb8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"norm(X @ betah - yh, 2) == 0 -> \", norm(X @ betah - yh, 2))  # == 0\n",
    "\n",
    "print(\"||e||_2^2:\")\n",
    "print(norm(X @ betah - y, 2) ** 2)\n",
    "print(norm(e, 2) ** 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c054fd5b-4bc9-4bea-8068-a1b46b804974",
   "metadata": {},
   "source": [
    "### Ridge Regression / Regularization in SVD Domain\n",
    "\n",
    "The minimization problem\n",
    "\n",
    "$$\\mathrm{min}_\\beta \\left(||\\mathbf{X} \\beta - \\mathbf{y}||_2^2 + \\alpha^2 ||\\beta||_2^2 \\right)$$\n",
    "\n",
    "is known as regression, aka Tikhonov regularization, aka regression with penalty on $||\\beta||_2^2$ using the ridge coefficient coefficient $\\alpha^2$ (this is a hyper parameter in model training, often variable $\\lambda = \\alpha^2$ is used).\n",
    "\n",
    "For limit $\\alpha^2=0$ this minimization problem is identical to above standard least squares solution.\n",
    "\n",
    "The analytical solution is well known as\n",
    "\n",
    "$$\\hat{\\beta}^\\mathrm{Ridge} = \\mathbf{X}^{+\\mathrm{L,Ridge}}  \\mathbf{y}$$\n",
    "\n",
    "with the ridge left inverse matrix given in SVD fashion \n",
    "\n",
    "$$\\mathbf{X}^{+\\mathrm{L,Ridge}} = \\mathbf{V} \\left((\\mathbf{S}^\\mathrm{H}\\mathbf{S} + \\alpha^2 \\mathbf{I})^{-1} \\mathbf{S}^\\mathrm{H}\\right) \\mathbf{U}^\\mathrm{H}$$\n",
    "\n",
    "Considering only one singular value $\\sigma_i$ and its corresponding left/right singular vectors, the left inverse $\\mathbf{S}^\\mathrm{+L,Ridge} = \\left((\\mathbf{S}^\\mathrm{H}\\mathbf{S} + \\alpha^2 \\mathbf{I})^{-1} \\mathbf{S}^\\mathrm{H}\\right)$ reduces to\n",
    "\n",
    "$$\\frac{\\sigma_i}{\\sigma_i^2 + \\alpha^2},$$\n",
    "\n",
    "which can be discussed conveniently with below plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f68be2ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 1 / 10\n",
    "lmb = alpha**2\n",
    "\n",
    "singval = np.logspace(-4, 4, 2**6)\n",
    "# ridge regression\n",
    "inv_singval = singval / (singval**2 + alpha**2)\n",
    "\n",
    "plt.plot(singval, 1 / singval, label=\"no penalty\")\n",
    "plt.plot(singval, inv_singval, label=\"penalty\")\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xticks(10.0 ** np.arange(-4, 5))\n",
    "plt.yticks(10.0 ** np.arange(-4, 5))\n",
    "plt.axis(\"equal\")\n",
    "plt.xlabel(r\"$\\sigma_i$\")\n",
    "plt.ylabel(r\"$\\sigma_i \\,\\,\\,/\\,\\,\\, (\\sigma_i^2 + \\alpha^2)$\")\n",
    "plt.title(r\"ridge penalty $\\alpha =$\" + str(alpha))\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "print(\"alpha =\", alpha, \"alpha^2 = lambda =\", lmb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "205e6667",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(1)\n",
    "mean, stdev = 0, 10\n",
    "M, N = 3, 3\n",
    "A_tmp = rng.normal(mean, stdev, [M, N])\n",
    "[U_tmp, s_tmp, Vh_tmp] = svd(A_tmp)\n",
    "V_tmp = Vh_tmp.conj().T\n",
    "s_tmp = [10, 8, 0.5]  # create sing vals\n",
    "S_tmp = diagsvd(s_tmp, M, N)\n",
    "\n",
    "# create full rank square matrix to work with (no nullspaces except 0-vectors!)\n",
    "A = U_tmp @ S_tmp @ Vh_tmp\n",
    "[U, s, Vh] = svd(A)\n",
    "print(\"A\\n\", A)\n",
    "print(\"rank of A: \", matrix_rank(A))\n",
    "print(\"sigma\", s)\n",
    "S = diagsvd(s, M, N)\n",
    "V = Vh.conj().T\n",
    "\n",
    "# b as column space linear combination\n",
    "b = 1 * U[:, 0] + 1 * U[:, 1] + 1 * U[:, 2]\n",
    "\n",
    "xh = inv(A) @ b\n",
    "print(\"xh =\", xh, \"\\nA xh =\", A @ xh, \"\\nb =\", b)\n",
    "# == (because in b all U weighted with unity gain)\n",
    "print(\"inverted sigma no penalty: \", 1 / s)\n",
    "print(\"||xh||_2^2 =\", norm(xh, 2))\n",
    "print(\"norm of vec: inverted sigma no penalty: \", norm(1 / s, 2))\n",
    "\n",
    "lmb = 2\n",
    "Sli_ridge = inv(S.conj().T @ S + lmb * np.eye(3)) @ S.conj().T\n",
    "Ali_ridge = V @ Sli_ridge @ U.conj().T\n",
    "xh_ridge = Ali_ridge @ b\n",
    "print(\"xh_ridge =\", xh_ridge, \"\\nA xh_ridge =\", A @ xh_ridge, \"\\nb = \", b)\n",
    "# == (because in b all U weighted with unity gain)\n",
    "print(\"inverted sigma with penalty: \", s / (s**2 + lmb))\n",
    "print(\"||xh_ridge||_2^2 =\", norm(xh_ridge, 2))\n",
    "print(\"norm of vec: inverted sigma with penalty: \", norm(s / (s**2 + lmb), 2))\n",
    "\n",
    "fig1 = plt.figure(figsize=(5, 5))\n",
    "ax = plt.axes(projection=\"3d\")\n",
    "w = Vh @ xh\n",
    "wr = Vh @ xh_ridge\n",
    "for n in range(3):\n",
    "    ax.plot(\n",
    "        [0, w[n] * V[0, n]],\n",
    "        [0, w[n] * V[1, n]],\n",
    "        [0, w[n] * V[2, n]],\n",
    "        color=\"C\" + str(n),\n",
    "        lw=1,\n",
    "        ls=\":\",\n",
    "        label=r\"$\\hat{x}$@$v_i$, no penalty\",\n",
    "    )\n",
    "    ax.plot(\n",
    "        [0, wr[n] * V[0, n]],\n",
    "        [0, wr[n] * V[1, n]],\n",
    "        [0, wr[n] * V[2, n]],\n",
    "        color=\"C\" + str(n),\n",
    "        lw=3,\n",
    "        ls=\"-\",\n",
    "        label=r\"$\\hat{x}$@$v_i$, penalty\",\n",
    "    )\n",
    "\n",
    "ax.plot(\n",
    "    [0, xh[0]], [0, xh[1]], [0, xh[2]], \"black\", label=r\"$\\hat{x}$, no penalty\"\n",
    ")\n",
    "ax.plot(\n",
    "    [0, xh_ridge[0]],\n",
    "    [0, xh_ridge[1]],\n",
    "    [0, xh_ridge[2]],\n",
    "    \"C7\",\n",
    "    label=\"$\\hat{x}$, penalty\",\n",
    ")\n",
    "ax.set_xlabel(r\"$x$\")\n",
    "ax.set_ylabel(r\"$y$\")\n",
    "ax.set_zlabel(r\"$z$\")\n",
    "lim = 1\n",
    "ax.set_xlim(-lim, lim)\n",
    "ax.set_ylim(-lim, lim)\n",
    "ax.set_zlim(-lim, lim)\n",
    "ax.set_title(\"V / row space\")\n",
    "plt.legend()\n",
    "\n",
    "\n",
    "fig2 = plt.figure(figsize=(5, 5))\n",
    "ax = plt.axes(projection=\"3d\")\n",
    "w = Vh @ xh\n",
    "wr = Vh @ xh_ridge\n",
    "for n in range(3):\n",
    "    ax.plot(\n",
    "        [0, U[0, n]],\n",
    "        [0, U[1, n]],\n",
    "        [0, U[2, n]],\n",
    "        color=\"C\" + str(n),\n",
    "        lw=2,\n",
    "        ls=\"-\",\n",
    "        label=r\"$u_i$\",\n",
    "    )\n",
    "    ax.plot(\n",
    "        [0, s[n] * U[0, n]],\n",
    "        [0, s[n] * U[1, n]],\n",
    "        [0, s[n] * U[2, n]],\n",
    "        color=\"C\" + str(n),\n",
    "        lw=1,\n",
    "        ls=\":\",\n",
    "        label=r\"$\\sigma_i \\cdot u_i$\",\n",
    "    )\n",
    "\n",
    "ax.plot([0, b[0]], [0, b[1]], [0, b[2]], \"black\", lw=1, label=r\"$b$\")\n",
    "ax.set_xlabel(r\"$x$\")\n",
    "ax.set_ylabel(r\"$y$\")\n",
    "ax.set_zlabel(r\"$z$\")\n",
    "lim = 5\n",
    "ax.set_xlim(-lim, lim)\n",
    "ax.set_ylim(-lim, lim)\n",
    "ax.set_zlim(-lim, lim)\n",
    "ax.set_title(\"U / row space\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05837f9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "if matplotlib_widget_flag:\n",
    "    plt.close(fig1)\n",
    "    plt.close(fig2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b96bff94-6ef3-40e3-8273-2aa218d2edb7",
   "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": 5
}