{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython import display"
   ]
  },
  {
   "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: Singular Value Decomposition (SVD)\n",
    "\n",
    "## Objectives\n",
    "\n",
    "- SVD\n",
    "- Four subspaces of a matrix in SVD domain\n",
    "- Rank-1 matrix superposition\n",
    "\n",
    "## Special Python Packages\n",
    "\n",
    "Some convenient functions are found in `scipy.linalg`, some in `numpy.linalg` \n",
    "\n",
    "## Highly Recommended Material\n",
    "- Kenji Hiranabe's [Graphic notes on Gilbert Strang's \"Linear Algebra for Everyone\"](https://github.com/kenjihiranabe/The-Art-of-Linear-Algebra/blob/main/The-Art-of-Linear-Algebra.pdf)\n",
    "\n",
    "- Gilbert Strang's\n",
    "    - [A 2020 Vision of Linear Algebra, Zoom Notes, 2021](https://ocw.mit.edu/courses/res-18-010-a-2020-vision-of-linear-algebra-spring-2020/resources/zoomnotes_18-010/) and/or\n",
    "    - [A 2020 Vision of Linear Algebra, Slides, 2020](https://ocw.mit.edu/courses/res-18-010-a-2020-vision-of-linear-algebra-spring-2020/resources/mitres_18_010s20_la_slides/)\n",
    "    - [Lecture Notes for Linear Algebra Sample Sections](https://math.mit.edu/~gs/LectureNotes/samples.pdf) and/or\n",
    "    - [brilliant textbooks on linear algebra](https://math.mit.edu/~gs/) and/or\n",
    "    - [video lecture in course 18.06](https://www.youtube.com/watch?v=TX_vooSnhm8) and/or\n",
    "    - [video lecture in course 18.065](https://www.youtube.com/watch?v=rYz83XPxiZo)\n",
    "- Steven L. Brunton's [video lecture(s) on SVD](https://www.youtube.com/watch?v=nbBvuuNVfco)"
   ]
  },
  {
   "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=3, floatmode=\"fixed\", suppress=True)\n",
    "\n",
    "rng = np.random.default_rng(1234)\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": [
    "## Generate a Matrix with desired rank and dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 7  # number of rows\n",
    "N = 4  # number of cols\n",
    "# set desired rank, here N-1 for a rank deficient matrix,\n",
    "# i.e. neither full column nor full row rank\n",
    "rank = min(M, N) - 1\n",
    "\n",
    "if False:  # another small toy example\n",
    "    M = 4  # number of rows\n",
    "    N = 3  # number of cols\n",
    "    rank = 2\n",
    "\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.integers(-3, 3, M) + 1j * rng.integers(-3, 3, M)\n",
    "        row = rng.integers(-3, 3, N) + 1j * rng.integers(-3, 3, N)\n",
    "        A += np.outer(col, row)  # superposition of rank-1 matrices\n",
    "else:\n",
    "    dtype = \"int\"\n",
    "    A = np.zeros([M, N], dtype=dtype)\n",
    "    for i in range(rank):\n",
    "        col = rng.integers(-3, 3, M)\n",
    "        row = rng.integers(-3, 3, N)\n",
    "        A += np.outer(col, row)  # superposition of rank-1 matrices\n",
    "# check if rng produced desired matrix rank\n",
    "print(\"        rank of A:\", matrix_rank(A))\n",
    "print(\"A =\\n\", A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could otherwise consider the matrix **A** on slide 4/30 from the brilliant [A 2020 Vision of Linear Algebra, Slides, 2020](https://ocw.mit.edu/courses/res-18-010-a-2020-vision-of-linear-algebra-spring-2020/resources/mitres_18_010s20_la_slides/). This is the case we have discussed in detail in the 2nd tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if True:\n",
    "    A = np.array([[1, 4, 5], [3, 2, 5], [2, 1, 3]])\n",
    "    M, N = A.shape[0], A.shape[1]\n",
    "    rank = matrix_rank(A)\n",
    "    print(\"rank r =\", rank)\n",
    "    print(\"A=\\n\", A)\n",
    "    [U, s, Vh] = svd(A)\n",
    "    V = Vh.T.conj()\n",
    "\n",
    "    print(\n",
    "        \"column space is spanned by plane in 3D space with these two vectors (1st two columns):\"\n",
    "    )\n",
    "    print(A[:, :rank])\n",
    "    print(\n",
    "        \"only two independent columns (the rank tells us this already), because\"\n",
    "    )\n",
    "    print(\n",
    "        \"A[:,0,None] + A[:,1,None] == A[:,2,None]:\\n\",\n",
    "        A[:, 0, None] + A[:, 1, None] == A[:, 2, None],\n",
    "    )\n",
    "\n",
    "    print(\"left null space is spanned by line in 3D space with this vector:\")\n",
    "    lns = np.array([[-1], [7], [-10]])\n",
    "    print(lns)\n",
    "    print(\n",
    "        \"check that we get the zero vector for left null space stuff:\\n\",\n",
    "        A.T @ lns,\n",
    "    )\n",
    "\n",
    "    print(\"null space is spanned by line in 3D space with this vector:\")\n",
    "    ns = np.array([[-1], [-1], [1]])\n",
    "    print(ns)\n",
    "    print(\"check that we get the zero vector for null space stuff:\\n\", A @ ns)\n",
    "\n",
    "    print(\"row space is spanned by plane in 3D space with two vectors\")\n",
    "    rs = np.array([[1, 0], [0, 1], [1, 1]])\n",
    "    print(\"rs=\\n\", rs)\n",
    "    print(\"project rs into V space to get weights for the V vectors\")\n",
    "    print(Vh @ rs)  # third weights are zero, as these belong to nullspace\n",
    "    print(\"linear combinations of 1st and 2nd V vectors yield 2 rs vectors:\")\n",
    "    print(V @ (Vh @ rs)[:, 0, None])\n",
    "    print(V @ (Vh @ rs)[:, 1, None])\n",
    "\n",
    "    print(\"we could check that the null space, the left null space,\")\n",
    "    print(\"the column space and the row space that we came up here\")\n",
    "    print(\"is also (and actually in a nicer way, because\")\n",
    "    print(\n",
    "        \"SVD orthonormality holds) spanned via the V and U matrices of the SVD\"\n",
    "    )\n",
    "    print(\"see code below\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SVD of Matrix A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[U, s, Vh] = svd(A)\n",
    "S = diagsvd(s, M, N)  # for full SVD the matrix S has same dim as A\n",
    "V = Vh.conj().T\n",
    "\n",
    "print(\"U =\\n\", U)\n",
    "# number of non-zero sing values along diag must match rank\n",
    "print(\"\\nnon-zero singular values: \", s[:rank], \"\\n\")\n",
    "print(\"S =\\n\", S)\n",
    "print(\"V =\\n\", V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "\n",
    "# image from the brilliant project Kenji Hiranabe: Graphic notes on\n",
    "# Gilbert Strang's \"Linear Algebra for Everyone\"\n",
    "# found at https://github.com/kenjihiranabe/The-Art-of-Linear-Algebra\n",
    "# CC0-1.0 License\n",
    "Image(\n",
    "    \"https://github.com/kenjihiranabe/The-Art-of-Linear-Algebra/raw/2357e987993ba88eb34bbe16e038ce1b150c4878/SVD.png\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Four Subspaces in the SVD Domain\n",
    "\n",
    "We denote \n",
    "- column space $C(\\mathbf{A})$, aka *image*, *range*\n",
    "- left null space $N(\\mathbf{A}^\\mathrm{H})$\n",
    "- row space $C(\\mathbf{A}^\\mathrm{H})$\n",
    "- null space $N(\\mathbf{A})$, aka *kernel*\n",
    "\n",
    "The vectors of column space $C(\\mathbf{A})$ and left null space $N(\\mathbf{A}^\\mathrm{H})$ are related to the matrix $\\mathbf{U}$ (matrix output).\n",
    "\n",
    "The vectors of row space $C(\\mathbf{A}^\\mathrm{H})$ and null space $N(\\mathbf{A})$ are related to the matrix $\\mathbf{V}$ (matrix input).\n",
    "\n",
    "Since $\\mathbf{U}$ and $\\mathbf{V}$ are unitary matrices, they fulfill **orthonormality** and therefore the properties of orthogonal subspaces\n",
    "\n",
    "$$C(\\mathbf{A}) \\perp N(\\mathbf{A}^\\mathrm{H})$$\n",
    "\n",
    "$$C(\\mathbf{A}^\\mathrm{H}) \\perp N(\\mathbf{A})$$\n",
    "\n",
    "immediately can be deduced. In other words, since we know by the property of SVD, that the matrices $\\mathbf{U}$ and $\\mathbf{V}$ are orthonormal, we see that\n",
    "\n",
    "- column space $C(\\mathbf{A})$ is orthogonal to left null space $N(\\mathbf{A}^\\mathrm{H})$, since both spaces are spanned by the dedicated vectors in $\\mathbf{U}$\n",
    "- row space $C(\\mathbf{A}^\\mathrm{H})$ is orthogonal to null space $N(\\mathbf{A})$, since both spaces are spanned by the dedicated vectors in $\\mathbf{V}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# image from the brilliant project Kenji Hiranabe: Graphic notes on\n",
    "# Gilbert Strang's \"Linear Algebra for Everyone\"\n",
    "# found at https://github.com/kenjihiranabe/The-Art-of-Linear-Algebra\n",
    "# CC0-1.0 License\n",
    "Image(\n",
    "    \"https://github.com/kenjihiranabe/The-Art-of-Linear-Algebra/raw/2357e987993ba88eb34bbe16e038ce1b150c4878/4-Subspaces.png\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "also see the nice graphics of the 4 subspaces with indicated SVD vectors in this wiki http://mlwiki.org/index.php/Four_Fundamental_Subspaces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all stuff that is in matrix U (matrix output)\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 the 0 vector is the left nullspace\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 (matrix input)\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 the 0 vector is the null space\n",
    "print(\"null space (orthogonal to row space):\")\n",
    "print(V[:, rank:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we use null_space() function\n",
    "\n",
    "# this very often yields the same vectors as doing the SVD manually:\n",
    "print(\"null space: \\n\", null_space(A))\n",
    "\n",
    "# this might be different from U[:,rank:], but spans the same space:\n",
    "print(\"left null space: \\n\", null_space(A.conj().T))\n",
    "# we might want to check this (how should we do this actually?!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pure Row Space --> Column Space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# x as linear combination\n",
    "# of the first two right singular vectors\n",
    "# i.e. from the row space:\n",
    "x = 1 * V[:, 0, None] + 2 * V[:, 1, None]\n",
    "# let matrix A act on x\n",
    "print(A @ x)\n",
    "# result must be identical with the linear combination\n",
    "# of the first two left singular vectors\n",
    "# i.e. from the column space, weighted by their dedicated singular values and\n",
    "# same weights 1,2 as above\n",
    "print(1 * S[0, 0] * U[:, 0, None] + 2 * S[1, 1] * U[:, 1, None])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Row Space + Null Space --> Column Space + Null Vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# x vector as a linear combination of row space and null space\n",
    "x_r = 3 * V[:, 0, None]  # from row space\n",
    "x_n = 4 * V[:, rank, None]  # from null space\n",
    "print(\"x from pure null space must yield zero vector in the left nullspace:\")\n",
    "print(\"A @ x_n =\\n\", A @ x_n)\n",
    "x = x_r + x_n  # linear combination\n",
    "print(\"let matrix A act on x:\")\n",
    "print(A @ x)\n",
    "print(\n",
    "    \"note: we transpose the following identical results to have convenient row display\"\n",
    ")\n",
    "print((A @ x).T)\n",
    "print((A @ x_r).T)\n",
    "print((3 * S[0, 0] * U[:, 0, None]).T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Superposition of Rank-1 Matrices\n",
    "\n",
    "Eckhart-Young theorem, cf. https://en.wikipedia.org/wiki/Low-rank_approximation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ar = np.zeros([M, N], dtype=\"complex128\")\n",
    "for r in range(rank):\n",
    "    print(\"\\nSum of\", r + 1, \"rank-1 matrices -> check (A-Ar)\")\n",
    "    Ar += S[r, r] * np.outer(U[:, r], V[:, r].conj())\n",
    "    print(\n",
    "        \"Frobenius norm (i.e. root(sum squared sing val)): \",\n",
    "        norm(A - Ar, \"fro\"),\n",
    "    )\n",
    "    print(\"Spectral norm / L2 norm (i.e. sigma1)\", norm(A - Ar, 2))\n",
    "    print(\"allclose(A, Ar):\", np.allclose(A, Ar))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the last case `Ar` is fully reconstructed from rank-1 matrices, yielding Frobenius and spectral norm for `A-Ar` of 0 (numerical precision does not give us zero). We can also check with allclose()"
   ]
  },
  {
   "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
}