exercise04_rightinv.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 Right Matrix Inverse\n",
    "\n",
    "## Objectives\n",
    "\n",
    "- Matrix A of dimension (M x N)\n",
    "- SVD for a **flat/fat** matrix, we assume a matrix with **full row rank** $r=M$\n",
    "- Four subspaces in SVD domain\n",
    "- Projection matrices\n",
    "- Right 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 Flat/Fat, Full Row Rank Matrix A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 3  # number of rows\n",
    "N = 7  # number of cols\n",
    "\n",
    "rank = min(M, N)  # set desired rank == full row 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(\"flat/fat matrix with full row rank\")\n",
    "print(\"-> matrix U contains only the column space\")\n",
    "print(\"-> left 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 right of diagonal part\n",
    "print(\"V =\\n\", V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Four Subspaces in SVD Domain"
   ]
  },
  {
   "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:])  # for full row rank this is only the zero vector\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:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Right 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 bottom zero block\")\n",
    "print(\"Si =\\n\", Si)\n",
    "# right inverse using 'inverse' SVD:\n",
    "Ari = V @ Si @ U.conj().T\n",
    "# right inverse using a dedicated pinv algorithm\n",
    "# proper choice is done by pinv() itself\n",
    "Ari_pinv = pinv(A)\n",
    "print(\"pinv() == right inverse via SVD?\", np.allclose(Ari, Ari_pinv))\n",
    "print(\"S @ Si = \\n\", S @ Si, \"\\nyields MxM identity matrix\")\n",
    "print(\"A @ Ari = \\n\", A @ Ari, \"\\nyields MxM identity matrix\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Projection Matrices for the Right Inverse Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v_row = 1 * V[:, 0, None]\n",
    "v_null = 2 * V[:, rank, None]\n",
    "v_tmp = v_row + v_null\n",
    "\n",
    "# projection onto row space\n",
    "P_CAH = Ari @ A\n",
    "print(\n",
    "    \"projection matrix P_CAH projects V-space stuff to row space:\\n\",\n",
    "    \"P_CAH @ v_tmp == v_row:\",\n",
    "    np.allclose(P_CAH @ v_tmp, v_row),\n",
    ")\n",
    "\n",
    "# projection onto column space\n",
    "# full rank and identity since we don't have a left null space\n",
    "# so each column space vector is projected onto itself\n",
    "P_CA = A @ Ari  # = I_MxM\n",
    "print(\n",
    "    \"projection matrix P_CA projects a column space vector onto itself:\\n\",\n",
    "    \"P_CA @ U[:, 0] == U[:, 0]:\",\n",
    "    np.allclose(P_CA @ U[:, 0], U[:, 0]),\n",
    ")\n",
    "\n",
    "# projection onto null space, for flat/fat this space might be very large\n",
    "P_NA = np.eye(N, N) - P_CAH\n",
    "print(\n",
    "    \"projection matrix P_NA projects V-space stuff to null space:\\n\",\n",
    "    \"P_NA @ v_tmp == v_null:\",\n",
    "    np.allclose(P_NA @ v_tmp, v_null),\n",
    ")\n",
    "\n",
    "# projection onto left null space\n",
    "P_NAH = np.eye(M, M) - P_CA  # == null matrix\n",
    "print(\"projection matrix P_NAH is the MxM zero matrix\\n\", P_NAH)"
   ]
  },
  {
   "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
}