exercise04_leftinv.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",
    "# 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",
   "metadata": {},
   "source": [
    "# Exercise 4: SVD and Left Matrix Inverse\n",
    "\n",
    "## Objectives\n",
    "\n",
    "- Matrix A of dimension (M x N)\n",
    "- SVD for a **tall/thin** matrix, we assume a matrix with **full column rank** $r=N$\n",
    "- Four subspaces in SVD domain\n",
    "- Projection matrices\n",
    "- Least squares / normal equation(s)\n",
    "- Left inverse\n",
    "\n",
    "## Special Python Packages\n",
    "\n",
    "Some convenient functions are found in `scipy.linalg`, some in `numpy.linalg` "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some Initial Python Stuff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.linalg import svd, diagsvd, inv, pinv, null_space, norm\n",
    "from numpy.linalg import matrix_rank\n",
    "\n",
    "np.set_printoptions(precision=2, floatmode=\"fixed\", suppress=True)\n",
    "\n",
    "rng = np.random.default_rng(1234)\n",
    "mean, stdev = 0, 1\n",
    "\n",
    "# we might convince ourselves that all works for complex data as well\n",
    "# then the ^H operator (conj().T) needs to be used instead of just .T\n",
    "use_complex = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SVD of Tall/Thin, Full Column Rank Matrix A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 7  # number of rows\n",
    "N = 3  # number of cols\n",
    "\n",
    "rank = min(M, N)  # set desired rank == full column rank == independent columns\n",
    "print(\"desired rank of A:\", rank)\n",
    "\n",
    "if use_complex:\n",
    "    dtype = \"complex128\"\n",
    "    A = np.zeros([M, N], dtype=dtype)\n",
    "    for i in range(rank):\n",
    "        col = rng.normal(mean, stdev, M) + 1j * rng.normal(mean, stdev, M)\n",
    "        row = rng.normal(mean, stdev, N) + 1j * rng.normal(mean, stdev, N)\n",
    "        A += np.outer(col, row)  # superposition of rank-1 matrices\n",
    "else:\n",
    "    dtype = \"float64\"\n",
    "    A = np.zeros([M, N], dtype=dtype)\n",
    "    for i in range(rank):\n",
    "        col = rng.normal(mean, stdev, M)\n",
    "        row = rng.normal(mean, stdev, N)\n",
    "        A += np.outer(col, row)  # superposition of rank-1 matrices\n",
    "# check if rng produced desired rank\n",
    "print(\"        rank of A:\", matrix_rank(A))\n",
    "print(\"tall/thin matrix with full column rank\")\n",
    "print(\"-> matrix V contains only the row space\")\n",
    "print(\"-> null space is only the zero vector\")\n",
    "print(\"A =\\n\", A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[U, s, Vh] = svd(A, full_matrices=True)\n",
    "S = diagsvd(s, M, N)  # for full SVD the matrix S has same dim as A\n",
    "V = Vh.conj().T\n",
    "Uh = U.conj().T\n",
    "\n",
    "print(\"U =\\n\", U)\n",
    "# number of non-zero sing values along diag must match rank\n",
    "print(\"non-zero singular values: \", s[:rank])\n",
    "print(\"S =\\n\", S)  # contains 0 Matrix below diagonal part\n",
    "print(\"V =\\n\", V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Four Subspaces in SVD Domain\n",
    "\n",
    "The null space $N(\\mathbf{A})$ of the tall/thin, full column rank matrix $\\mathbf{A}$ is only $\\mathbf{0}$, i.e. $N(\\mathbf{A})=\\mathbf{0}$. Except for $\\mathbf{x}=\\mathbf{0}$, all other $\\mathbf{x}$ are mapped to the column space $C(\\mathbf{A})$. This, however, requires, that the $\\mathbf{V}$ matrix completely spans the row space and no $\\mathbf{v}$ vectors span a dedicated null space.\n",
    "\n",
    "The tall/thin, full column rank matrix $\\mathbf{A}$ spans a rather large left null space $N(\\mathbf{A}^\\mathrm{H})$ with dimension $M-\\mathrm{rank}(A)$.\n",
    "\n",
    "We, therefore, here deal with a linear set of equations with **more equations than unknowns** ($M>N$, more rows than columns) , i.e. the **over-determined** case. For this case, we can find a solution in the **least-squares error** sense by help of the **left inverse** as discussed below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all stuff that is in matrix U\n",
    "print(\"U =\\n\", U)\n",
    "\n",
    "# column space C(A)\n",
    "print(\"\\ncolumn space (orthogonal to left null space):\")\n",
    "print(U[:, :rank])\n",
    "\n",
    "# left null space, if empty only 0 vector\n",
    "print(\"left null space (orthogonal to column space):\")\n",
    "print(U[:, rank:])\n",
    "\n",
    "print(\"###\")\n",
    "\n",
    "# all stuff that is in matrix V\n",
    "print(\"\\nV =\\n\", V)\n",
    "\n",
    "# row space\n",
    "print(\"\\nrow space (orthogonal to null space):\")\n",
    "print(V[:, :rank])\n",
    "\n",
    "# null space N(A), if empty only 0 vector\n",
    "print(\"null space (orthogonal to row space):\")\n",
    "print(V[:, rank:])  # for full column rank this is only the zero vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Left Inverse via SVD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[U, s, Vh] = svd(A, full_matrices=True)\n",
    "V = Vh.conj().T\n",
    "Uh = U.conj().T\n",
    "\n",
    "Si = diagsvd(1 / s, N, M)  # works if array s has only non-zero entries\n",
    "print(\"Inverse singular value matrix with right zero block\")\n",
    "print(\"Si =\\n\", Si)\n",
    "# left inverse using 'inverse' SVD:\n",
    "Ali = V @ Si @ Uh\n",
    "# left inverse using a dedicated pinv algorithm\n",
    "# proper choice is done by pinv() itself\n",
    "Ali_pinv = pinv(A)\n",
    "print(\"pinv == left inverse via SVD?\", np.allclose(Ali, Ali_pinv))\n",
    "print(\"Si @ S\\n\", Si @ S, \"\\nyields NxN identity matrix\")\n",
    "print(\"Ali @ A = \\n\", Ali @ A, \"\\nyields NxN identity matrix\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Projection Matrices for the Left Inverse Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_col = U[:, 0, None]\n",
    "u_left_null = U[:, rank, None]\n",
    "u_tmp = u_col + u_left_null\n",
    "\n",
    "# projection onto row space == I_NxN\n",
    "# full rank and identity since we don't have a null space\n",
    "# so each vector of the row space is projected onto itself\n",
    "P_CAH = Ali @ A\n",
    "print(\n",
    "    \"projection matrix P_CAH projects a row space vector onto itself:\\n\",\n",
    "    \"P_CAH @ V[:,0] == V[:,0]:\",\n",
    "    np.allclose(P_CAH @ V[:, 0], V[:, 0]),\n",
    ")\n",
    "\n",
    "# projection onto column space\n",
    "P_CA = A @ Ali\n",
    "print(\n",
    "    \"projection matrix P_CA projects U-space stuff to column space:\\n\",\n",
    "    \"P_CA @ (u_tmp) == u_col:\",\n",
    "    np.allclose(P_CA @ (u_tmp), u_col),\n",
    ")\n",
    "\n",
    "# projection onto null space == null matrix\n",
    "P_NA = np.eye(N, N) - P_CAH\n",
    "print(\"projection matrix P_NA is a null matrix\\nP_NA=\\n\", P_NA)\n",
    "\n",
    "# projection onto left null space\n",
    "P_NAH = np.eye(M, M) - P_CA\n",
    "print(\n",
    "    \"projection matrix P_NAH projects U-space stuff to left null space:\\n\",\n",
    "    \"P_NAH @ (u_tmp)==u_lnull\",\n",
    "    np.allclose(P_NAH @ (u_tmp), u_left_null),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# design a vector:\n",
    "# one entry from column space + one entry from left null space\n",
    "# so we take some special left singular vectors:\n",
    "b = U[:, 0, None] + U[:, rank, None]  # same as u_tmp above\n",
    "print(\"b==u_tmp\", np.allclose(b, u_tmp))\n",
    "\n",
    "# the vector b is a linear combination and lives in column+left null spaces\n",
    "# with above introduced projection matrices we can project b\n",
    "# (i) onto column space\n",
    "bhat = P_CA @ b  # project b to column space C(A)\n",
    "print(\"bhat==U[:, 0, None]:\", np.allclose(bhat, U[:, 0, None]))\n",
    "# and\n",
    "# (ii) onto left null space\n",
    "e = P_NAH @ b  # project b to left null space N(A^H)\n",
    "print(\"e==U[:, 0, None]:\", np.allclose(e, U[:, rank, None]))\n",
    "\n",
    "# to find x_hat, i.e. the LS solution of the inverse problem\n",
    "# we bring b to row space via left inverse:\n",
    "# only the column space part of b (i.e. bhat) is brought to the row space\n",
    "# we can never map back a 'zero' (i.e. e)\n",
    "print(\"x_hat = Ali @ b == Ali @ bhat:\", np.allclose(Ali @ b, Ali @ bhat))\n",
    "# we expect that this is the scaled first right singular vector\n",
    "print(\"x_hat=\\n\", V[:, 0, None] / S[0, 0])\n",
    "print(Ali @ bhat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Left Inverse from Normal Equation(s) / Least Squares Error Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Derivation 1\n",
    "\n",
    "Vector addition from example above with an error term $\\mathbf{e}$\n",
    "\n",
    "$\\hat{\\mathbf{b}} + \\mathbf{e} = \\mathbf{b} \\rightarrow \\mathbf{e} = \\mathbf{b} - \\hat{\\mathbf{b}}$\n",
    "\n",
    "We know, that pure row space $\\hat{\\mathbf{x}}$ maps to pure column space $\\hat{\\mathbf{b}}$\n",
    "\n",
    "$\\hat{\\mathbf{b}} = \\mathbf{A} \\hat{\\mathbf{x}}$\n",
    "\n",
    "Inserting this\n",
    "\n",
    "$\\mathbf{e} = \\mathbf{b} - \\hat{\\mathbf{b}} = \\mathbf{b} - \\mathbf{A} \\hat{\\mathbf{x}}$\n",
    "\n",
    "The vector $\\mathbf{A} \\mathbf{x}$ is always living in the column space no matter what $\\mathbf{x}$ is constructed of.\n",
    "\n",
    "The vector $\\mathbf{e}$ is orthogonal to column space (since it lives in left null space).\n",
    "\n",
    "So, we know that the inner product must solve to zero\n",
    "\n",
    "$(\\mathbf{A} \\mathbf{x})^\\mathrm{H} \\mathbf{e} = 0 \\rightarrow \\mathbf{x}^\\mathrm{H} \\mathbf{A}^\\mathrm{H} (\\mathbf{b} - \\mathbf{A} \\hat{\\mathbf{x}}) = 0$\n",
    "\n",
    "Rearranging yields\n",
    "\n",
    "$\\mathbf{x}^\\mathrm{H} \\mathbf{A}^\\mathrm{H} \\mathbf{b} = \\mathbf{x}^\\mathrm{H} \\mathbf{A}^\\mathrm{H} \\mathbf{A} \\hat{\\mathbf{x}}$\n",
    "\n",
    "and by canceling $\\mathbf{x}^\\mathrm{H} $ the famous normal equation is obtained\n",
    "\n",
    "$\\mathbf{A}^\\mathrm{H} \\mathbf{b} = \\mathbf{A}^\\mathrm{H} \\mathbf{A} \\hat{\\mathbf{x}}$\n",
    "\n",
    "This can be solved using the inverse of $\\mathbf{A}^\\mathrm{H} \\mathbf{A}$ (this matrix is full rank and therefore invertible)\n",
    "\n",
    "$(\\mathbf{A}^\\mathrm{H} \\mathbf{A})^{-1} \\mathbf{A}^\\mathrm{H} \\mathbf{b} = (\\mathbf{A}^\\mathrm{H} \\mathbf{A})^{-1} (\\mathbf{A}^\\mathrm{H} \\mathbf{A}) \\hat{\\mathbf{x}}$\n",
    "\n",
    "Since $(\\mathbf{A}^\\mathrm{H} \\mathbf{A})^{-1} (\\mathbf{A}^\\mathrm{H} \\mathbf{A}) = \\mathbf{I}$ holds, we get the solution in least-squares error sense for $\\mathbf{x}$ in the row space of $\\mathbf{A}$\n",
    "\n",
    "$(\\mathbf{A}^\\mathrm{H} \\mathbf{A})^{-1} \\mathbf{A}^\\mathrm{H} \\mathbf{b} = \\hat{\\mathbf{x}}$\n",
    "\n",
    "We find the **left inverse** of $\\mathbf{A}$ as\n",
    "\n",
    "$\\mathbf{A}^{+L} = (\\mathbf{A}^\\mathrm{H} \\mathbf{A})^{-1} \\mathbf{A}^\\mathrm{H}$\n",
    "\n",
    "such that left multiplying\n",
    "\n",
    "$\\mathbf{A}^{+L} \\mathbf{A} = \\mathbf{I}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Derivation 2\n",
    "\n",
    "Here, the origin of the naming *least squares error* is made more obvious.\n",
    "\n",
    "We set up an optimization problem defining **least** amount of **squared** length of the **error** vector\n",
    "\n",
    "$\\mathrm{min}_\\mathbf{x} ||\\mathbf{e}||^2_2 = \\mathrm{min}_\\mathbf{x} ||\\mathbf{b} - \\mathbf{A} {\\mathbf{x}}||_2^2$\n",
    "\n",
    "We could solve this with help of calculus. But we have a nice tool, i.e. the properties of subspaces, that not requires pages of calculation:\n",
    "\n",
    "We must find the minimum distance from $\\mathbf{b}$ to the column space of $\\mathbf{A}$.\n",
    "\n",
    "This can be only achieved if the error vector $\\mathbf{e}=\\mathbf{b} - \\mathbf{A} {\\mathbf{x}}$ is orthogonal to the column space of $\\mathbf{A}$. \n",
    "\n",
    "This in turn means that $\\mathbf{e}$ must live in the left null space of $\\mathbf{A}$, i.e. $\\mathbf{e} \\in N(\\mathbf{A}^\\mathrm{H})$. \n",
    "\n",
    "By definition of left nullspace we have $\\mathbf{A}^\\mathrm{H} \\mathbf{e} = \\mathbf{0}$. \n",
    "\n",
    "So, let us insert $\\mathbf{e}$ into $\\mathbf{A}^\\mathrm{H} \\mathbf{e} = \\mathbf{0}$:\n",
    "\n",
    "$\\mathbf{A}^\\mathrm{H} (\\mathbf{b} - \\mathbf{A} {\\mathbf{x}}) = \\mathbf{0}$\n",
    "\n",
    "$\\mathbf{A}^\\mathrm{H} \\mathbf{b} = \\mathbf{A}^\\mathrm{H} \\mathbf{A} {\\mathbf{x}}$\n",
    "\n",
    "The optimum $\\mathbf{x}$ which solves the problem is just as above in derivation 1\n",
    "\n",
    "$(\\mathbf{A}^\\mathrm{H} \\mathbf{A})^{-1} \\mathbf{A}^\\mathrm{H} \\mathbf{b} = \\hat{\\mathbf{x}}$\n",
    "\n",
    "And again, we find the **left inverse** of $\\mathbf{A}$ as\n",
    "\n",
    "$\\mathbf{A}^{+L} = (\\mathbf{A}^\\mathrm{H} \\mathbf{A})^{-1} \\mathbf{A}^\\mathrm{H}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ali_normaleq = inv(A.conj().T @ A) @ A.conj().T\n",
    "\n",
    "xhat = Ali_normaleq @ b  # LS solution in row space\n",
    "print(\"xhat = \", xhat)\n",
    "bhat = A @ xhat  # map to column space\n",
    "\n",
    "# thus this is the projection matrix that maps b (column + left nullspace stuff)\n",
    "# to the column space of A in the LS sense\n",
    "P_CA2 = A @ Ali_normaleq\n",
    "\n",
    "print(\"P_CA == P_CA2 ?\", np.allclose(P_CA, P_CA2))  # check with above result\n",
    "print(\n",
    "    \"P_CA2 @ b == bhat:\", np.allclose(P_CA2 @ b, bhat)\n",
    ")  # check that both outputs are equal"
   ]
  },
  {
   "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
}