just_the_code/roch_mmids_chap_svd_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7a24eca6",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "remove-cell",
     "colab-keep"
    ]
   },
   "source": [
    "***\n",
    "\n",
    "*Course:* [Math 535](https://people.math.wisc.edu/~roch/mmids/) - Mathematical Methods in Data Science (MMiDS)  \n",
    "*Chapter:* 3-Singular value decomposition   \n",
    "*Author:* [Sebastien Roch](https://people.math.wisc.edu/~roch/), Department of Mathematics, University of Wisconsin-Madison  \n",
    "*Updated:* March 4, 2024   \n",
    "*Copyright:* &copy; 2024 Sebastien Roch\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ae295bc",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# You will need the files:\n",
    "#     * mmids.py\n",
    "#     * h3n2-snp.csv\n",
    "#     * h3n2-other.csv \n",
    "#     * advertising.csv \n",
    "# from https://github.com/MMiDS-textbook/MMiDS-textbook.github.io/tree/main/utils\n",
    "#\n",
    "# IF RUNNING ON GOOGLE COLAB (RECOMMENDED):\n",
    "# \"Upload to session storage\" from the Files tab on the left\n",
    "# Alternative instructions: https://colab.research.google.com/notebooks/io.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a623de9",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# PYTHON 3\n",
    "import numpy as np\n",
    "from numpy import linalg as LA\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import mmids\n",
    "seed = 535\n",
    "rng = np.random.default_rng(seed)\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd15d93f",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "## Motivating example: visualizing viral genomes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45824cf6",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We consider an application of dimensionality reduction in biology. We will look at SNP data from viruses. A little background first. From [Wikipedia](https://en.wikipedia.org/wiki/Single-nucleotide_polymorphism):\n",
    "\n",
    "> A single-nucleotide polymorphism (SNP) is a substitution of a single nucleotide that occurs at a specific position in the genome, where each variation is present at a level of more than 1% in the population. For example, at a specific base position in the human genome, the C nucleotide may appear in most individuals, but in a minority of individuals, the position is occupied by an A. This means that there is a SNP at this specific position, and the two possible nucleotide variations -- C or A -- are said to be the alleles for this specific position."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ea1f2b4",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Quoting [Jombart et al., BMC Genetics (2010)](https://bmcgenet.biomedcentral.com/articles/10.1186/1471-2156-11-94), we analyze:\n",
    "\n",
    "> the population structure of seasonal influenza A/H3N2 viruses using hemagglutinin (HA) sequences. Changes in the HA gene are largely responsible for immune escape of the virus (antigenic shift), and allow seasonal influenza to persist by mounting yearly epidemics peaking in winter. These genetic changes also force influenza vaccines to be updated on a yearly basis. [...] Assessing the genetic evolution of a pathogen through successive epidemics is of considerable epidemiological interest. In the case of seasonal influenza, we would like to ascertain how genetic changes accumulate among strains from one winter epidemic to the next."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7067e95a",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Some details about the Jombart et al. dataset:\n",
    "\n",
    "> For this purpose, we retrieved all sequences of H3N2 hemagglutinin (HA) collected between 2001 and 2007 available from Genbank. Only sequences for which a location (country) and a date (year and month) were available were retained, which allowed us to classify strains into yearly winter epidemics. Because of the temporal lag between influenza epidemics in the two hemispheres, and given the fact that most available sequences were sampled in the northern hemisphere, we restricted our analysis to strains from the northern hemisphere (latitudes above 23.4°north). The final dataset included 1903 strains characterized by 125 SNPs which resulted in a total of 334 alleles. All strains from 2001 to 2007 were classified into six winter epidemics (2001-2006). This was done by assigning all strains from the second half of the year with those from the first half of the following year. For example, the 2005 winter epidemic comprises all strains collected between the 1st of July 2005 and the 30th of June 2006."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7cdf489",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We load a dataset, which contains a subset of strains from the dataset mentioned above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b83b6d8",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('h3n2-snp.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b3af8db",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "The first five rows are the following."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f02958a",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c2a2e8a",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Overall it contains $1642$ strains. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1206abaf",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "df.shape[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bca88c91",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "The data lives in a $318$-dimensional space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eda16f79",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "df.shape[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a75fef28",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Obviously, vizualizing this data is not straighforward. How can we make sense of it? More specifically, how can we explore any underlying structure it might have. Quoting [Wikipedia](https://en.wikipedia.org/wiki/Exploratory_data_analysis):\n",
    "\n",
    "> In statistics, exploratory data analysis (EDA) is an approach of analyzing data sets to summarize their main characteristics, often using statistical graphics and other data visualization methods. [...] Exploratory data analysis has been promoted by John Tukey since 1970 to encourage statisticians to explore the data, and possibly formulate hypotheses that could lead to new data collection and experiments."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30ec2372",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "In this chapter we will encounter an importatn mathematical technique for dimension reduction, which allow us to explore this data -- and find interesting structure -- in $2$ (rather than $318$!) dimensions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e61c63c",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "## Background: matrix rank; review of eigenvalues and eigenvectors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8565843d",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "**NUMERICAL CORNER:** In Numpy, one can compute the rank of a matrix using the function [`numpy.linalg.matrix_rank`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.matrix_rank.html). We will see later in the course how to compute it using the singular value decomposition (which is how `LA.matrix_rank` does it).\n",
    "\n",
    "Let's try the example above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8446ecbb",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "w1 = np.array([1., 0., 1.])\n",
    "w2 = np.array([0., 1., 1.])\n",
    "w3 = np.array([1., -1., 0.])\n",
    "A = np.stack((w1, w2, w3), axis=-1)\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "033e47b1",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We compute the rank of `A`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b1621c8",
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "LA.matrix_rank(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d86837c",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We take only the first two columns of `A` this time to form `B`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bb4c733",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "B = np.stack((w1, w2),axis=-1)\n",
    "print(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7089ca1",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "LA.matrix_rank(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da7a1de1",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Recall that, in Numpy, `@` is used for matrix product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f39d413f",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "C = np.array([[1., 0., 1.],[0., 1., -1.]])\n",
    "print(C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ad45c68",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "LA.matrix_rank(C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "973126ba",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(B @ C)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6c5c659",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "$\\unlhd$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e802a8ec",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "**NUMERICAL CORNER:** In Numpy, the outer product is computed using [`numpy.outer`](https://numpy.org/doc/stable/reference/generated/numpy.outer.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b679c5e4",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "u = np.array([0., 2., -1.])\n",
    "v = np.array([3., -2.])\n",
    "Z = np.outer(u, v)\n",
    "print(Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8b93100",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(LA.matrix_rank(Z))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7288ecde",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "$\\unlhd$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34705776",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "**NUMERICAL CORNER:** In Numpy, the eigenvalues and eigenvectors of a matrix can be computed using [`numpy.linalg.eig`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.eig.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c4ae28e",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "A = np.array([[2.5, -0.5], [-0.5, 2.5]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17f59e11",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "w, v = LA.eig(A)\n",
    "print(w)\n",
    "print(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82850d26",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "$\\unlhd$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15b90f00",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "**NUMERICAL CORNER:** Hence, we can check whether a matrix is positive semidefinite by computing its eigenvalues using [`numpy.linalg.eig`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.eig.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d64a7056",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "A = np.array([[1, -1], [-1, 1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f25e6205",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "w, v = LA.eig(A)\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71d19844",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "B = np.array([[1, -2], [-2, 1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3d36639",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "z, u = LA.eig(B)\n",
    "print(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3077c710",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "$\\unlhd$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53d5aafa",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "## Power iteration"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f279c91",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "**NUMERICAL CORNER:** We implement the algorithm suggested by the *Power Iteration Lemma*. That is, we compute $B^{k} \\mathbf{x}$, then normalize it. To obtain the corresponding singular value and left singular vector, we use that $\\sigma_1 = \\|A \\mathbf{v}_1\\|$ and $\\mathbf{u}_1 = A \\mathbf{v}_1/\\sigma_1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b8c58e3",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def topsing(A, maxiter=10):\n",
    "    x = rng.normal(0,1,np.shape(A)[1])\n",
    "    B = A.T @ A\n",
    "    for _ in range(maxiter):\n",
    "        x = B @ x\n",
    "    v = x / LA.norm(x)\n",
    "    s = LA.norm(A @ v)\n",
    "    u = A @ v / s\n",
    "    return u, s, v"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d51ccc5c",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We will apply it to our previous two-cluster example. The necessary functions are in [mmids.py](https://raw.githubusercontent.com/MMiDS-textbook/MMiDS-textbook.github.io/main/utils/mmids.py), which is available on the [GitHub of the notes](https://github.com/MMiDS-textbook/MMiDS-textbook.github.io/tree/main). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4530a72d",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "d, n, w = 10, 100, 3.\n",
    "X1, X2 = mmids.two_clusters(d, n, w)\n",
    "X = np.concatenate((X1, X2), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71159161",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111,aspect='equal')\n",
    "ax.scatter(X[:,0], X[:,1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad9dba90",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Let's compute the top singular vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c3f4c70",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "u, s, v = topsing(X)\n",
    "print(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50a0f729",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "This is approximately $-\\mathbf{e}_1$. We get roughly the same answer (possibly up to sign) from Python's [`numpy.linalg.svd`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.svd.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f3039b3",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "u, s, vh = LA.svd(X)\n",
    "print(vh.T[:,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d0a8394",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Recall that, when we applied $k$-means clustering to this example with $d=1000$ dimension, we obtained a very poor clustering. Let's try again after projecting onto the top singular vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d041745e",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "d, n, w = 1000, 100, 3.\n",
    "X1, X2 = mmids.two_clusters(d, n, w)\n",
    "X = np.concatenate((X1, X2), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7819215",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111,aspect='equal')\n",
    "ax.scatter(X[:,0], X[:,1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce54d693",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "assign = mmids.kmeans(X, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af1cb596",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111,aspect='equal')\n",
    "ax.scatter(X[:,0], X[:,1], c=assign)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "def009ee",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Let's try again, but after projecting on the top singular vector. Recall that this corresponds to finding the best one-dimensional approximating subspace. The projection can be computed using the truncated SVD $Z= U_{(1)} \\Sigma_{(1)} V_{(1)}^T$. We can interpret the rows of $U_{(1)} \\Sigma_{(1)}$ as the coefficients of each data point in the basis $\\mathbf{v}_1$. We will work in that basis. We need one small hack: because our implementation of $k$-means clustering expects data points in at least $2$ dimension, we add a column of $0$'s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bd8c83a",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "u, s, v = topsing(X)\n",
    "Xproj = np.stack((u*s, np.zeros(np.shape(X)[0])), axis=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa342428",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111,aspect='equal')\n",
    "ax.scatter(Xproj[:,0], Xproj[:,1])\n",
    "plt.ylim([-3,3])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3d7e5a5",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "There is a small -- yet noticeable -- gap around 0."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fd4f905",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "A histogram of the first component of `Xproj` gives a better sense of the density of points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23542848",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "plt.hist(Xproj[:,0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74aba306",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We run $k$-means clustering on the projected data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8f59e55",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "assign = mmids.kmeans(Xproj, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e8439fb",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111,aspect='equal')\n",
    "ax.scatter(X[:,0], X[:,1], c=assign)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2819ea34",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Much better. We will give an explanation of this outcome in an upcoming (optional) subsection. In essence, quoting [BHK, Section 7.5.1]:\n",
    "\n",
    "> [...] let's understand the central advantage of doing the projection to [the top $k$ right singular vectors]. It is simply that for any reasonable (unknown) clustering of data points, the projection brings data points closer to their cluster centers."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "233e420d",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Finally, looking at the top right singular vector (or its first ten entries for lack of space), we see that it does align quite well (but not perfectly) with the first dimension. In the next (optional) section, we try again with the top two singular vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17fd61ce",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(v[:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ba89244",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "$\\unlhd$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edd5d5d1",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "**NUMERICAL CORNER:** We implement this last algorithm. We will need our previous implementation of *Gram-Schimdt*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8231f3dc",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def svd(A, l, maxiter=100):\n",
    "    V = rng.normal(0,1,(np.size(A,1),l))\n",
    "    for _ in range(maxiter):\n",
    "        W = A @ V\n",
    "        Z = A.T @ W\n",
    "        V, R = mmids.gramschmidt(Z)\n",
    "    W = A @ V\n",
    "    S = [LA.norm(W[:, i]) for i in range(np.size(W,1))]\n",
    "    U = np.stack([W[:,i]/S[i] for i in range(np.size(W,1))],axis=-1)\n",
    "    return U, S, V"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4a4baa4",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Note that above we avoided forming the matrix $A^T A$. With a small number of iterations, that approach potentially requires fewer arithmetic operations overall and it allows to take advantage of the possible sparsity of $A$ (i.e. the fact that it may have many zeros)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2191dc47",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We apply it again to our two-cluster example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "530b8f27",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "d, n, w = 1000, 100, 3.\n",
    "X1, X2 = mmids.two_clusters(d, n, w)\n",
    "X = np.concatenate((X1, X2), axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dc409a2",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Let's try again, but after projecting on the top two singular vectors. Recall that this corresponds to finding the best two-dimensional approximating subspace. The projection can be computed using the truncated SVD $Z= U_{(2)} \\Sigma_{(2)} V_{(2)}^T$. We can interpret the rows of $U_{(2)} \\Sigma_{(2)}$ as the coefficients of each data point in the basis $\\mathbf{v}_1,\\mathbf{v}_2$. We will work in that basis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dd0cb69",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "U, S, V = svd(X, 2)\n",
    "Xproj = np.stack((U[:,0]*S[0], U[:,1]*S[1]), axis=-1)\n",
    "assign = mmids.kmeans(Xproj, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20dba73c",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111,aspect='equal')\n",
    "ax.scatter(X[:,0], X[:,1], c=assign)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cb68508",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Finally, looking at the first two right singular vectors, we see that the first one does align quite well with the first dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ced7bfce",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(np.stack((V[:,0], V[:,1]), axis=-1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "641029b9",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "$\\unlhd$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bb42ceb",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "## Application to dimensionality reduction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7c684e8",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "**NUMERICAL CORNER:** Having established a formal connection between PCA and SVD, we implement PCA using our SVD algorithm (in [mmids.py](https://raw.githubusercontent.com/MMiDS-textbook/MMiDS-textbook.github.io/main/utils/mmids.py)). We perform mean centering (now is the time to read that quote about the importance of mean centering again), but not the optional standardization. We use the fact that, in Numpy, subtracting a matrix by a vector whose dimension matches the number of columns performs row-wise subtraction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b09d0c2",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def pca(X, l, maxiter=100):\n",
    "    mean = np.mean(X, axis=0)  # Compute mean of each column\n",
    "    Y = X - mean # Mean center each column\n",
    "    U, S, V = mmids.svd(Y, l, maxiter=maxiter)\n",
    "    return U[:,:l] @ np.diag(S[:l])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e75d5b8",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We apply it to the Gaussian Mixture Model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "278c27ed",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "d, n, w = 1000, 100, 3.\n",
    "X = mmids.two_mixed_clusters(d, n, w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c20ec7c",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "T = pca(X, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "415b757b",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Plotting the result, we see that PCA does succeed in finding the main direction of variation, with the first principal component aligning well with the first coordinate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57b32a54",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111,aspect='equal')\n",
    "ax.scatter(T[:,0], T[:,1], s=10)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "275541b2",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Note however that the first two principal components (especially the second one) in fact \"capture more noise\" than what can be seen in the first two coordinates, a form of overfitting. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2233dcd",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111,aspect='equal')\n",
    "ax.scatter(X[:,0], X[:,1], s=10)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66fe24c6",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "**TRY IT!** Compute the first two right singular vectors $\\mathbf{v}_1$ and $\\mathbf{v}_2$ of $X$ after mean centering. Do they align well with the first and second standard basis vectors $\\mathbf{e}_1$ and $\\mathbf{e}_2$? Why or why not? ([Open in Colab](https://colab.research.google.com/github/MMiDS-textbook/MMiDS-textbook.github.io/blob/main/just_the_code/roch_mmids_chap_svd_notebook.ipynb))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "380b49d3",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "$\\unlhd$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b8df0e8",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We return to our motivating example. We apply PCA to our genetic dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "986c897b",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We load the dataset again and examine its first rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97d4dbaa",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('h3n2-snp.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cea95e3",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d239c820",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Recall that it contains $1642$ strains and lives in a $318$-dimensional space. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "783d4f92",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "df.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bdfb609",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Our goal is to find a \"good\" low-dimensional representation of the data. Ten dimensions will do here. We use PCA. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ebb77ef",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "A = df[[df.columns[i] for i in range(1,len(df.columns))]].to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfc8b8d3",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "n_dims = 10\n",
    "T = pca(A, n_dims)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a0afd96",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We plot the first two principal components, and see what appears to be some potential structure. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fb29e88",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111,aspect='equal')\n",
    "ax.scatter(T[:,0], T[:,1], s=10)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "464048db",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "There seems to be some reasonably well-defined clusters in this projection. We use $k$-means to identiy clusters. We take advantage of the implementation in scikit-learn, [sklearn.cluster.KMeans](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html). By default, it finds $8$ clusters. The clusters can be extracted from the attribute `labels_`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4637703",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c76d532c",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "n_clusters = 8\n",
    "kmeans = KMeans(n_clusters=n_clusters, init='k-means++', n_init=10).fit(T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6f66de3",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "assign = kmeans.labels_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92ecdcb8",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "To further reveal the structure, we look at our the clusters spread out over the years. That information is in a separate file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcc45855",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "dfoth = pd.read_csv('h3n2-other.csv')\n",
    "dfoth.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04f30947",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "year = dfoth['year'].to_numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67806035",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "For each cluster, we plot how many of its data points come from a specific year. Each cluster has a different color."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90dcc23d",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "for i in range(n_clusters):\n",
    "    unique, counts = np.unique(year[assign == i], return_counts=True)\n",
    "    ax.bar(unique, counts, label=i)\n",
    "\n",
    "ax.set(xlim=(2001, 2007), xticks=np.arange(2002, 2007))\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ceb2157d",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Remarkably, we see that each cluster comes mostly from one year or two consecutive ones. In other words, the clustering in this low-dimensional projection captures some true underlying structure that is not explicitly in the genetic data on which it is computed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d922ea0",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Going back to the first two principal components, we color the points on the scatterplot by year. (We use [`legend_elements()`](https://matplotlib.org/stable/api/collections_api.html#matplotlib.collections.PathCollection.legend_elements) for automatic legend creation.) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95ccfc8c",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, aspect='equal')\n",
    "scatter = ax.scatter(T[:,0], T[:,1], s=10,  c=year, label=year)\n",
    "plt.legend(*scatter.legend_elements())\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4afb347a",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "To some extent, one can \"see\" the virus evolving from year to year. The $x$-axis in particular seems to correlate strongly with the year, in the sense that samples from later years tend to be towards one side of the plot. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73b39208",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "To further quantify this observation, we use [numpy.corrcoef](https://numpy.org/doc/stable/reference/generated/numpy.corrcoef.html) to compute the correlation coefficients between the year and the first $10$ principal components. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a5094f6",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "corr = np.zeros(n_dims)\n",
    "for i in range(n_dims):\n",
    "    corr[i] = np.corrcoef(np.stack((T[:,i], year)))[0,1]\n",
    "\n",
    "print(corr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea4e97db",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Indeed, we see that the first three or four principal components correlate well with the year."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d3d97ed",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "## Further applications of the SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1005447",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "**NUMERICAL CORNER:** In Numpy, the Frobenius norm of a matrix can be computed using the default of the function `numpy.linalg.norm` while the induced norm can be computed using the same function with [`ord` parameter set to `2`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.norm.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63cb0cc2",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "A = np.array([[1., 0.],[0., 1.],[0., 0.]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8d46a15",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "LA.norm(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1a5984f",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "LA.norm(A, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f99994b",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "$\\unlhd$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7919f520",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "**NUMERICAL CORNER:** In Numpy, the pseudoinverse of a matrix can be computed using the function [`numpy.linalg.pinv`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.pinv.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0376ed6d",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "M = np.array([[1.5, 1.3], [1.2, 1.9], [2.1, 0.8]])\n",
    "print(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54103fb2",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "Mp = LA.pinv(M)\n",
    "print(Mp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1dbd5d6d",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "Mp @ M"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49bb1726",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Let's try our previous example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9437b01a",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "A = np.array([[1., 0.], [-1., 0.]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e8fb307",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "Ap = LA.pinv(A)\n",
    "print(Ap)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cff02f51",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "$\\unlhd$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bb3d45d",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "**NUMERICAL CORNER:** In Numpy, the condition number of a matrix can be computed using the function [`numpy.linalg.cond`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.cond.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7827bc5f",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "For example, orthogonal matrices have condition number $1$, the lowest possible value for it (Why?). That indicates that orthogonal matrices have good numerical properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63ecf9ca",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "q = 1/np.sqrt(2)\n",
    "Q = np.array([[q, q], [q, -q]])\n",
    "print(Q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93866e57",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "LA.cond(Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbc80e6e",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "In contrast, matrices with nearly linearly dependent columns have large condition numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec379d77",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "eps = 1e-6\n",
    "A = np.array([[q, q], [q, q+eps]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e43a3d44",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "LA.cond(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c264aa7",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "u, s, vh = LA.svd(A)\n",
    "print(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "799fc641",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We compute the solution to $A \\mathbf{x} = \\mathbf{b}$ when $\\mathbf{b}$ is the left singular vector of $A$ corresponding to the largest singular value. Recall that in the proof of the *Conditioning of Matrix-Vector Multiplication Theorem*, we showed that the worst case bound is achieved when $\\mathbf{z} = \\mathbf{b}$ is right singular vector of $M= A^{-1}$ corresponding to the lowest singular value. In a previous example, given a matrix $A = \\sum_{j=1}^n \\sigma_j \\mathbf{u}_j \\mathbf{v}_j^T$ in compact SVD form, we derived a compact SVD for the inverse as\n",
    "\n",
    "$$\n",
    "A^{-1} = \\sigma_n^{-1} \\mathbf{v}_n \\mathbf{u}_n^T + \\sigma_{n-1}^{-1} \\mathbf{v}_{n-1} \\mathbf{u}_{n-1}^T + \\cdots + \\sigma_1^{-1} \\mathbf{v}_1 \\mathbf{u}_1^T.\n",
    "$$\n",
    "\n",
    "Here, compared to the SVD of $A$, the order of the singular values is reversed and the roles of the left and right singular vectors are exchanged. So we take $\\mathbf{b}$ to be the top left singular vector of $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76c174d1",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "b = u[:,0]\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b4a3d9b",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x = LA.solve(A,b)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1914b6c",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We make a small perturbation in the direction of the second right singular vector. Recall that in the proof of the *Conditioning of Matrix-Vector Multiplication Theorem*, we showed that the worst case is achieved when $\\mathbf{d} = \\delta\\mathbf{b}$ is top right singular vector of $M = A^{-1}$. By the argument above, that is the left singular vector of $A$ corresponding to the lowest singular value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f39d3bd",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "delta = 1e-6\n",
    "bp = b + delta*u[:,1]\n",
    "print(bp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4257fea",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "The relative change in solution is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65881e7e",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "xp = LA.solve(A,bp)\n",
    "print(xp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "459f6fff",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "(LA.norm(x-xp)/LA.norm(x))/(LA.norm(b-bp)/LA.norm(b))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31ea0275",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Note that this is exactly the condition number of $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d8789f3",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "$\\unlhd$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb4c04f3",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "**NUMERICAL CORNER:** We give a quick example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d0441cf",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "A = np.array([[1., 101.],[1., 102.],[1., 103.],[1., 104.],[1., 105]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d929a03",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "LA.cond(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ded71737",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "LA.cond(A.T @ A)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbca2bb1",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "This observation -- and the resulting increased numerical instability -- is one of the reasons we previously developed an alternative approach to the least-squares problem. Quoting [Sol, Section 5.1]:\n",
    "\n",
    "> Intuitively, a primary reason that $\\mathrm{cond}(A^T A)$ can be large is that columns of $A$ might\n",
    "look “similar” [...] If two columns $\\mathbf{a}_i$ and $\\mathbf{a}_j$ satisfy $\\mathbf{a}_i \\approx \\mathbf{a}_j$, then the least-squares residual length $\\|\\mathbf{b} - A \\mathbf{x}\\|_2$ will not suffer much if we replace multiples of $\\mathbf{a}_i$ with multiples of $\\mathbf{a}_j$ or vice versa. This wide range of nearly—but not completely—equivalent solutions yields poor conditioning. [...] To solve such poorly conditioned problems, we will employ an alternative technique with closer attention to the column space of $A$ rather than employing row operations as in Gaussian elimination. This strategy identifies and deals with such near-dependencies explicitly, bringing about greater numerical stability."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dd89f3c",
   "metadata": {
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "$\\unlhd$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31ea96c6",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "**NUMERICAL CORNER:** Here is a numerical example taken from [[TB](https://books.google.com/books/about/Numerical_Linear_Algebra.html?id=JaPtxOytY7kC), Lecture 19]. We will approximate the following function with a polynomial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "785f6d55",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "n = 100 \n",
    "t = np.arange(n)/(n-1)\n",
    "b = np.exp(np.sin(4 * t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4859aba8",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "plt.plot(t, b)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "faf3899b",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We use a [Vandermonde matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix), which can be constructed using [`numpy.vander`](https://numpy.org/doc/stable/reference/generated/numpy.vander.html), to perform polynomial regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3afd4580",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "m = 17\n",
    "A = np.vander(t, m, increasing=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a76dca64",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "The condition numbers of $A$ and $A^T A$ are both high in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79031573",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print(LA.cond(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcc8b214",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "print(LA.cond(A.T @ A))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33d3760f",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We first use the normal equations and plot the residual vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3ad28f4",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "xNE = LA.solve(A.T @ A, A.T @ b)\n",
    "print(LA.norm(b - A@xNE))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66091345",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "plt.plot(t, b - A@xNE)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e34377d8",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We then use `numpy.linalg.qr` to compute the QR solution instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb05e9e0",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "Q, R = LA.qr(A)\n",
    "xQR = mmids.backsubs(R, Q.T @ b)\n",
    "print(LA.norm(b - A@xQR)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ca88395",
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "plt.plot(t, b - A@xNE)\n",
    "plt.plot(t, b - A@xQR)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b752d1d",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "$\\unlhd$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "038b8bd5",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "$\\newcommand{\\horz}{\\rule[.5ex]{2.5ex}{0.5pt}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6847f5a",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "## More advanced material: low-rank approximation; why project; additional proofs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a3b4e63",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "**NUMERICAL CORNER:** We return to our example with the two Gaussian clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6cc6179",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "d, n, w = 1000, 100, 3.\n",
    "X1, X2 = mmids.two_clusters(d, n, w)\n",
    "X = np.concatenate((X1, X2), axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4548288",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "In reality, we cannot compute the matrix norms of $X-C$ and $X_k-C$ as the true centers are not known. But, because this is simulated data, we happen to know the truth and we can check the validity of our results in this case. The centers are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3e8ebd5",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "C1 = np.stack([np.concatenate(([-w], np.zeros(d-1))) for _ in range(n)])\n",
    "C2 = np.stack([np.concatenate(([w], np.zeros(d-1))) for _ in range(n)])\n",
    "C = np.concatenate((C1, C2), axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3c28a57",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "We use [`numpy.linalg.svd`](https://numpy.org/doc/stable/reference/generated/numpy.linalg.svd.html) function to compute the norms from the formulas in the *Matrix Norms and Singular Values Lemma*. First, we observe that the singular values of $X-C$ are decaying slowly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b932ce4",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "uc, sc, vhc = LA.svd(X-C)\n",
    "plt.plot(sc)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15326f6f",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "The $k$-means objective with respect to the true centers under the full-dimensional data is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13875f5d",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "frob = np.sum(sc**2)\n",
    "print(frob)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cea7ad8a",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "while the square of the top singular value of $X-C$ is only: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "153bc8b3",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "top_sval_sq = sc[0]**2\n",
    "print(top_sval_sq)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0afca03",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "Finally, we compute the $k$-means objective with respect to the true centers under the projected one-dimensional data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2852846b",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "u, s, vh = LA.svd(X)\n",
    "frob_proj = np.sum((s[0] * np.outer(u[:,0],vh[:,0]) - C)**2)\n",
    "print(frob_proj)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2655a1ca",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": [
     "colab-keep"
    ]
   },
   "source": [
    "$\\unlhd$"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}