06a_latent_models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Music Machine Learning - Latent models and clustering (Full Optional)\n",
    "\n",
    "### Author: Philippe Esling (esling@ircam.fr)\n",
    "\n",
    "In this course we will cover\n",
    "1. Unsupervised learning and [clustering](#clustering).\n",
    "2. Motivating the need for [latent variable](#latent) models.\n",
    "2. Approaching the problem naïvely with the [k-Means](#kmeans) algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unsupervised learning and clustering\n",
    "\n",
    "Up to now, we have dealt with *supervised* problems (classification, regression), which implies that we have *paired* data. This implies that we both have data examples $\\{\\mathbf{x} \\in \\mathbb{R}^{n}\\}$ and their corresponding labels $\\{\\mathbf{y}_{i} \\in \\mathbb{R}^{n}\\}$, and we are trying to find the function to link these two sets, or equivalently, the joint distribution $p(\\mathbf{y} \\vert \\mathbf{x})$.\n",
    "\n",
    "Now imagine that we only have access to $\\{\\mathbf{x} \\in \\mathbb{R}^{n}\\}$ (a set of data without particular annotations). For instance, consider the following set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helper_plot import hdr_plot_style\n",
    "from sklearn.datasets.samples_generator import make_blobs\n",
    "import matplotlib.pyplot as plt\n",
    "hdr_plot_style();\n",
    "X, y_true = make_blobs(n_samples=400, centers=4, cluster_std=0.60, random_state=0)\n",
    "X = X[:, ::-1]; plt.figure(figsize=(10, 8));\n",
    "plt.scatter(X[:, 0], X[:, 1], c='b', s=60, edgecolor='w'); plt.grid(True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we can clearly see that there are some *clusters* underlying our data. This is a form of *hidden structure*, which we might want to model. This type of learning is called *unsupervised learning* (here specifically, clustering). We might want to model the complete distribution $\\mathbf{x} \\sim p(\\mathbf{x})$. However, this might be very hard to define directly, as we know that there is some *hidden information* (data that we did not observe, such as the cluster identity here)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latent'> </a>\n",
    "## Latent variable model\n",
    "\n",
    "A latent variable is simply a variable that you never observe. In our previous example, these variables would be the *membership* of different points to clusters. These latent variables can be introduced to\n",
    "- Quantify uncertainty\n",
    "- Handle missing values\n",
    "\n",
    "In our case, we introduce a latent variable $\\bz \\in \\mathbb{R}^{n}$, which corresponds to the cluster identity (or membership of each point to a given cluster), where $z_{i} = C_{k}$, and $C_k \\in \\mathbb{N}$ is a cluster identifier.\n",
    "\n",
    "The goal of *unsupervised learning* and *latent variable models* is to model the joint distribution $p(\\mathbf{x}, \\mathbf{z})$, allowing to obtain $p(\\mathbf{z}\\vert\\mathbf{x})$ (infering the hidden structure) and even $p(\\mathbf{x} \\vert \\mathbf{z})$ (generating new examples from a given structure). Here, we will try to find the hidden *membership* of points to different clusters. We will start with a naïve and deterministic approach to this problem.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true
   },
   "source": [
    "## Simple clustering : k-Means\n",
    "\n",
    "We recall here the basic behind the k-Means algorithm. Given a set of observations ($\\mathbf{x}_{1}$, $\\mathbf{x}_{2}$, ..., $\\mathbf{x}_{n}$), where each observation is a $d$-dimensional real vector, k-means clustering aims to partition the $n$ observations into $k \\leq n$ sets $S = \\{S_{1}, S_{2}, \\cdots, S_{k}\\}$ so as to minimize the within-cluster sum of squares (WCSS) (sum of distance functions of each point in the cluster to the $K$ center). In other words, its objective is to find:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "{\\underset {\\mathbf {S} }{\\operatorname {arg\\,min} }}\\sum _{i=1}^{k}\\sum _{\\mathbf {x} \\in S_{i}}\\left\\|\\mathbf {x} -{\\boldsymbol {\\mu }}_{i}\\right\\|^{2} \n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Given an initial set of $k$ means $\\{m_{1}^{(0)},\\cdots,m_{k}^{(0)}\\}$, the algorithm proceeds by alternating between two steps:\n",
    "\n",
    "**Assignment step**: Assign each observation to the cluster whose mean yields the least within-cluster sum of squares (WCSS). Since the sum of squares is the squared Euclidean distance, this is intuitively the \"nearest\" mean. (Mathematically, this means partitioning the observations according to the Voronoi diagram generated by the means).\n",
    "\n",
    "$$\n",
    "S_{i}^{(t)}={\\big \\{}x_{p}:{\\big \\|}x_{p}-m_{i}^{(t)}{\\big \\|}^{2}\\leq {\\big \\|}x_{p}-m_{j}^{(t)}{\\big \\|}^{2}\\ \\forall j,1\\leq j\\leq k{\\big \\}}\n",
    "$$\n",
    "\n",
    "where each $x_{p}$ is assigned to exactly one $S^{(t)}$, even if it could be assigned to two or more of them.  \n",
    "\n",
    "**Update step**: Calculate the new means to be the centroids of the observations in the new clusters.\n",
    "\n",
    "$$\n",
    "m_{i}^{(t+1)}={\\frac {1}{|S_{i}^{(t)}|}}\\sum _{x_{j}\\in S_{i}^{(t)}}x_{j}\n",
    "$$\n",
    "\n",
    "Since the arithmetic mean is a least-squares estimator, this also minimizes the within-cluster sum of squares (WCSS) objective.\n",
    "The algorithm has converged when the assignments no longer change. Since both steps optimize the WCSS objective, and there only exists a finite number of such partitionings, the algorithm must converge to a (local) optimum. There is no guarantee that the global optimum is found using this algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Observing k-Means with `scikit-learn`\n",
    "\n",
    "As usual, we can first witness the effect of the kMeans algorithm using scikit. In the following code, we instantiate a `KMeans` object, which will allow to cluster our dataset of points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Plot the data with K Means Labels\n",
    "from sklearn.cluster import KMeans\n",
    "kmeans = KMeans(4, random_state=2, max_iter=10, n_init=1, init='random')\n",
    "labels = kmeans.fit(X).predict(X)\n",
    "plt.figure(figsize=(10,8)); plt.grid(True);\n",
    "plt.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='magma', edgecolor='w');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true
   },
   "source": [
    "From an intuitive standpoint, we might expect that the clustering assignment for some points is more certain than others: for example, there appears to be a very slight overlap between the two middle clusters, such that we might not have complete confidence in the cluster assigment of points between them.\n",
    "Unfortunately, the *k*-means model has no intrinsic measure of probability or uncertainty of cluster assignments (although it may be possible to use a bootstrap approach to estimate this uncertainty).\n",
    "For this, we must think about generalizing the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementing the k-Means algorithm\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "**Exercise**  \n",
    "\n",
    "  1. Update the ''kmeans'' function to implement the clustering algorithm.\n",
    "  2. Perform the plot using the spread function to display the results of clustering.\n",
    "  3. Compare the results depending on the number of clusters (example are displayed below).\n",
    "  4. What observations can you make on the quality of these clusters ?\n",
    "  5. Compare your results with the `scikit` kMeans function.\n",
    "  \n",
    "***\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def kmeans(X, k):\n",
    "    # X: d x n data matrix\n",
    "    # k: number of seeds\n",
    "    \n",
    "    ######################\n",
    "    # YOUR CODE GOES HERE\n",
    "    ######################\n",
    "\n",
    "    return label\n",
    "\n",
    "labels = kmeans(X.transpose(), 4)\n",
    "plt.figure(figsize=(10,8)); plt.grid(True);\n",
    "plt.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='magma', edgecolor='w');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Limitations of k-Means\n",
    "One way to think about the *k*-means model is that it places a circle (or, in higher dimensions, a hyper-sphere) at the center of each cluster, with a radius defined by the most distant point in the cluster.\n",
    "This radius acts as a hard cutoff for cluster assignment within the training set: any point outside this circle is not considered a member of the cluster.\n",
    "We can visualize this cluster model with the following function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true
   },
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "from scipy.spatial.distance import cdist\n",
    "\n",
    "def plot_kmeans(kmeans, X, n_clusters=4, rseed=0, ax=None):\n",
    "    labels = kmeans.fit_predict(X)\n",
    "    plt.figure(figsize=(10,8))\n",
    "    # plot the input data\n",
    "    ax = ax or plt.gca()\n",
    "    ax.axis('equal')\n",
    "    ax.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='magma', zorder=2, edgecolor='w')\n",
    "\n",
    "    # plot the representation of the KMeans model\n",
    "    centers = kmeans.cluster_centers_\n",
    "    radii = [cdist(X[labels == i], [center]).max()\n",
    "             for i, center in enumerate(centers)]\n",
    "    for c, r in zip(centers, radii):\n",
    "        ax.add_patch(plt.Circle(c, r, fc='#AAAAAA', lw=3, alpha=0.5, zorder=1))\n",
    "\n",
    "kmeans = KMeans(n_clusters=4, random_state=0)\n",
    "plot_kmeans(kmeans, X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true
   },
   "source": [
    "An important observation for *k*-means is that these cluster models *must be circular*: *k*-means has no built-in way of accounting for oblong or elliptical clusters.\n",
    "So, for example, if we take the same data and transform it, the cluster assignments end up becoming muddled:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "editable": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "rng = np.random.RandomState(13)\n",
    "X_stretched = np.dot(X, rng.randn(2, 2))\n",
    "kmeans = KMeans(n_clusters=4, random_state=0)\n",
    "plot_kmeans(kmeans, X_stretched)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true
   },
   "source": [
    "By eye, we recognize that these transformed clusters are non-circular, and thus circular clusters would be a poor fit.\n",
    "Nevertheless, *k*-means is not flexible enough to account for this, and tries to force-fit the data into four circular clusters.\n",
    "This results in a mixing of cluster assignments where the resulting circles overlap: see especially the bottom-right of this plot.\n",
    "\n",
    "These two disadvantages of *k*-means—its lack of flexibility in cluster shape and lack of probabilistic cluster assignment—mean that for many datasets (especially low-dimensional datasets) it may not perform as well as you might hope."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probabilistic clustering\n",
    "\n",
    "We might imagine addressing the weaknesses of the *k*-means model. For example, you could measure uncertainty in cluster assignment by comparing the distances of each point to *all* cluster centers, rather than focusing on just the closest. You might also imagine allowing the cluster boundaries to be ellipses rather than circles, so as to account for non-circular clusters.\n",
    "It turns out these are two essential components of a *probabilistic clustering*\n",
    "\n",
    "- Usual clustering (kMeans) is done with a _hard_ decision threshold\n",
    "    * Each point belongs to one cluster\n",
    "    * For instance the basic k-Means algorithm\n",
    "$$\n",
    "C_{id} = f(\\mathbf{x})\n",
    "$$\n",
    "- We might want to do probabilistic (soft) clustering\n",
    "    * Each point belongs more or less to all clusters\n",
    "    * Defines a probability of belonging\n",
    "$$\n",
    "p(C_{id} \\vert x)\n",
    "$$\n",
    "\n",
    "Using probabilistic models has several interests, as it allows to provide finer evaluation of the clusters and their parameters. However, this will require to estimate both the parameters and cluster identity altogether. This is the goal of the **Expectation-Maximization** algorithm, that we will see in the next notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus - Hierarchical clustering for audio thumbnailing\n",
    "\n",
    "we will perform a simple shot at the problem of *audio structure discovery* and *audio thumbnailing*. The idea is to try to automatically infer the structure of a piece of music from its inner similarities in an *unsupervised way*. As previously, we will rely on our helper functions to import some audio dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from helper_data import import_dataset, compute_transforms, compute_features\n",
    "class_path = 'data/music_speech'\n",
    "# 0.1 - Import the classification dataset\n",
    "data_struct = import_dataset(class_path, 'music-speech')\n",
    "# 0.2 - Pre-process the audio to obtain spectral transforms \n",
    "data_struct = compute_transforms(data_struct)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We briefly recall here that the principle of *hierarchical agglomerative clustering* is to start with a singleton cluster, and clusters are iteratively merged until one single cluster remains. This results in a \"cluster tree,\" which is also called dendrogram. The opposite approach (starting with one cluster and divide into clusters until only singleton clusters remain) is called *divisive hierarchical clustering*. The algorithm can be summarized via the following pseudocode\n",
    "\n",
    "**1**: Compute a distance or similarity matrix.  \n",
    "**2**: Each data point is represented as a singleton cluster.  \n",
    "**3**: Repeat  \n",
    "**4**: $\\;\\;$ Merge two closest clusters (e.g., based on distance between most similar or dissimilar members).  \n",
    "**5**: $\\;\\;$ Update the distance (or similarity) matrix.  \n",
    "**6**: Until one single cluster remains.  \n",
    "\n",
    "As the algorithm is very easy to implement, we will learn how to apply it on more complex problems. The idea here is to use a *smoothed version* (time-wise) of audio tracks and try to find the *structure* of this music in an unsupervised way. Therefore, we will try to find similarities. To do so, rely on the documentation for the `cluster` and `linkage` function to find a way to perform hierarchical clustering on the set of spectrogram windows.\n",
    "\n",
    "***\n",
    "\n",
    "**Exercise**  \n",
    "<div markdown=\"1\">  \n",
    "\n",
    "  1. Update the loop to perform hierarchical clustering.\n",
    "  2. Evaluate different distances and tracks to see the effects.\n",
    "  3. Implement your **own distance function** between spectrograms.\n",
    "  4. Find a way to plug your distance matrix inside the clustering.\n",
    "  \n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.cluster.hierarchy import dendrogram, linkage\n",
    "import matplotlib.gridspec as gridspec\n",
    "song_ex = 11\n",
    "nbClusters = 20\n",
    "# Extract the Constant-Q\n",
    "curCQT = data_struct[\"spectrum_CQT\"][song_ex];\n",
    "nb_points = curCQT.shape[1]\n",
    "# Perform a smoothed version\n",
    "smooth_target = 200\n",
    "smooth_wins = np.floor(nb_points / smooth_target) * 2\n",
    "smoothCQT = np.zeros((smooth_target, curCQT.shape[0]))\n",
    "# Prepare set of windows\n",
    "first_win = smooth_wins / 2\n",
    "last_win = nb_points - (smooth_wins / 2)\n",
    "win_set = np.round(np.linspace(first_win, last_win, smooth_target))\n",
    "win_starts = (win_set - first_win)\n",
    "win_ends = win_starts + (smooth_wins)\n",
    "# Go through the points\n",
    "for t in range(smooth_target):\n",
    "    winCQT = curCQT[:, int(win_starts[t]):int(win_ends[t])]\n",
    "    smoothCQT[t, :] = np.mean(winCQT, axis=1)\n",
    "\n",
    "######################\n",
    "# YOUR CODE GOES HERE\n",
    "######################\n",
    "    \n",
    "# Plot a stylish dendrogram\n",
    "fig = plt.figure(figsize=(12,8))\n",
    "gs1 = gridspec.GridSpec(3, 1, height_ratios=[3, 1, 2])\n",
    "gs1.update(wspace=0.025, hspace=0.01)\n",
    "plt.subplot(gs1[0])\n",
    "Rdict = dendrogram(Z, leaf_rotation=90., leaf_font_size=8., p=nbClusters, count_sort=False, distance_sort=False, color_threshold=0.8)\n",
    "ax = plt.gca()\n",
    "ax.axis('off')\n",
    "plt.subplot(gs1[1])\n",
    "ax = plt.gca()\n",
    "ax.axis('off')\n",
    "for i in range(len(Rdict[\"ivl\"])):\n",
    "    plt.plot([i, Rdict[\"ivl\"][i]], [1, 0])\n",
    "plt.autoscale(enable=True, axis='x', tight=True)\n",
    "plt.subplot(gs1[2])\n",
    "plt.imshow(np.flipud(smoothCQT.transpose()), aspect='auto')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "em_assignment.ipynb",
   "provenance": [],
   "version": "0.3.2"
  },
  "kernelspec": {
   "display_name": "Python 3.10.4 ('atiam_ml')",
   "language": "python",
   "name": "python3"
  },
  "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.4"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "vscode": {
   "interpreter": {
    "hash": "19d703e906ada558fae81aed9eb14243028cc1b26072879ea22a3f4f6919ed00"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}