{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# TOPIC 0 \n", "\n", "# Introduction\n", "\n", "## 2 Clustering: an objective and an algorithm\n", "\n", "***\n", "*Course:* [Math 535](http://www.math.wisc.edu/~roch/mmids/) - Mathematical Methods in Data Science (MMiDS) \n", "*Author:* [Sebastien Roch](http://www.math.wisc.edu/~roch/), Department of Mathematics, University of Wisconsin-Madison \n", "*Updated:* Sep 8, 2020 \n", "*Copyright:* © 2020 Sebastien Roch\n", "***" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Consider the following fundamental problem in data science. We are given $n$ vectors $\\mathbf{x}_1,\\ldots,\\mathbf{x}_n$ in $\\mathbb{R}^d$. Our goal is to find a good clustering: loosely speaking, we want to partition these data points into $k$ disjoint subsets -- or clusters -- with small pairwise distances within clusters and large pairwise distances across clusters. To make this rather imprecise problem more precise, we consider a specific objective function known as the $k$-means objective." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "![clustering](https://upload.wikimedia.org/wikipedia/commons/thumb/c/c8/Cluster-2.svg/320px-Cluster-2.svg.png)\n", "\n", "[(Source)](https://commons.wikimedia.org/wiki/File:Cluster-2.svg)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Fix a number of clusters $k$. Formally, we think of a clustering as a partition.\n", "\n", "**Definition (Partition):** A partition of $[n] = \\{1,\\ldots,n\\}$ of size $k$ is a collection of non-empty subsets $C_1,\\ldots,C_k \\subseteq [n]$ that:\n", "- are pairwise disjoint, i.e., $C_i \\cap C_j = \\emptyset$, $\\forall i \\neq j$\n", "- cover all of $[n]$, i.e., $\\cup_{i=1}^k C_i = [n]$. $\\lhd$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### 2.1 The $k$-means objective" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Given $n$ vectors $\\mathbf{x}_1,\\ldots,\\mathbf{x}_n$ in $\\mathbb{R}^d$ and a partition $C_1, \\ldots, C_k \\subseteq [n]$, it will be useful to have notation for the corresponding cluster assignment: we define $c(i) = j$ if $i \\in C_j$. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Under the $k$-means objective, the cost of $C_1,\\ldots,C_k$ is then defined as\n", "\n", "$$\n", "\\mathcal{G}(C_1,\\ldots,C_k) = \\min_{\\boldsymbol{\\mu}_1,\\ldots,\\boldsymbol{\\mu}_k \\in \\mathbb{R}^d} \n", "G(C_1,\\ldots,C_k; \\boldsymbol{\\mu}_1, \\ldots, \\boldsymbol{\\mu}_k)\n", "$$\n", "\n", "where\n", "\n", "\n", "\\begin{align}\n", "G(C_1,\\ldots,C_k; \\boldsymbol{\\mu}_1, \\ldots, \\boldsymbol{\\mu}_k)\n", "&= \\sum_{j=1}^k \\sum_{i \\in C_j} \\|\\mathbf{x}_i - \\boldsymbol{\\mu}_j\\|^2\n", "= \\sum_{i=1}^n \\|\\mathbf{x}_i - \\boldsymbol{\\mu}_{c(i)}\\|^2.\n", "\\end{align}\n", "\n", "\n", "Here $\\boldsymbol{\\mu}_i \\in \\mathbb{R}^d$ is the representative -- or center -- of cluster $C_i$. Note that $\\boldsymbol{\\mu}_i$ need not be one of the $\\mathbf{x}_j$'s." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Our goal is to find a partition that minimizes $\\mathcal{G}(C_1,\\ldots,C_k)$, i.e., to solve the problem:\n", "\n", "$$\n", "\\min_{C_1,\\ldots,C_k} \\mathcal{G}(C_1,\\ldots,C_k)\n", "$$\n", "\n", "over all partitions of $[n]$ of size $k$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "To quote [Wikipedia](https://en.wikipedia.org/wiki/Cluster_analysis#Centroid-based_clustering):\n", "\n", "> In centroid-based clustering, clusters are represented by a central vector, which may not necessarily be a member of the data set. When the number of clusters is fixed to k, k-means clustering gives a formal definition as an optimization problem: find the k cluster centers and assign the objects to the nearest cluster center, such that the squared distances from the cluster are minimized." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "![k-means](https://upload.wikimedia.org/wikipedia/commons/thumb/d/d1/KMeans-density-data.svg/239px-KMeans-density-data.svg.png)\n", "\n", "[(Source)](https://commons.wikimedia.org/wiki/File:KMeans-density-data.svg)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "In general, the problem is NP-hard, that is, no fast algorithm is expected to exist to solve it. The $k$-means algorithm is a popular heuristic. It is based on the idea that the following two sub-problems are easy to solve: \n", "\n", "1. finding the optimal representatives for a fixed partition\n", "2. finding the optimal partition for a fixed set of representatives.\n", "\n", "One then alternates between the two (perhaps until progress falls below a tolerance). To elaborate on the first step, we review an elementary fact about quadratic functions." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "#### 2.1.1 Preliminary result: minimizing a quadratic function" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "A key step of the algorithm involves minimizing a quadratic function. Consider the function\n", "\n", "$$q(x) = a x^2 + b x + c.$$\n", "\n", "When $a > 0$, $q$ has a unique minimum. \n", "\n", "***\n", "**Lemma (Minimum of Quadratic Function):** Let $q(x) = a x^2 + b x + c$ where $a > 0$ and $x \\in \\mathbb{R}$. The unique minimum of $q$ is achieved at $x^* = -b/2a$.\n", "***" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "*Proof:* By the *First-Order Necessary Condition*, a global minimizer of $q$ (which is necessarily a local minimizer) satisfies the condition\n", "\n", "$$\n", "\\frac{\\mathrm{d}}{\\mathrm{d}x} q(x) = 2 a x + b = 0,\n", "$$\n", "\n", "whose unique solution is \n", "\n", "$$\n", "x^*= -\\frac{b}{2a}.\n", "$$ " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "To see that $x^*$ is indeed a global minimizer, we re-write $q$ as\n", "\n", "\n", "\\begin{align}\n", "q(x) \n", "&= a \\left(x^2 + 2 \\left[\\frac{b}{2a}\\right] x\\right) + c\\\\\n", "&= a \\left(x^2 + 2 \\left[\\frac{b}{2a}\\right] x + \\left[\\frac{b}{2a}\\right]^2\\right) - a \\left[\\frac{b}{2a}\\right]^2 + c\\\\\n", "&= a (x - x^*)^2 + \\left[c - \\frac{b^2}{4a}\\right].\n", "\\end{align}\n", " \n", "\n", "Clearly, any other $x$ gives a higher value for $q$. $\\square$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "**NUMERICAL CORNER** Here's a numerical example." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "# Julia version: 1.5.1\n", "using Plots, LinearAlgebra" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "q(x) = a*x^2 + b*x + c;\n", "\n", "x = LinRange(-2, 2, 100); # 100 equally spaced points between -2 and 2\n", "\n", "a, b, c = 2, 4, -1;\n", "y₁ = q.(x); # the . applies the function q element-wise to the 1d array x\n", "plot(x, y₁, label=\"y1\", lw=2)\n", "\n", "a, b, c = 2, -4, 4;\n", "y₂ = q.(x);\n", "plot!(x, y₂, label=\"y2\", lw=2) # the ! indicates that it modifies the existing plot\n", "\n", "a, b, c = -2, 0, 4;\n", "y₃ = q.(x);\n", "plot!(x, y₃, label=\"y3\", lw=2)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "#### 2.1.2 Finding the optimal representatives" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "***\n", "\n", "Lemma (Optimal Representatives): Fix a partition $C_1,\\ldots,C_k$. The optimal representatives under $G$ are the centroids\n", "\n", "$$\n", "\\boldsymbol{\\mu}_i^* = \\frac{1}{|C_i|} \\sum_{j\\in C_i} \\mathbf{x}_j.\n", "$$\n", "\n", "***\n", "\n", "*Proof idea:* The objective $G$ can be written as a sum, where each term is a quadratic function in one component of one of the $\\boldsymbol{\\mu}_i$'s. Each of these terms is minimized by the average of the corresponding components of the $\\mathbf{x}_j$'s belonging $C_i$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "*Proof:* Note that we can expand the $k$-means objective as\n", "\n", "\n", "\\begin{align}\n", "\\sum_{i=1}^k \\sum_{j \\in C_i} \\|\\mathbf{x}_j - \\boldsymbol{\\mu}_i\\|^2\n", "&= \\sum_{i=1}^k \\sum_{j \\in C_i} \\sum_{m=1}^d (x_{j,m} - \\mu_{i,m})^2\\\\\n", "&= \\sum_{m=1}^d \\sum_{i=1}^k \\left[\\sum_{j \\in C_i} (x_{j,m} - \\mu_{i,m})^2\\right].\n", "\\end{align}\n", "\n", "\n", "The expression in square brackets is a quadratic function in $\\mu_{i,m}$:\n", "\n", "$$\n", "q_{i,m}(\\mu_{i,m})\n", "= \\left\\{\\sum_{j \\in C_i} x_{j,m}^2\\right\\} + \\left\\{- 2 \\sum_{j \\in C_i} x_{j,m}\\right\\} \\mu_{i,m} + \\left\\{|C_i| \\right\\} \\mu_{i,m}^2,\n", "$$\n", "\n", "and therefore, by the formula for the minimum of a quadratic function, \n", "is minimized at \n", "\n", "$$\n", "\\mu_{i,m}^* \n", "= - \\frac{- 2 \\sum_{j \\in C_i} x_{j,m}}{2 |C_i|}\n", "= \\frac{1}{|C_i|} \\sum_{j \\in C_i} x_{j,m}.\n", "$$\n", "\n", "Since each term in the sum making up $G$ is strictly minimized at $\\boldsymbol{\\mu}_1^*,\\ldots, \\boldsymbol{\\mu}_k^*$, so is $G$. $\\square$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "#### 2.1.3 Finding the optimal partition" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "***\n", "\n", "Lemma (Optimal Clustering): Fix the representatives $\\boldsymbol{\\mu}_1,\\ldots,\\boldsymbol{\\mu}_k$. An optimal partition under $G$ is obtained as follows. For each $j$, find the $\\boldsymbol{\\mu}_i$ that minimizes $\\|\\mathbf{x}_j - \\boldsymbol{\\mu}_i\\|^2$ (picking one arbitrarily in the case of ties) and assign $\\mathbf{x}_j$ to $C_i$. \n", "\n", "***\n", "\n", "*Proof:* By definition, each term in\n", "\n", "$$\n", "G(C_1,\\ldots,C_k; \\boldsymbol{\\mu}_1, \\ldots, \\boldsymbol{\\mu}_k)\n", "= \\sum_{i=1}^n \\|\\mathbf{x}_i - \\boldsymbol{\\mu}_{c(i)}\\|^2\n", "$$\n", "\n", "is, when the $\\boldsymbol{\\mu}_j$'s are fixed, minimized by the assignment in the statement. Hence so is $G$. $\\square$\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### 2.2 The $k$-means algorithm\n", "\n", "We are now ready to describe the $k$-means algorithm. We start from a random assignment of clusters. We then alternate between the optimal choices in the lemmas. In lieu of pseudo-code, we write out the algorithm in Julia. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "The input X is assumed to be a collection of $n$ vectors $\\mathbf{x_1}, \\ldots, \\mathbf{x}_n \\in \\mathbb{R}^d$ stacked into a matrix. The other input, k, is the desired number of clusters. There is an optional input maxiter for the maximum number of iterations, which is set to $10$ by default.\n", "\n", "We first define separate functions for the two main steps. To find the minimum of an array, we use the function [findmin](https://docs.julialang.org/en/v1/base/collections/#Base.findmin)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "text/plain": [ "opt_reps (generic function with 1 method)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function opt_reps(X, k, assign)\n", " n, d = size(X)\n", " reps = zeros(Float64, k, d) # rows are representatives\n", " for j = 1:k\n", " in_j = [i for i=1:n if assign[i] == j] \n", " reps[j,:] = sum(X[in_j,:],dims=1) ./ length(in_j)\n", " end\n", " return reps\n", "end" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "opt_clust (generic function with 1 method)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function opt_clust(X, k, reps)\n", " n, d = size(X) # n=number of rows, d=number of columns \n", " dist = zeros(Float64, n) # distance to rep\n", " assign = zeros(Int64, n) # cluster assignments\n", " for i = 1:n\n", " dist[i],assign[i] = findmin([norm(X[i,:].-reps[j,:]) for j=1:k])\n", " end\n", " @show G = sum(dist.^2)\n", " return assign\n", "end" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "The main function follows. Below, rand(1:k) is a uniformly chosen integer between 1 and k (inclusive)." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": true, "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "text/plain": [ "mmids_kmeans (generic function with 1 method)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function mmids_kmeans(X, k; maxiter=10)\n", " n, d = size(X)\n", " assign = [rand(1:k) for i=1:n] # start with random assignments\n", " reps = zeros(Int64, k, d) # initialization of reps\n", " for iter = 1:maxiter\n", " # Step 1: Optimal representatives for fixed clusters\n", " reps = opt_reps(X, k, assign) \n", " # Step 2: Optimal clusters for fixed representatives\n", " assign = opt_clust(X, k, reps) \n", " end\n", " return assign\n", "end" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "The $k$-means algorithm is only a heuristic. In particular, it is not guaranteed to find the global minimum of the $k$-means objective. However, it is guaranteed to improve the objective at every iteration, or more precisely, not to make it worse.\n", "\n", "***\n", "\n", "**Theorem (Convergence of $k$-means):** The sequence of objective function values produced by the $k$-means algorithm is non-increasing.\n", "\n", "***\n", "\n", "*Proof idea:* By the *Optimal Representatives* and *Optimal Clustering* lemmas, each step does not increase the objective." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "*Proof:* Let $C_1',\\ldots,C_k'$ be the current clusters, with representatives $\\boldsymbol{\\mu}_1',\\ldots,\\boldsymbol{\\mu}_k'$. After Step 1, the new representatives are $\\boldsymbol{\\mu}_1'',\\ldots,\\boldsymbol{\\mu}_k''$. By the *Optimal Representatives Lemma*, they satisfy\n", "\n", "$$\n", "\\sum_{i=1}^k \\sum_{j \\in C_i'} \\|\\mathbf{x}_j - \\boldsymbol{\\mu}_i''\\|^2\n", "\\leq \\sum_{i=1}^k \\sum_{j \\in C_i'} \\|\\mathbf{x}_j - \\boldsymbol{\\mu}_i'\\|^2.\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "After Step 2, the new clusters are $C_1'',\\ldots,C_k''$. By the *Optimal Clustering Lemma*, they satisfy\n", "\n", "$$\n", "\\sum_{i=1}^k \\sum_{j \\in C_i''} \\|\\mathbf{x}_j - \\boldsymbol{\\mu}_i''\\|^2\n", "\\leq \\sum_{i=1}^k \\sum_{j \\in C_i'} \\|\\mathbf{x}_j - \\boldsymbol{\\mu}_i''\\|^2.\n", "$$\n", "\n", "Combining these two inequalities gives\n", "\n", "$$\n", "\\sum_{i=1}^k \\sum_{j \\in C_i''} \\|\\mathbf{x}_j - \\boldsymbol{\\mu}_i''\\|^2\n", "\\leq \\sum_{i=1}^k \\sum_{j \\in C_i'} \\|\\mathbf{x}_j - \\boldsymbol{\\mu}_i'\\|^2,\n", "$$\n", "\n", "as claimed. $\\square$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "The sequence of objective values is monotone and bounded from below by $0$. [Hence it converges](https://en.wikipedia.org/wiki/Monotone_convergence_theorem#Convergence_of_a_monotone_sequence_of_real_numbers)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "*Exercise:* Modify mmids_kmeans to take a tolerance tol as input and stop when the improvement in objective value G falls below the tolerance. $\\lhd$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## 2.3 Simulated data\n", "\n", "We will test our implementation of $k$-means on a simple simulated dataset." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Recall that a standard Normal variable $X$ has PDF\n", "\n", "$$\n", "f_X(x) \n", "= \\frac{1}{\\sqrt{2 \\pi}}\n", "\\exp\\left(\n", "- x^2/2\n", "\\right).\n", "$$\n", "\n", "Its mean is $0$ and its variance is $1$. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "This is what its PDF looks like:\n", "\n", "![Normal PDF](https://upload.wikimedia.org/wikipedia/commons/thumb/4/44/Standard_Normal_Distribution.png/640px-Standard_Normal_Distribution.png)\n", "([Source](https://upload.wikimedia.org/wikipedia/commons/thumb/4/44/Standard_Normal_Distribution.png/1599px-Standard_Normal_Distribution.png))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "To construct a $d$-dimensional version, we take $d$ independent standard Normal variables $X_1, X_2, \\ldots, X_d$ and form the vector $\\mathbf{X} = (X_1,\\ldots,X_d)$. We will say that $\\mathbf{X}$ is a standard Normal $d$-vector. By independence, its joint PDF is given by the product of the PDFs of the $X_i$'s, that is,\n", "\n", "\n", "\\begin{align}\n", "f_{\\mathbf{X}}(\\mathbf{x})\n", "&= \\prod_{i=1}^d \\frac{1}{\\sqrt{2 \\pi}}\n", "\\exp\\left(\n", "- x_i^2/2\n", "\\right)\\\\\n", "&= \\frac{1}{\\prod_{i=1}^d \\sqrt{2 \\pi}}\n", "\\exp\\left(\n", "- \\sum_{i=1}^d x_i^2/2\n", "\\right)\\\\\n", "&= \\frac{1}{(2 \\pi)^{d/2}} \n", "\\exp(-\\|\\mathbf{x}\\|^2/2).\n", "\\end{align}\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "We can also shift and scale it.\n", "\n", "**Definition (Spherical Gaussian):** Let $\\mathbf{Z}$ be a standard Normal $d$-vector,\n", "let $\\boldsymbol{\\mu} \\in \\mathbb{R}^d$ and let $\\sigma \\in \\mathbb{R}_+$. Then we will refer to the transformed random variable $\\mathbf{X} = \\boldsymbol{\\mu} + \\sigma \\mathbf{Z}$ as a spherical Gaussian with mean $\\boldsymbol{\\mu}$ and variance $\\sigma^2$. We use the notation $\\mathbf{Z} \\sim N_d(\\boldsymbol{\\mu}, \\sigma^2 I)$. $\\lhd$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "The following function generates $n$ data points from two spherical $d$-dimensional Gaussians with variance $1$, one with mean $-w\\mathbf{e}_1$ and one with mean $w \\mathbf{e}_1$. Below, randn(d) generates a d-dimensional spherical Gaussian." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "two_clusters (generic function with 1 method)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function two_clusters(d, n, w)\n", " X1 = reduce(hcat, [vcat(-w, zeros(d-1)) .+ randn(d) for i=1:n])'\n", " X2 = reduce(hcat, [vcat(w, zeros(d-1)) .+ randn(d) for i=1:n])'\n", " return X1, X2\n", "end" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "We will mix these two datasets to form an interesting case for clustering." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### 2.3.1 Two dimension\n", "\n", "We start with $d=2$." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d, n, w = 2, 100, 3.\n", "X1, X2 = two_clusters(d, n, w)\n", "X = vcat(X1, X2)\n", "scatter(X[:,1], X[:,2], legend=false, size=(700,350))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "In the figure above, we observe two fairly well-separated clusters. Let's run $k$-means on this dataset using $k=2$. " ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "G = sum(dist .^ 2) = 1642.660544012831\n", "G = sum(dist .^ 2) = 412.35579408034704\n", "G = sum(dist .^ 2) = 412.35579408034704\n", "G = sum(dist .^ 2) = 412.35579408034704\n", "G = sum(dist .^ 2) = 412.35579408034704\n", "G = sum(dist .^ 2) = 412.35579408034704\n", "G = sum(dist .^ 2) = 412.35579408034704\n", "G = sum(dist .^ 2) = 412.35579408034704\n", "G = sum(dist .^ 2) = 412.35579408034704\n", "G = sum(dist .^ 2) = 412.35579408034704\n" ] } ], "source": [ "assign = mmids_kmeans(X, 2);" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "Our default of $10$ iterations seem to have been enough for the algorithm to converge. We can visualize the result by [coloring](https://docs.juliaplots.org/latest/attributes/) the points according to the assignment. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scatter(X[:,1], X[:,2], marker_z=assign, legend=false, size=(700,350))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### 2.3.2 General dimension\n", "\n", "Let's see what happens in higher dimension. We repeat our experiment with $d=1000$." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d, n, w = 1000, 100, 3\n", "X1, X2 = two_clusters(d, n, w)\n", "X = vcat(X1, X2)\n", "scatter(X[:,1], X[:,2], legend=false, size=(700,350))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "This dataset is in $1000$ dimensions, but we've plotted the data in only the first two dimensions. We observe well-separated clusters. If we plot in any two dimensions not including the first one instead, we see only one cluster. " ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scatter(X[:,2], X[:,3], legend=false, size=(350,350))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "Let's see how $k$-means fares on this dataset." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "G = sum(dist .^ 2) = 199775.6175486944\n", "G = sum(dist .^ 2) = 199754.31646446598\n", "G = sum(dist .^ 2) = 199717.27042900465\n", "G = sum(dist .^ 2) = 199717.27042900465\n", "G = sum(dist .^ 2) = 199717.27042900465\n", "G = sum(dist .^ 2) = 199717.27042900465\n", "G = sum(dist .^ 2) = 199717.27042900465\n", "G = sum(dist .^ 2) = 199717.27042900465\n", "G = sum(dist .^ 2) = 199717.27042900465\n", "G = sum(dist .^ 2) = 199717.27042900465\n" ] } ], "source": [ "assign = mmids_kmeans(X, 2);" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "scatter(X[:,1], X[:,2], marker_z=assign, legend=false, size=(700,350))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Our attempt at clustering does not appear to have been successful. What happened? While these clusters are easy to tease apart if we know to look at the first coordinate only, in the full space the within-cluster and between-cluster distances become harder to distinguish: the noise overwhelms the signal. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "The function below plots the histograms of within-cluster and between-cluster distances for a sample of size $n$ in $d$ dimensions with a given offset. As $d$ increases, the two distributions become increasingly indistinguishable. Later in the course, we will develop dimension-reduction techniques that help deal with this issue." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "data": { "text/plain": [ "highdim_2clusters (generic function with 1 method)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function highdim_2clusters(d, n, w)\n", " # generate datasets\n", " X1, X2 = two_clusters(d, n, w)\n", " \n", " # within-cluster distances for X1\n", " intra = vec([norm(X1[i,:] .- X1[j,:]) for i=1:n, j=1:n if j>i])\n", " h = histogram(intra, normalize=:probability,\n", " label=\"within-cluster\", alpha=0.75, title=\"dim=$d\") # alpha=transparency\n", " \n", " # between-cluster distances\n", " inter = vec([norm(X1[i,:] .- X2[j,:]) for i=1:n, j=1:n])\n", " histogram!(inter, normalize=:probability,\n", " label=\"between-cluster\", alpha=0.75)\n", "end" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "Below we plot the results for dimensions$d=3, 100, 1000$. What do you observe?" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "highdim_2clusters(3, 100, 3)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "highdim_2clusters(100, 100, 3)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "highdim_2clusters(1000, 100, 3)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "*Exercise:* Let$Z_1 \\sim \\mathrm{N}(\\mu_1, \\sigma_1^2)$and$Z_2 \\sim \\mathrm{N}(\\mu_2, \\sigma_2^2)$be independent Normal variables with mean$\\mu_1$,$\\mu_2$and variance$\\sigma_1^2$and$\\sigma_2^2$respectively. Recall that$Z_1 + Z_2$is still Normal.\n", "\n", "(a) What are the mean and variance of$Z_1 + Z_2$?\n", "\n", "(b) Compute$\\mathbb{E}[\\|\\mathbf{X}_1 - \\mathbf{X}_2\\|^2]$and$\\mathbb{E}[\\|\\mathbf{Y}_1 - \\mathbf{Y}_2\\|^2]$explicitly in the claim above.$\\lhd$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "### 2.3.3 Proof of the claim [optional]\n", "\n", "In this optional section, we give a formal statement of the phenomenon described in the previous subsection." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "***\n", "**Claim:** Let$\\mathbf{X}_1, \\mathbf{X}_2, \\mathbf{Y}_1$be independent spherical$d$-dimensional Gaussians with mean$-w_d \\mathbf{e}_1$and variance$1$, where$\\{w_d\\}$is a monotone sequence in$d$. Let$\\mathbf{Y}_2$be an indepedent spherical$d$-dimensional Gaussian with mean$w_d \\mathbf{e}_1$and variance$1$. Then, letting$\\Delta_d = \\|\\mathbf{Y}_1 - \\mathbf{Y}_2\\|^2 - \\|\\mathbf{X}_1 - \\mathbf{X}_2\\|^2$, as$d \\to +\\infty$\n", "\n", "$$\n", "\\frac{\\mathbb{E}[\\Delta_d]}{\\sqrt{\\mathrm{Var}[\\Delta_d]}} \\to\n", "\\begin{cases}\n", "0, & \\text{if w_d \\ll d^{1/4}}\\\\\n", "+\\infty, & \\text{if w_d \\gg d^{1/4}}\n", "\\end{cases}\n", "$$\n", "\n", "where$w_d \\ll d^{1/4}$means$w_d/d^{1/4} \\to 0$.\n", "\n", "***\n", "\n", "The ratio is the statement is referred to as the [signal-to-noise ratio](https://en.wikipedia.org/wiki/Signal-to-noise_ratio). The implication of this theorem is easier to understand in a simulation. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "We will need the following lemmas. Recall first that the covariance of real-valued random variables$W_1$and$W_2is \n", "\n", "\n", "\\begin{align}\n", "\\mathrm{Cov}[W_1,W_2] \n", "&= \\mathbb{E}[(W_1 - \\mathbb{E}[W_1])(W_2 - \\mathbb{E}[W_2])]\\\\\n", "&= \\mathbb{E}[W_1 W_2] - \\mathbb{E}[W_1] \\mathbb{E}[W_2].\n", "\\end{align}\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "***\n", "**Lemma (Variance of a Sum):** IfW_1$and$W_2are real-valued random variables, then\n", "\n", "$$\n", "\\mathrm{Var}[W_1 + W_2]\n", "= \\mathrm{Var}[W_1] + \\mathrm{Var}[W_2] + 2\\,\\mathrm{Cov}[W_1,W_2].\n", "$$\n", "\n", "***\n", "\n", "*Proof:* By definition of the variance and linearity of expectation,\n", "\n", "\n", "\\begin{align}\n", "\\mathrm{Var}[W_1 + W_2] \n", "&= \\mathbb{E}[(W_1 + W_2 - \\mathbb{E}[W_1 + W_2])^2]\\\\\n", "&= \\mathbb{E}[(\\{W_1 - \\mathbb{E}[W_1]\\} + \\{W_2 - \\mathbb{E}[W_2]\\})^2].\n", "\\end{align}\n", "\n", "\n", "Expanding the square gives the claim.\\square$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "***\n", "**Lemma:** Let$W$be a real-valued random variable symmetric about zero, that is, such that$W$and$-W$are identically distributed. Then for all odd$k$,$\\mathbb{E}[W^k] = 0$. \n", "***\n", "\n", "*Proof:* By the symmetry, \n", "\n", "$$\n", "\\mathbb{E}[W^k] = \\mathbb{E}[(-W)^k] = \\mathbb{E}[(-1)^k W^k] = - \\mathbb{E}[W^k].\n", "$$\n", "\n", "The only way to satisfy this equation is to have$\\mathbb{E}[W^k] = 0$.$\\square$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "*Proof idea (Claim):* The only coordinate contributing to$\\mathbb{E}[\\Delta_d]$is the first one by linearity of expectation, while all coordinates contribute to$\\mathrm{Var}[\\Delta_d]$. More specifically, a calculation shows that the former is$c_0 w^2$while the latter is$c_1 w^2 + c_2 d$, where$c_0, c_1, c_2$are constants.\n", "\n", "*Proof (Claim):* Write$w := w_d$and$\\Delta := \\Delta_d$to simplify the notation. There are two steps:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*(1) Expectation of$\\Delta$:* By defintion, the random variables$X_{1,i} - X_{2,i}$,$i = 1,\\ldots, d$, and$Y_{1,i} - Y_{2,i}$,$i = 2,\\ldots, d, are identically distributed. So, by linearity of expectation, \n", "\n", "\n", "\\begin{align}\n", "\\mathbb{E}[\\Delta]\n", "&= \\sum_{i=1}^d \\mathbb{E}[(Y_{1,i} - Y_{2,i})^2] - \\sum_{i=1}^d \\mathbb{E}[(X_{1,i} - X_{2,i})^2]\\\\\n", "&= \\mathbb{E}[(Y_{1,1} - Y_{2,1})^2] - \\mathbb{E}[(X_{1,1} - X_{2,1})^2].\n", "\\end{align}\n", "\n", "\n", "Further, we can writeY_{1,1} - Y_{1,2} \\sim (Z_1 -w) - (Z_2+w)$where$Z_1, Z_2 \\sim N(0,1)$are independent, where here$\\simindicates equality in distribution. Hence, we have \n", "\n", "\n", "\\begin{align}\n", "\\mathbb{E}[(Y_{1,1} - Y_{2,1})^2]\n", "&= \\mathbb{E}[(Z_1 - Z_2 - 2w)^2]\\\\\n", "&= \\mathbb{E}[(Z_1 - Z_2)^2] \n", "- 4w \\,\\mathbb{E}[Z_1 - Z_2]\n", "+ 4 w^2.\n", "\\end{align}\n", "\n", "\n", "Similarly,X_{1,1} - X_{1,2} \\sim Z_1 - Z_2$so$\\mathbb{E}[(X_{1,1} - X_{2,1})^2] = \\mathbb{E}[(Z_1 - Z_2)^2]$. Since$\\mathbb{E}[Z_1 - Z_2] = 0$, we finally get$\\mathbb{E}[\\Delta] = 4 w^2$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "*(2) Variance of$\\Delta:* Using the observations from (1) and the independence of the coordinates we get\n", "\n", "\n", "\\begin{align}\n", "\\mathrm{Var}[\\Delta]\n", "&= \\sum_{i=1}^d \\mathrm{Var}[(Y_{1,i} - Y_{2,i})^2] + \\sum_{i=1}^d \\mathrm{Var}[(X_{1,i} - X_{2,i})^2]\\\\\n", "&= \\mathrm{Var}[(Z_1 - Z_2 - 2w)^2] \n", "+ (2d-1) \\,\\mathrm{Var}[(Z_1 - Z_2)^2].\n", "\\end{align}\n", "\n", "\n", "By the *Variance of a Sum*, \n", "\n", "\n", "\\begin{align}\n", "\\mathrm{Var}[(Z_1 - Z_2 - 2w)^2]\n", "&= \\mathrm{Var}[(Z_1 - Z_2)^2 - 4w(Z_1 - Z_2) + 4w^2]\\\\ \n", "&= \\mathrm{Var}[(Z_1 - Z_2)^2 - 4w(Z_1 - Z_2)]\\\\\n", "&= \\mathrm{Var}[(Z_1 - Z_2)^2] \n", "+ 16 w^2 \\mathrm{Var}[Z_1 - Z_2]\\\\\n", "&\\quad - 8w \\,\\mathrm{Cov}[(Z_1 - Z_2)^2, Z_1 - Z_2].\n", "\\end{align}\n", "\n", "\n", "BecauseZ_1$and$Z_2$are independent,$\\mathrm{Var}[Z_1 - Z_2]\n", "= \\mathrm{Var}[Z_1] + \\mathrm{Var}[Z_2] = 2$. Moreover, the random variable$(Z_1 - Z_2)is symmetric, so\n", "\n", "\n", "\\begin{align}\n", "\\mathrm{Cov}[(Z_1 - Z_2)^2, Z_1 - Z_2]\n", "&= \\mathbb{E}[(Z_1 - Z_2)^3]\\\\ \n", "& \\quad - \\mathbb{E}[(Z_1 - Z_2)^2] \\,\\mathbb{E}[Z_1 - Z_2]\\\\\n", "&= 0.\n", "\\end{align}\n", "\n", "\n", "Finally,\n", "\n", "$$\n", "\\mathrm{Var}[\\Delta]\n", "= 32 w^2 \n", "+ 2d \\,\\mathrm{Var}[(Z_1 - Z_2)^2]\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "*Putting everything together:*\n", "\n", "$$\n", "\\frac{\\mathbb{E}[\\Delta]}{\\sqrt{\\mathrm{Var}[\\Delta]}}\n", "=\n", "\\frac{4 w^2}{\\sqrt{32 w^2 \n", "+ 2d \\,\\mathrm{Var}[(Z_1 - Z_2)^2]}}.\n", "$$\n", "\n", "Takingd \\to +\\infty$gives the claim.$\\square\$" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.5.1", "language": "julia", "name": "julia-1.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.5.1" } }, "nbformat": 4, "nbformat_minor": 2 }