Intro
History
1950, Turing test.
1951-1955, Perceptron.
- 1951, Marvin Minsky designed SNARC, the first artificial neural network;
- 1951, Christopher Strachey wrote a checkers program and Dietrich Prinz wrote one for chess;
- 1955, Allen Newell and Herbert Simon created “Logic Theorilst” for
proving theorems.
1956, the birth of AI. Dartmouth conference, John McCarthy coined the name "Artificial Intelligence".
1956-1974, golden years. People were optimistic about AI, thinking that it would be solved in a few years.
- Search algorithms via path elimination;
- NLP: First chatterbot ELIZA;
- Robotics: stack blocks, walk with lower limbs.
1974-1980, first AI winter. People realized that AI was much harder than they thought, and funding was cut.
1980-1987, AI Boom. Neural network revives.
- Expert system: answer questions about a specific domain of knowledge, using logical rules from experts;
- 1982, John Hopfield proposes Hopfield net;
- 1982, Geoffrey Hinton and David Rumelhart proposes backpropagation.
1987-1993, second AI winter. Most AI projects are not that useful.
1993-2011, steady development. However, AI researchers call AI other names. Neural network were dead again. SVM, graphical models, and reinforcement learning were popular.
- 1997, Deep blue from IBM beats world chess champion;
- 2005, Stanford robot won DARPA grand challenge by driving autonomously for 131 miles;
2012-now, deep learning.
- 2012, Hinton’s group use deep learning for imagenet, improved accuracy by 10%, which is seen as the beginning of deep learning era;
- 2013, Deepmind beats human on Atari;
- 2016, Deepmind (AlphaGo) beats human on Go;
- 2018-, Foundation models (CLIP, ChatGPT, Midjourey, etc.)
Framework
Supervised Learning: Given
- Domain set \(\mathcal{X}\), which is the set of all possible inputs;
- Label set \(\mathcal{Y}\), which is the set of all possible outputs (\(\mathcal{Y}=\{0,1\}\) for classification);
- Training set \(S=((x_1,y_1),...,(x_n,y_n))\), where \(x_i\in\mathcal{X}\) and \(y_i\in\mathcal{Y}\). Denote \(X_{\text{train}}=(x_1,x_2,...,x_n)\) and \(Y_{\text{train}}=(y_1,y_2,...,y_n)\).
You want to learn a function (called predictor, hypothesis or classifier) \(f:\mathcal{X}\to\mathcal{Y}\) such that \(f(x_i) \approx y_i\) for all \(i\).
In order to evaluate the performance of the learned function, we use a loss function \(l(f, x_i, y_i)\) that quantifies the error between the predicted value \(f(x_i)\) and the true value \(y_i\).
- For categorical target (classification), letting \(l(f,x_i,y_i)=[f(x_i)\neq y_i]\) is good but not differentiable.
- For real number target (regression), let \(l(f,x_i,y_i)=\text{dist}(f(x_i)-y_i).\)
Then the training loss is the average of all individual losses
Similarly, we have test loss, validation loss, and population loss.
Step by step process:
- Identify the task you want to solve;
- Create a dataset, containing thousands, or millions of examples;
- Define a loss function to evaluate the performance of the learned function;
- Learn a function to minimize the loss function;
- Evaluate the learned function on a test set to see how well it generalizes.
Optimization: find \(f\) to minimize the loss \(L\).
Generalization: the learned function \(f\) should perform well on unseen data, not just the training data.
Sometimes, your function may overfit or underfit:
- Overfitting: Your function had too much representation power. It captures the noise in the training data, leading to poor performance on unseen data.
- Underfitting: Your function does not have enough representation power. It fails to capture the underlying patterns in the training data, leading to poor performance on both training and unseen data.
In a classical view, to avoid overfitting, we should restrict the representation power, which is called regularization. However, in a modern view, there are implicit regularizations to prevent overfitting for neural networks, so that overfitting almost never happens.
Unsupervised learning: Given input \(X=(x_1,x_2,...,x_n)\), learn a function \(f\) such that \(f(x_i)\) captures the underlying structure of the data.
Common unsupervised learning tasks:
- Clustering: Group similar data points together.
- Principle component analysis (PCA) / Dimension reduction: Reduce the dimensionality of the data while preserving its variance.
- Generative model: Learn the underlying distribution of the data to generate new samples that resemble the training data. Usually done by mapping a Gaussian to the target distribution.
- Anomaly detection: Identify data points that are significantly different from the rest of the data.
- ...
In pratice, semi-supervised learning is often used, where a small amount (say 10%) of labeled data is used to guide the learning process.
(Implicitly), we need some assumptions.
- Continuity assumption: Points which are close to each other are more likely to share a label;
- Manifold assumption: The data lie approximately on a manifold of much lower dimension than the input space.
Optimization
- Zero-order methods: only has information about \(f\). Hard to optimize;
- First-order methods: has information about \(f\) and \(\nabla f\);
- Second-order methods: has information about \(f\), \(\nabla f\), and \(\nabla^2 f\). Hessian matrix has size \(O(d^2)\), is time-consuming.
Almost universal optimization algorithm: (Stochastic) Gradient Descent.
Gradient Descent
The intuition is, gradient gives the direction of fastest descent of the function.
where \(\eta_n\) is the step size (learning rate) at iteration \(n\).
Smoothness
We assume that \(f\in \mathscr{C}^2\) in the following.
Consider the first order Taylor expansion of \(f\):
Now we would like to bound the term \(g(\boldsymbol{\xi})\).
Smoothness means that the function is bounded by a quadratic function at every point. We will show that it is equivalent to the Lipschitz continuity of the gradient and the norm of the Hessian matrix being bounded.
- \(f\) is \(L\)-smooth;
- \(\nabla f\) is \(L\)-Lipschitz continuous;
- \(\|\nabla^2 f(\boldsymbol{x})\|\le L\) for all \(\boldsymbol{x}\), where \(\|\cdot\|\) is the Spectral Norm of the matrix.
\(1 \Rightarrow 3\): Consider its remainder term in the Lagrange form
Then \(f\) is \(L\)-smooth iff for all \(\boldsymbol{x},\boldsymbol{y}\),
Suppose \(\boldsymbol{y}=t\boldsymbol{y'}\), where \(\|\boldsymbol{y'}\|=1\). Then we have
Let \(t\to 0^+\), it gives
Notice that \(\nabla^2 f(\boldsymbol{x})\) is self-adjoint, so
\(3\Rightarrow 2\):
\(2\Rightarrow 1\):
Back to GD, we can show that
Thus we should set \(\eta=O\left(\frac{1}{L}\right)\).
Convexity
Similarly, we have
- \(f\) is convex;
- \(f(\boldsymbol{x'})\ge f(\boldsymbol{x})+\nabla f(\boldsymbol{x})^\top(\boldsymbol{x'}-\boldsymbol{x})\) for all \(\boldsymbol{x},\boldsymbol{x'}\);
- \(\nabla^2 f(\boldsymbol{x})\) is positive semi-definite for all \(\boldsymbol{x}\).
Convexity assumes lower bound of Hessian eigenvalues, while smoothness assumes upper bound of Hessian eigenvalues.
Since \(f(\boldsymbol{x_t})\le f(\boldsymbol{x_{t-1}})\le ...\le f(\boldsymbol{x_0})\), we have
Therefore, we only need \(T=\frac{L\|\boldsymbol{x_0}-\boldsymbol{x^*}\|^2}{2\epsilon}\) steps to find a solution \(f(\boldsymbol{x_T})\le f(\boldsymbol{x^*})+\epsilon\).
Thus
We can see that strong convexity ensures a linear convergence rate.
Stochastic Gradient Descent
GD has two limitations:
- Each iteration requires computing the full gradient, which is expensive for large datasets;
- It may get stuck at stationary points (local minima, saddle points) for non-convex functions.
The solution is to add randomness. Stochastic Gradient Descent (SGD) is a variant of GD that only uses a small random subset of the data to compute the gradient at each iteration. Formally,
Usually, \(|S|\) (batch size) is 64, 128, or 256.
- If \(|S|\) is too small, the variance of \(\boldsymbol{G_n}\) is large;
- If \(|S|\) is too large, \(\boldsymbol{G_n}\) is slow to compute.
Thus
Therefore, letting
gives \(\mathbb{E}[f(\boldsymbol{\overline{x_T}})]\le f(\boldsymbol{x^*})+\epsilon\).
We can see that the convergence rate of SGD is \(O\left(\frac{1}{\sqrt{n}}\right)\), which is slower than GD's \(O\left(\frac{1}{n}\right)\) rate. This is because of the variance introduced by the stochastic gradient. Furthermore, in strongly convex case, SGD achieves \(O\left(\frac{1}{n}\right)\) rate.
In practice, noise is not that bad, and we do not care about training loss but population loss. What's more, we can use various tricks to reduce the variance.
Stochastic Variance Reduced Gradient
The main idea is, maintaining some gradient statistics so that the expected gradient is still correct, but the variance is reduced.
- For \(s=1,2,...\)
- \(\widetilde{\boldsymbol{x}}=\boldsymbol{\widetilde{x}_{s-1}}\)
- \(\widetilde{\boldsymbol{u}}=\frac{1}{N}\sum_{i=1}^N\nabla l_i(\widetilde{\boldsymbol{x}})=\nabla f(\widetilde{\boldsymbol{x}})\)
- \(\boldsymbol{x_0}=\widetilde{\boldsymbol{x}}\)
- For \(t=1,2,...,m\)
- Randomly pick a sample \(i_t\) from \(\{1,2,...,N\}\)
- \(\boldsymbol{x_t}=\boldsymbol{x_{t-1}}-\eta(\nabla l_{i_t}(\boldsymbol{x_{t-1}})-\nabla l_{i_t}(\widetilde{\boldsymbol{x}})+\boldsymbol{\widetilde{u}})\)
- Two update rules: \(\boldsymbol{\widetilde{x}_s}=\boldsymbol{x_m}\) or \(\boldsymbol{x_t}\) for randomly chosen \(t\in \{0,1,...,m-1\}\).
In practice, we usually choose the first rule, while the second rule is used for theoretical analysis.
Immediately we have
Taking expectation over \(i\), we have
Denote the update term as \(\boldsymbol{v_t}=\nabla l_{i_t}(\boldsymbol{x_{t-1}})-\nabla l_{i_t}(\boldsymbol{\widetilde{x}})+\nabla f(\widetilde{\boldsymbol{x}})\). Then
and
Next, we consider a small step \(t\),
Summing over \(t=1,2,...,m\), we have
Finally we have
This theorem shows that SVRG has a linear convergence rate in the strongly convex case. Remember that SGD only achieves \(O\left(\frac{1}{n}\right)\) rate in this case.
Linear Coupling
Mirror Descent
Here's another way to explain GD. By Taylor expansion, \(f(\boldsymbol{x_n})+\nabla f(\boldsymbol{x_n})^\top(\boldsymbol{x}-\boldsymbol{x_n})\) can approximate \(f(\boldsymbol{x})\) to a certain extent. However, we can't directly minimize this approximation, so a regularization term \(\frac{1}{2\eta}\|\boldsymbol{x}-\boldsymbol{x_n}\|^2\) is added to ensure the solution not too far from \(\boldsymbol{x_n}\). That is,
which is equivalent to \(\boldsymbol{x_{n+1}}=\boldsymbol{x_n}-\eta\nabla f(\boldsymbol{x_n})\), the update rule of GD.
Now we want to generalize this term.
Assume \(w\) is \(1\)-strictly convex, then \(V_{\boldsymbol{x}}(\boldsymbol{y})\ge \frac{1}{2}\|\boldsymbol{y}-\boldsymbol{x}\|^2\), which is the original term.
We can see that it's actually GD on mirror space. One can verify that the mirror map is \(\nabla w\), so
Thus
Immediatey we have
We can see that MD has a convergence rate of \(O\left(\frac{1}{\sqrt{n}}\right)\).
MD performs well when \(\rho\) (or \(\|\nabla f(\boldsymbol{x})\|\)) is small, while GD performs well when \(\|\nabla f(\boldsymbol{x})\|\) is large (recall lemma). A straightforward idea is to balance these two algorithms.
Linear Coupling
This leads to the linear coupling method. For each step, we run MD and GD in parallel, and then combine the results.
The first term:
The second term:
Combining them gives
Summing over \(n\) gives
Thus
Setting \(T=4\sqrt{\frac{L\Theta}{d}}\) gives
Now we get a new parameter \(d'=\frac{d}{2}\) and update all the hyperparameters. Then do the same thing again with \(T'=4\sqrt{\frac{L\Theta}{d'}}\), and so on. Finally, we want \( f(\boldsymbol{\overline{x_T}})-f(\boldsymbol{x^*})\le \epsilon\). So the total number of iterations is
It achieves a convergence rate of \(O\left(\frac{1}{n^2}\right)\), which is proved to be the best possible under such assumptions.
Non-convex Optimization
Matrix Completion
Given a matrix \(A\) with some entries missing, complete the matrix by filling in the missing entries. Assume
- \(A\) is low-rank;
- The known entries of \(A\) are uniformly distributed;
- (Incoherence) Let the SVD of \(A\) be \(A = U \Sigma V^\top\), where \(U\in\mathbb{R}^{n\times r}, V\in\mathbb{R}^{m\times r},\Sigma\in\mathbb{R}^{r\times r}\). There is a constant \(u\) such that for all \(i\in[n],j\in[m]\),
Approach: Find low-rank \(U\) and \(V\), so that \(UV^\top\) matches the known entries of \(A\). Formally,
where \(\Omega\) is the set of indices of the known entries of \(A\) and \(P_\Omega\) is the projection onto these entries.
Unlike before, this loss function is not convex, one can find multiple local minima. However, we notice that if we fix \(V\), the loss function is convex in \(U\) and vice versa. Therefore, we can use alternating minimization to find a local minimum. That is, for \(t\ge 0,1,...,T-1\), do
So that each subproblem is convex and can be solved efficiently.
Escaping From Saddle Points
In fact, non-convex objective is common in ML, since this assumption: no spurious local minima (i.e., all the local minima are equally good, it suffices to find one of them). Now the only remaining question is how to escape saddle points.
For stationary points \(\nabla L(\boldsymbol{x})=0\):
- If \(\nabla^2 L(\boldsymbol{x})\succ 0\) , then \(\boldsymbol{x}\) is a local minimum;
- If \(\nabla^2 L(\boldsymbol{x})\prec 0\), then \(\boldsymbol{x}\) is a local maximum;
- If \(\nabla^2 L(\boldsymbol{x})\) has both positive and negative eigenvalues, then \(\boldsymbol{x}\) is a strict saddle point;
- If \(\nabla^2 L(\boldsymbol{x})\succeq 0\), then \(\boldsymbol{x}\) is a local minimum or a flat saddle point.
Assume the loss function is strict saddle, i.e., it does not contain any flat saddle points (which is reasonable). We want to show that SGD can (efficietly) escape saddle points for these loss functions.
-
\(\| f(\boldsymbol{x_0})\|\) is big. By smoothness, there exists some constant \(c>0\) such that $$\mathbb{E}\left[f(\boldsymbol{x_1})\right]\le f(\boldsymbol{x_0})-c.$$
-
Otherwise, \(\boldsymbol{x_0}\) is close to a stationary point.
- If \(\boldsymbol{x_0}\) is close to a local minimum, then we will not get out with large probability.
- Otherwise, \(\boldsymbol{x_0}\) is close to a strict saddle point. Then there is a negative direction. Then a random perturbation (in SGD) will give positive projection on the negative direction. We will follow this direction and escape the saddle. Since the Hessian is smooth, the needed number of steps is polynomial.
Generalization
The data is sampled from a distribution \(\mathcal{D}\) on \(\mathcal{X}\times \mathcal{Y}\), which is unknown to the learner. The learner has access to a training set \(S\) of size \(m\), which is sampled from \(\mathcal{D}^m\). The learner can calculate the empirical loss on the training set:
where \(l\) is the loss function.
However, the generalization error of a classifier \(h:\mathcal{X}\to\mathcal{Y}\) is defined as the expected loss with respect to the distribution \(\mathcal{D}\):
where \(l\) is the loss function. For reliability, we assume that there exists a labeling function \(f:\mathcal{X}\to\mathcal{Y}\) such that \(L_\mathcal{D}(f)=0.\)
We would study how to bound the generalization error \(L_\mathcal{D}(h)\) based on the empirical loss \(L_S(h)\).
The No-Free-Lunch Theorem
First, we show that there is no universal learner that can perform well on all possible distributions. In other words, prior knowledge about \(\mathcal{D}\) is necessary for it to be learnable. This is known as the No-Free-Lunch theorem.
- There exists a function \(f:\mathcal{X}\to\{0,1\}\) with \(L_\mathcal{D}(f)=0\);
-
\[\mathbb{E}_{S\sim \mathcal{D}^m}\left[L_{\mathcal{D}}(A(S))\right]\ge\frac{1}{4}. \]
We will show that for every algorithm \(A\), that receives a training set \(S\) of size \(m\) from \(\mathcal{X}\times \{0,1\}\) and returns a function \(A(S):\mathcal{X}\to\{0,1\}\), there exists a distribution \(\mathcal{D}_i\) such that \(\mathbb{E}_{S\sim \mathcal{D}^m}\left[L_{\mathcal{D}}(A(S))\right]\ge\frac{1}{4}.\) Equivalently, we will show that
There are \(k=(2m)^m\) possible sequences of \(m\) examples from \(\mathcal{X}\), denote them by \(S_1,\ldots,S_k\). For each \(S_j=(x_1,\ldots,x_m)\), let \(S^i_j=\left((x_1,f_i(x_1)),\ldots,(x_m,f_i(x_m))\right)\). If the distribution is \(\mathcal{D}_i\), then \(S_1^i,\ldots,S_k^i\) are the possible training sets.
Now we fix some \(j\in[k]\). Denote \(S_j=(x_1,\ldots,x_m)\) and \(\{v_1,...,v_m\}=\mathcal{X}-S_j\). For each function \(h:\mathcal{X}\to\{0,1\}\) and \(i\in[T]\) we have
Thus
Now we fix some \(r\in[p]\). We can partition all the function \(f_1,\ldots,f_T\) into \(\frac{T}{2}\) disjoint pairs, where for a pair \((f_i,f_{i'})\) we have that for every \(x\in\mathcal{X}\), \(f_i(x)\neq f_{i'}(x)\) if and only if \(c=v_r\). Since for such a pair we have that \(S_j^i=S_j^{i'}\), it follows that
Immediately,
Thus
PAC Learning
Empirical Risk Minimization (ERM) means that the learner chooses a hypothesis \(h\) that minimizes the empirical loss on the training set \(S\). But it may lead to overfitting. To avoid this, we may apply the ERM learning rule over a restricted search space. The PAC (Probably Approximately Correct) learning framework formalizes this idea.
We can decompose the generalization error of a hypothesis \(h\):
where \(\epsilon_{\text{app}}:=\min\limits_{h'\in\mathcal{H}}L_\mathcal{D}(h')\) is the approximation error and \(\epsilon_{\text{est}}:=L_\mathcal{D}(h)-\min\limits_{h'\in\mathcal{H}}L_\mathcal{D}(h')\) is the estimation error. The more we restrict the hypothesis class \(\mathcal{H}\), the larger the approximation error \(\epsilon_{\text{app}}\) is, but (maybe) the smaller the estimation error \(\epsilon_{\text{est}}\) is. The goal of PAC learning is to find a hypothesis class balancing these two errors.
i.e., \(\epsilon_{\text{app}}=0\).
We will show that if \(\mathcal{H}\) is finite, then it is PAC learnable by ERM, and the sample complexity \(m_\mathcal{H}(\epsilon,\delta)\le\log(|\mathcal{H}|/\delta)/\epsilon\).
What about infinite hypothesis class? We will show the class of threshold functions is PAC learnable, also by ERM. (A function \(f:\mathbb{R}\to\{0,1\}\) is a threshold function if there exists a threshold \(t\) such that \(f(x)=\mathbb{I}(x< t)\) for all \(x\in\mathbb{R}\).)
Given a training set \(S\), let \(b_0=\max\{x:(x,1)\in S\}\) and \(b_1=\min\{x:(x,0)\in S\}\). Let \(b_S\) be a threshold corresponding to an ERM hypothesis \(h_S\), which implies that \(b_S\in (b_0,b_1]\). Therefore, a sufficient condition for \(L_\mathcal{D}(h_S)\le \epsilon\) is that both \(b_0\ge a_0\) and \(b_1\le a_1\).
Thus
In pratice, realizable assumption may not hold. In this case, we can replace PAC learning by agnostic PAC learning.
VC Dimension
We have shown that threshold functions are PAC learnable. What about other infinite hypothesis classes? Generally, we have VC dimension to do the analysis.
The intuition is, recall No-Free-Lunch theorem, the adversary used a finite set \(\mathcal{C}\subseteq \mathcal{X}\) and concentrated on examples in \(\mathcal{C}\). Thus we should study the behavior of \(\mathcal{H}\) on \(\mathcal{C}\).
Recall the proof of No-Free-Lunch theorem, we have
This leads us directly to the definition of the VC dimension.
A direct consequence of the corollary above is therefore
- \(\mathcal{H}\) is agnostic PAC learnable with sample complexity
- \(\mathcal{H}\) is PAC learnable with sample complexity
Supervised Learning
Linear Methods
Linear Regression
Let \(\mathcal{X}\subset\mathbb{R}^{ d}\) and \(\mathcal{Y}=\mathbb{R}\). Given \(S=((x_1,y_1),...,(x_n,y_n))\), we would like to learn a linear function \(f:\mathbb{R}^d\to\mathbb{R}\) that best approximates the data in \(S\).
By definition, the hypothesis class of linear regression predictors is
Clearly, \(b\) can be ignored, since we can always add one more dimension with a constant value of \(1\) to \(x\), so that \(w^\topx+b=[w,b]^\top[x,1]\).
With the squared loss function $$L(y,x_i,y_i)=(f(x_i)-y_i)^2,$$ the exact minimizer \(w^*\) is given by
where \(X\) is the matrix with rows \(x_i^\top\) and \(Y\) is the vector with entries \(y_i\).
Since the loss function is convex, we can also use gradient descent to find \(w^*\).
Perceptron
What if the task is not regression, but classification? The first attempt is let \(f(x)=\text{sign}(w^\topx)\), where \(\text{sign}(z)=1\) if \(z\geq 0\) and \(-1\) otherwise.
This is called a perceptron, which can be seen as a simple form of neural network. However, the sign function is not differentiable, so we cannot use gradient descent but instead use the perceptron algorithm.
Repeat until convergence (all examples are classified correctly):
- Pick a random example \((x,y)\) from \(S\).
- If \(y=1\) and \(w^\topx<0\), then \(w\leftarrow w+x\).
- Else if \(y=-1\) and \(w^\topx\ge 0\), then \(w\leftarrow w-x\).
Immediately, \(w_{t+1}^\topw^*\ge t\gamma\).
On the other hand, since \(y_ix_i^\topw_t\le 0\), we have $$|w_{t+1}|2=|w_t+y_ix_i|2=|w_t|2+|y_ix_i|2+2y_ix_i^\topw_t\le |w_t|2+R2.$$
Combining these two inequalities gives
So \(t\le \frac{R^2}{\gamma^2}\), which proves the theorem.
Logistic Regression
Rather than using sign function, we can output a probability \(f(x)\in[0,1]\), which means the probability that \(y=-1\) given \(x\). Equivalently, we predict a distribution \((f(x), 1-f(x))\) over the two classes. And the target distribution is \((1,0)\) if \(y=-1\) and \((0,1)\) if \(y=1\). Finally, we can use L1, cross-entropy or other loss functions. The above approach is called logistic regression.
The hypothesis class associated with logistic regression is the composition of a logistic function over the class of linear functions, where the logistic function is generally the sigmoid function \(\sigma_{\text{sig}}:\mathbb{R}\to[0,1]\) defined as
Cross-entropy is a widely used way to compute loss between two distributions. Given a distribution \(p\) and a target distribution \(y\), the cross-entropy is defined as
In particular, for binary classification, we have
Ridge Regression
To avoid overfitting, we can add a regularization term to the loss function. For example, in linear regression, we can use the following loss function:
where \(\lambda>0\) is a hyperparameter that controls the strength of regularization. The second term is called L2 regularization or ridge regularization. From PAC learning perspective, it restricts the hypothesis class to \(\|w\|_2^2\le c\) for some \(c>0\).
Its gradient and Hessian are given by
Note that it is \(\lambda\)-strongly convex, so GD converges efficiently. In pratice, the update step is divided into two parts:
- \(\widehat{w_{t+1}}=w_t-\frac{\eta}{n}\sum_{i=1}^n (w_t^\topx_i-y_i)x_i\) (GD)
- \(w_{t+1}=(1-\eta\lambda)\widehat{w_{t+1}}\) (weight decay)
and \(\lambda\) is chosen by validation.
Lasso Regression
In some scenarios, we want to restrict the number of non-zero entries in \(w\), which is called feature selection. Unfortunately, \(\|w\|_0\) is not convex, so we use L1 regularization instead. Similar to ridge regression, the loss function is given by
where \(\|w\|_1=\sum_{j=1}^d |w_j|\). This is called Lasso.
Now it becomes convex. However, it is not differentiable, since
In practice, we do the following update:
- \(\widehat{w_{t+1}}=w_t-\frac{\eta}{n}\sum_{i=1}^n (w_t^\topx_i-y_i)x_i\)
- For every coordinate \(i\):
- If \(\widehat{w_{t+1}}_i> \eta\lambda\), then \(w_{t+1,i}=\widehat{w_{t+1}}_i-\eta\lambda\).
- Else if \(\widehat{w_{t+1}}_i< -\eta\lambda\), then \(w_{t+1,i}=\widehat{w_{t+1}}_i+\eta\lambda\).
- Else \(w_{t+1,i}=0\).
Compressed Sensing
Let \(\boldsymbol{x}\in\mathbb{R}^d\) be a "compressible" vector. We may pick a non-adaptive matrix \(A\in\mathbb{R}^{n\times d}\), where \(n\ll d\), as the measurement matrix, and obtain the measurement vector \(\boldsymbol{y}=A\boldsymbol{x}\).
The question is how to choose \(A\) such that we can recover \(\boldsymbol{x}\) from \(\boldsymbol{y}\).
Since \(A(\boldsymbol{\tilde{x}}-\boldsymbol{x})=0\), we have \(\|\boldsymbol{\tilde{x}}-\boldsymbol{x}\|_2^2=0\), which implies \(\boldsymbol{\tilde{x}}=\boldsymbol{x}\).
However, \(\arg\min_{\boldsymbol{z}:A\boldsymbol{z}=\boldsymbol{y}}\|\boldsymbol{z}\|_0\) is hard to optimize. Instead, we can use L1 minimization.
For simplicity, we will use the following notation below: Given a vector \(\boldsymbol{v}\) and a set of indices \(I\), we denote \(\boldsymbol{v_I}\) as the vector obtained by keeping only the entries of \(\boldsymbol{v}\) indexed by \(I\) and setting the rest to \(0\).
Let \(\boldsymbol{h}=\boldsymbol{\tilde{x}}-\boldsymbol{x}\). We want to bound
First, let's consider \(\left\|\boldsymbol{h_{T_{0,1}^c}}\right\|_2\).
By triangle inequality, we have
Since \(\boldsymbol{\tilde{x}}=\boldsymbol{x}+\boldsymbol{h}\) has minimal L1 norm,
Immediately gives
Combining these two inequalities gives the lemma:
Next, let's consider \(\left\|\boldsymbol{h_{T_{0,1}}}\right\|_2\).
Since \(A\boldsymbol{h_{T_{0,1}}}=A\boldsymbol{h}-A\boldsymbol{h_{T_{0,1}^c}}=-\sum_{j\ge 2}A\boldsymbol{h_{T_j}}\), we have
So
In Claim 1, we know
Thus
Since \(\left\|\boldsymbol{h_{T_0}}\right\|_2\le \left\|\boldsymbol{h_{T_{0,1}}}\right\|_2\), finally we have
Combining the two claims immediately proves the theorem:
If \(\boldsymbol{x}\) is sparse, then \(\boldsymbol{\tilde{x}}\approx\boldsymbol{x}\):
How can we find a RIP matrix? One way is to use random matrices.
Denote \(B(\boldsymbol{x},\epsilon)=\{\boldsymbol{y}\in\mathbb{R}^d:\|\boldsymbol{x}-\boldsymbol{y}\|\le \epsilon\}\). Note that \(B\left(\boldsymbol{x},\frac{\epsilon}{2}\right)\) are disjoint for different \(\boldsymbol{x}\in Q\), and
Thus, we have
where \(\text{vol}(B)\) is the volume of \(B\).
By Lemma, there exists a set \(Q\) of size \(|Q|\le \left(\frac{20}{\epsilon}\right)^s\) such that
By definition of \(U_I\), immediately gives
Apply JL Lemma on \(\{U_I\boldsymbol{v}:\boldsymbol{v}\in Q\}\), we know that for \(n\ge \frac{4\cdot 6\ln\left(\frac{2|Q|}{\epsilon}\right)}{\epsilon^2}=24\frac{\ln\left(\frac{2}{\delta}\right)+s\ln\left(\frac{20}{\epsilon}\right)}{\epsilon^2}\), with probability at least \(1-\delta\), we have
It implies that
and
Let \(a=\sup_{\boldsymbol{x}\in S}\frac{\|A\boldsymbol{x}\|}{\|\boldsymbol{x}\|}-1\). For any \(\boldsymbol{x}\in S\) of size \(\|\boldsymbol{x}\|=1\), there exists \(\boldsymbol{v}\in Q\) such that
Solving it for \(a\) gives
On the other hand,
So far, we have
Also implies
Finally, let \(\delta'=d^s\delta,\epsilon'=3\epsilon\) and apply union bound over all possible \(I\subseteq[d]\) of size \(s\), we know that with probability at least \(1-\delta'\), and
What if the sparsity is not on a linear basis? That is, the range of nonlinear function \(G\) that maps a low-dimensional space to a high-dimensional space. For example, in rectangle face generation, \(z\in\mathbb{R}^{12}\) draws a rectangle face and \(G(z)\) draws a realistic face based on the rectangle face. How can we compress \(G(z)\) so that we can still recover it from the compressed vector?
Support Vector Machine
Recall that perceptron is a linear classifier that tries to find a hyperplane that separates the data points of different classes. However, there are often many such hyperplanes that can separate the data points. The Support Vector Machine (SVM) aims to find the optimal hyperplane that maximizes the margin between the two classes, which helps in better generalization.
Lagrange Duality
The first question is how to define the margin. A straightforward way is to define the margin as \(\max_i\frac{|\boldsymbol{w}^\top\boldsymbol{x_i}|}{\|\boldsymbol{w}\|_2}\).
Hard SVM: Find \(\arg\min_{\boldsymbol{w}}\frac{\|\boldsymbol{w}\|^2_2}{2}\) subject to \(y_i\boldsymbol{w}^\top\boldsymbol{x}_i \geq 1\) for all \(i\).
Now, we can use Lagrange multipliers to solve this constrained optimization problem.
Its Lagrangian is given by
To solve the dual problem \(\max_{\alpha\ge 0}\min_{\boldsymbol{w}}\mathcal{L}(\boldsymbol{w},\alpha)\), we take the derivative:
Substituting this back into the Lagrangian gives
However, this method assumes that the data is linearly separable. If it is not, we can relax the hard constraint by using hinge loss \(\ell(\boldsymbol{x},y)=\max\left(0,1-y\boldsymbol{w}^\top\boldsymbol{x}\right).\) Equivalently, we can introduce slack variables \(\xi_i\) that allow some misclassification.
Soft SVM: Find \(\arg\min_{\boldsymbol{w}}\left(\frac{\|\boldsymbol{w}\|^2_2}{2} + \lambda\sum_i\xi_i\right)\) subject to \(y_i\boldsymbol{w}^\top\boldsymbol{x}_i \geq 1 - \xi_i\) and \(\xi_i \geq 0\) for all \(i\).
Its Lagrangian is given by
To solve the dual problem \(\max_{\alpha\ge 0,\kappa\ge 0}\min_{\boldsymbol{w},\xi}\mathcal{L}(\boldsymbol{w},\xi,\alpha,\kappa)\), we take the derivative:
Substituting this back into the Lagrangian gives
We can see that these two dual problems are consistent.
Kernel Method
In fact, there are some quadratic programming algorithms that can solve the dual faster than the primal. Then why solve the dual problem instead of the primal one? The main reason is that we can use the so-called "kernel method".
Sometimes the data is not linearly separable but can be separated using higher-order features. So we can map the data to a higher-dimensional space using a feature map \(\phi(\cdot)\), and then apply SVM in that space.
Notice that in the dual problem, the data points only appear in the form of inner products \(\boldsymbol{x_i}^\top\boldsymbol{x_j}\). Therefore, we can replace these inner products with a kernel function \(K(\boldsymbol{x_i},\boldsymbol{x_j})=\phi(\boldsymbol{x_i})^\top\phi(\boldsymbol{x_j})\). So we do not need to explicitly compute the mapping \(\phi(\cdot)\), which may be computationally expensive.
- Pick a kernel function \(K(\cdot,\cdot)\), which corresponds to some feature map \(\phi(\cdot)\).
- Solve the following quadratic program:\[\max_{\forall i,0\le \alpha_i\le \lambda}\left(\sum_i \alpha_i - \frac{1}{2}\sum_{i,j} \alpha_i \alpha_j y_i y_j K(\boldsymbol{x_i},\boldsymbol{x_j})\right). \]
- The classifier is given by \(f(\boldsymbol{x})=\text{sign}\left(\sum_i \alpha_i y_i K(\boldsymbol{x_i},\boldsymbol{x})\right)\).
The next question is how to choose the kernel function so that it corresponds to a valid feature map.
is positive semi-definite for any \(\{\boldsymbol{x_i}\}\), then there exists a feature map \(\phi(\cdot)\) such that \(K(\boldsymbol{x_i},\boldsymbol{x_j})=\phi(\boldsymbol{x_i})^\top\phi(\boldsymbol{x_j})\).
For example, the polynomial kernel \(K(\boldsymbol{x_i},\boldsymbol{x_j})=(\boldsymbol{x_i}^\top\boldsymbol{x_j}+c)^d\) and the Gaussian kernel \(K(\boldsymbol{x_i},\boldsymbol{x_j})=\exp\left(-\frac{\|\boldsymbol{x_i}-\boldsymbol{x_j}\|^2_2}{2\sigma^2}\right)\) are both valid kernels.
Decision Tree
SVM makes predictions by single decision boundary, while decision tree makes predictions by recursively splitting on different features according to a tree structure.
Boolean Functional Analysis
First, we need to define some concepts with respect to boolean functions, i.e., \(f:\{-1,1\}^n\to [0,1]\). Let \(\mathcal{D}\) be the uniform distribution over \(\{-1,1\}^n\) and \(\chi_S(\boldsymbol{x})=\prod_{i\in S} x_i\) for any \(S\subseteq [n]\).
Clearly, \(\{\chi_S\}_{S\subseteq [n]}\) is a random orthonormal family, called the Fourier basis.
Intuitively, the easier the function is to represent, the simpler the Fourier coefficients are, i.e., the lower the degree and higher the sparsity.
We would like to approximate a decision tree by a low-degree and sparse function.
Clearly, a decision tree with \(s\) leaves can be represented by union of \(s\) "AND" terms, denoted by \(f\). Since every "AND" term has \(L_1\le 1\) and at most \(\log\left(\frac{s}{\epsilon}\right)\) variables, so \(L_1(f)\le s\) and the degree of \(f\) is at most \(\log\left(\frac{s}{\epsilon}\right)\).
Finally, let \(h\) be the truncation of \(f\) by keeping only the Fourier coefficients with \(|\widehat{f}_S|\ge\frac{\epsilon}{L_1(f)}\). Then \(L_0(h)\le \frac{L_1(f)}{\frac{\epsilon}{L_1(f)}}\le\frac{s^2}{\epsilon}\). By Parseval's identity, the missing terms have contribution at most
Thus, \(h\) are \(\epsilon\)-close to \(f\). By triangle inequality, we have $$\mathbb{E}{\boldsymbol{x}\sim \mathcal{D}}[T(\boldsymbol{x})- h(\boldsymbol{x})]^2\le 2\mathbb{E}\sim \mathcal{D}}\left[(T(\boldsymbol{x})- f(\boldsymbol{x}))^2+(f(\boldsymbol{x})- h(\boldsymbol{x}))^2\right]\le 4\epsilon.$$
However, in pratice, we do not know the decision tree \(T\) but only have access to random samples. Under this setting, we can use the following methods.
- Take \(m\) uniformly random samples \(\{\boldsymbol{x_i},y_i\}_{i=1}^m\) for \(f\).
- For every \(S\subseteq [n]\) with degree at most \(\log\left(\frac{s}{\epsilon}\right)\), estimate \(\widehat{f}_S\) by \(\frac{1}{m}\sum_{i=1}^m y_i \chi_S(\boldsymbol{x_i})\).
Construct Decision Tree
Now, we consider how to construct a decision tree from a boolean function. Unfortunately, finding the optimal decision tree is \(\mathsf{NP}\)-complete. Therefore, greedy algorithms are often used in practice.
The key step is to decide which feature to split on at each node. To do this, we need to define a metric to measure the "impurity" (that is, how mixed the classes are) of a node. A common choice is the Gini impurity. Gini impurity measures the probability of two randomly sampled elements having the same label. Formally, for a node with \(K\) classes, let \(p_k\) be the proportion of class \(k\) in the node. Then the Gini impurity is defined as
- At each node, compute the Gini impurity of the current node. If it is zero, stop.
- For each feature, compute the Gini impurity of the child nodes after splitting on that feature.
- Choose the feature that minimizes the weighted average Gini impurity of the child nodes.
- Repeat the process for each child node.
Theoretically, if decision tree is super deep, it can fit anything. However, it may overfit. To avoid this, one may force the decision tree to only consider a subset of data. We can use bagging (or bootstrap aggregating) to achieve this. Formally, denote the dataset by \((x_1,y_1),...,(x_n,y_n)\). We sample \(n\) times with replacement from the original dataset to get a new dataset \((x_1',y_1'),...,(x_n',y_n')\). Then we train the decision tree on this new dataset. This procedure will be repeated multiple times to get an ensemble of decision trees, which is called a random forest.
Boosting
If we have lots of weak leaners, can we combine them to get a strong learner? A straightforward approach is stacking weak learners. That is, let \(h_1, h_2, \ldots, h_n\) be a sequence of weak learners, and we train the final learner on \(\{((h_1(x),h_2(x),...,h_n(x)),y)|(x,y)\in S\}\), where \(S\) is the original training set.
In fact, in practice, we have a transparent learning algorithm, which means that the learner on any training set can be adaptively acquired. Therefore, we can iteratively add weak learners to the current model to reduce the loss. Boosting does exactly this.
Gradient Boosting
We can view the problem as an optimization problem in function space. So we can use gradient descent to solve it. That is, at each iteration, we find a weak learner that fits the negative gradient of the loss function with respect to the current model.
For \(t=1\) to \(T\):
- For all \(i\), calculate gradients $$g(x_i)=\frac{\partial L(H)}{\partial H(x_i)}.$$
- Find the best weak learner \(h_t\) that fits the negative gradients, i.e., $$h_t=\arg\min \sum\limits_{i=1}^n\ell(-g(x_i),h(x_i)).$$
- Update the learner \(H\leftarrow H+\eta h_t\) for some step size \(\eta\).
Output the final learner \(H\).
For example, if we use squared loss \(\ell(y,h(x))=\frac{1}{2}(y-h(x))^2\), then \(-g(x_i)=y_i-H(x_i)\), which is exactly the residual.
AdaBoost
Another popular boosting algorithms is AdaBoost. Instead of fitting the negative gradient, AdaBoost reweights the training examples at each iteration, so that harder examples are given more weight.
For \(t=1\) to \(T\):
- Find the weak learner \(h_t\) that minimizes the weighted error $$\epsilon_t=\text{Pr}_{i\sim \mathcal{D}_t}[h_t(x_i)\ne y_i];$$
- Construct distribution \(\mathcal{D}_t\):
where \(\alpha_t=\frac{1}{2}\ln\left(\frac{1-\epsilon_t}{\epsilon_t}\right)\) and \(Z_t\) is the normalization factor.
Output the final learner $$H_{\text{final}}(x)=\text{sign}\left(\sum_{t=1}^\top \alpha_t h_t(x)\right).$$
By the way, one can verify that AdaBoost is a special case of gradient boosting with exponential loss \(\ell(y,h(x))=\exp(-y h(x))\).
The training error of \(H_{\text{final}}\) is
On the other hand,
Let \(\gamma_t=\frac{1}{2}-\epsilon_t\). Then $$Z_t=\sqrt{1-4\gamma_t^2}\le \exp(-2\gamma_t^2).$$
Immediately we have
Margin Theory
Boosting has a very critical and counterintuitive property: even after the training error has dropped to zero, continuing training can still improve the performance of the model on test set.
The explanation is that boosting increases the margin of the training examples. The margin of an example \((x,y)\) is defined as \(yf(x)\). Intuitively, the larger the margin, the more confident the prediction is. So even if the training error is zero, increasing the margin can still improve the generalization.
Substituting \(Z_t\) and \(\alpha_t\) gives
Assume that \(\theta<\frac{1}{2}-\epsilon_t\) for all \(t\). Then the expression inside the parentheses is less than 1. So the probability decreases exponentially with \(T\).
The significance of margin theory is that we can bound the generalization error by the margin distribution on the training set, instead of the training error.
Nearest Neighbor
KNN
Consider a classification task. When a new query \(x\) is given, we predict its label by looking at the labels of its nearest neighbors in the training set. Such algorithms is called the nearest neighbor algorithms. The most common version is the \(k\)-nearest neighbor algorithm.
Find \(k\) nearest neighbors of \(x^*\) in \(S\), denoted by \(\{(x_{i_1},y_{i_1}),...,(x_{i_k},y_{i_k})\}\).
Predict the label of \(x^*\) by majority vote:
Choosing \(k\) is critical. If \(k\) is too small, the model may be too sensitive to noise. If \(k\) is too large, the model may be too smooth and miss important patterns.
\(k\)NN is non-parametric, which means that there are no parameters we can tune to make the performance better. In other words, the model structure determined from the dataset. Thus, some people think kNN belongs to unsupervised learning.
However, exact nearest neighbor search is hard, so we consider approximation. Since "nearest search" is equivalent to "\(R\)-near search" due to binary search, we can relax "\(R\)-near search" to "\(c\)-approximate \(R\)-near search", which means that we only need to find a point within distance \(cR\) if there exists a point within distance \(R\).
LSH and Approximate Nearest Neighbor Search
- If \(\|x-y\|\leq R\), then \(\Pr_{h\in\mathcal{H}}[h(x)=h(y)]\geq P_1\);
- If \(\|x-y\|\geq cR\), then \(\Pr_{h\in\mathcal{H}}[h(x)=h(y)]\leq P_2\).
An LSH family can be used to design an efficient algorithm for approximate near neighbor search.
Construct \(L\) hash tables, where each table \(i\) contains the dataset points hashed using the function \(g_i\).
When a query \(x\) is given, for each \(j=1,2,\cdots,L\):
- Find the bucket \(g_j(x)\) in the \(j\)-th hash table and retrieve all points in that bucket.
- For each of the retrieved point, check if it is within distance \(cR\) from \(x\). If yes, return it.
- Stop as soon as we have checked \(L'\) points.
On the other hand,
Thus
Thus, the probability that \(x^*\) is found in the first \(2L+1\) points is at least \(\frac{1}{2}-\frac{1}{e}\).
Finally, the problem is how to construct LSH families. For \(\ell_2\) space, we can use the following LSH Library:
where \(r\sim\mathcal{N}(0,I_d)\), \(b\sim\text{Unif}[0,w)\) and \(w>0\) is a parameter. Let \(c=\|x-y\|_2\), we have
where \(f_p\) is the PDF of \(\mathcal{N}(0,1)\). It can be shown that \(P(c)\) is a monotonically decreasing function of \(c\).
Metric Learning
Searching nearest neighbor in the original space \(\mathbb{R}^d\) may not be the best choice. We can use neural networks to learn a good feature space so that inside the space we have better nearest neighbor structure. This is called metric learning.
The key point is how to define the loss function.
NCA: $$L=\sum_i\sum\limits_{j:c_j=c_i}p_{i,j}$$
where \(p_{ii}=0\) and for \(i\neq j\),
LMNN: $$L=\sum_{i,j\in N_i,y_k\neq y_i}\max(0,|f(x_i)-f(x_j)|_2-|f(x_i)-f(x_k)|_2+1)$$
where \(N_i\) is the set of target neighbors of \(x_i\).
Unsupervised Learning
PCA
PCA is a method for finding the directions in high-dimensional data that are most informative, i.e., that capture the most variance. Formally, denote the data matrix by \(X=[x_1,x_2,\cdots,x_n]\in\mathbb{R}^{d\times n}\), we want to find
Equivalently, we can write it as
Clearly, the solution is the top eigenvector of \(XX^\top\). The \(k\)-PCA is to find the top \(k\) eigenvectors of \(XX^\top\). To compute PCA, we can use the power method:
For \(t=1,2,\cdots,T\):
- \(b_t = \frac{(XX^\top) b_{t-1}}{\|(XX^\top) b_{t-1}\|}.\)
Return \(b_T\).
To analyze the convergence, let \(v_1,v_2,\cdots,v_d\) be the eigenvectors of \(XX^\top\) with eigenvalues \(\lambda_1\geq \lambda_2\geq \cdots\geq \lambda_d\geq 0\). We can write \(b_0=\sum_{i=1}^d \alpha_i v_i\). Then
where \(Z_t\) is a normalization factor. Thus, \(b_t\) converges to \(v_1\) exponentially fast with rate \(\lambda_2/\lambda_1\). Once we know \(v_1\), we can remove \(v_1\) component and repeat the process to find \(v_2\) and so on.
Clustering
Clustering is the task of grouping a set of objects in such a way that objects in the same group (called a cluster) are more similar to each other than to those in other clusters.
K-means & Lloyd's Algorithm
One of the most popular clustering algorithms is K-means.
Find \(S=\{S_1,S_2,\cdots,S_k\}\) to minimize
where \(\mu_i = \frac{1}{|S_i|}\sum_{x\in S_i} x\) is the center of cluster \(S_i\).
However, it is \(\mathsf{NP}\)-hard to find the optimal solution. Lloyd's algorithm is an heuristic algorithm that converges quickly in practice.
Repeat until convergence:
- Assign each point to the nearest cluster center:\[S_i = \{x_j: i = \arg\min_{l} \|x_j - \mu_l\|^2\}, \quad \forall i=1,2,\cdots,k. \]
- Update each cluster center:\[\mu_i = \frac{1}{|S_i|}\sum_{x\in S_i} x, \quad \forall i=1,2,\cdots,k. \]
The above two algorithms are clustering under \(\ell_p\), but sometimes this is not the best metric.
Spectral Clustering
Assume that we can define "similarity" between any two points \(x_i\) and \(x_j\), denoted by \(A_{ij}\ge 0\). i.e., \(A\) is the adjacency matrix of a weighted graph. Then we can define the degree matrix \(D\) as a diagonal matrix with \(D_{ii}=\sum_{j=1}^n A_{ij}\), and the Laplacian matrix \(L=D-A\).
- The multiplicity of eigenvalue \(0\) of \(L\) equals the number of connected components \(S_1,S_2,...,S_k\) in \(G\);
- The eigenspace of eigenvalue \(0\) is spanned by the indicator vectors \(\mathbf{1}_{S_1},\mathbf{1}_{S_2},...,\mathbf{1}_{S_k}\).
Assume \(S_1,S_2,...,S_k\) are the connected components of \(G\). One can easily verify that \(L\mathbf{1}_{S_i}=0\) for all \(i=1,2,...,k\). Thus, the multiplicity of eigenvalue \(0\) is at least \(k\). On the other hand, if \(Lv=0\), then \(v^\top L v=0\). Since \(A_{ij}\ge 0\), we must have \(v_i=v_j\) for any \((i,j)\) such that \(A_{ij}>0\). Thus, \(v\) must be a linear combination of \(\mathbf{1}_{S_1},\mathbf{1}_{S_2},...,\mathbf{1}_{S_k}\). Therefore, the multiplicity of eigenvalue \(0\) is exactly \(k\).
However, in general, the graph is connected and we cannot directly use the above lemma. Instead, we compute the first \(k\) eigenvectors of \(L\) and use them as features for clustering.
Construct the similarity graph and compute the adjacency matrix \(A\), the degree matrix \(D\) and the Laplacian matrix \(L=D-A\).
Compute the first \(k\) eigenvectors (corresponding to smallest \(k\) eigenvalues) \(v_1,v_2,\cdots,v_k\) of \(L\).
Let \(U=[v_1,v_2,\cdots,v_k]\in\mathbb{R}^{n\times k}\), and let \(y_i\in\mathbb{R}^k\) be the \(i\)-th row of \(V\) for \(i=1,2,\cdots,n\).
Cluster \((y_1,y_2,\cdots,y_n)\) using K-means into \(k\) clusters \(C_1,C_2,\cdots,C_k\).
Output the clusters \(S_i=\{x_j:y_j\in C_i\}\) for \(i=1,2,\cdots,k\).
How to find the smallest eigenvalues? We already know that the power method can find the largest eigenvalue \(\lambda_{max}\) of \(L\). Since \(A\) is positive semi-definite, \(B=A-\lambda_{max}I\) is negative semi-definite, and the largest eigenvalue of \(B\) is \(\lambda_{min}-\lambda_{max}\). Thus, we can use the power method to find the largest eigenvalue of \(B\) and get the smallest eigenvalue of \(A\).
To understand why spectral clustering works, we need to introduce the RatioCut objective function. Given a clustering \(C_1,C_2,...,C_k\), define
where \(\text{cut}(C_i,\bar{C_i}) = \sum_{x\in C_i,y\in \bar{C_i}} A_{xy}\) is the total weight of edges between \(C_i\) and its complement \(\bar{C_i}\). However, minimizing RatioCut is \(\mathsf{NP}\)-hard.
Thus, let \(H_i=(H_{1i},H_{2i},...,H_{ni})\), we have
Immediately,
Now, we can relax \(H\) to be any orthogonal matrix and solve the following optimization problem:
Standard trace minimization tells us we should pick smallest \(k\) eigenvectors of \(L\) as columns of \(H\). This is exactly what spectral clustering does.
Since \(H^\top H=I\), we have \(\sum_{j=1}^n c_{ij}^2=1\) for all \(i\) and \(\sum_{i=1}^k c_{ij}^2\leq 1\) for all \(j\). That is, \((c_{ij}^2)\) forms the first \(k\) rows of a doubly stochastic matrix \(B\). By Birkhoff's theorem, \(B\) is a convex combination of permutation matrices.
Notice that all the feasible \((c_{ij}^2)\) form a compact and convex set, so the minimizer must occur at an extreme point, which is a permutation matrix. Thus, we should pick \(c_{ij}^2=1\) for \(j=i\) and \(0\) otherwise, which corresponds to the smallest \(k\) eigenvalues of \(L\).
