13. Basics of Machine Learning

Chapter PDF · Complete book

Imaging problems often concern the reconstruction or processing of an individual signal or image, whereas machine learning studies patterns across collections of observations. Their goals and performance measures differ, but both use tools such as linear models, regularization, and convex optimization.

13.1 Unsupervised Learning¶

In unsupervised learning, we observe $n$ points $(x_i)_{i=1}^n$.

The aim is to identify structure in these observations, for example to visualize them or group them into clusters.

We assume for simplicity that the observations lie in Euclidean space, $x_i \in \mathbb{R}^p$, where $p$ is the number of features. We store the observations as the rows of a matrix $X\in\mathbb{R}^{n\times p}$.

13.1.1 Dimensionality Reduction and PCA¶

Dimensionality reduction produces a compact representation of the data for visualization or subsequent analysis. By retaining informative features, it can also reduce computational cost and improve prediction.

Principal component analysis (PCA) projects the data orthogonally onto the principal axes of the covariance matrix.

Presentation of the method.

The empirical mean is defined as

and the empirical covariance as

With the centered data matrix $\tilde X \mathrel{:=}X-\mathbb{1}_n\hat m^*$, the covariance is $\hat C=\tilde X^* \tilde X/n$.

Suppose the points $(x_i)_i$ are i.i.d. with a finite second moment, and let $\mathbf{x}$ denote a random variable with their common distribution. The law of large numbers gives the following almost sure limits as $n \rightarrow +\infty$:

Let $\mu$ be the distribution on $\mathbb{R}^p$ of $\mathbf{x}$. Then these limits can also be written as

The PCA basis consists of the right singular vectors of the centered data matrix. Proposition 9.1 gives the reduced singular value decomposition

where $U \in \mathbb{R}^{n \times r}$ and $V \in \mathbb{R}^{p \times r}$ have orthonormal columns, and $r=\mathop{\mathrm{rank}}(\tilde X) \leqslant\min(n,p)$.

Write $V=(v_k)_{k=1}^r$ for the orthonormal columns, which are eigenvectors of $\hat C = V \mathop{\mathrm{diag}}(\sigma^2) V^\top$, with $v_k \in \mathbb{R}^p$. These vectors describe the principal directions of variation in the point cloud $(x_i)_i$ in $\mathbb{R}^p$.

We order the singular values as $\sigma_1 \geqslant\ldots \geqslant\sigma_r$, so the leading components account for the greatest variance.

Figure 13.1 displays an empirical covariance matrix and its spectrum $(\sigma_k^2)_k$ for the Iris dataset (Iris flower data set) introduced by Fisher. It contains 50 samples from each of three iris species: Iris setosa, Iris virginica, and Iris versicolor. The four features are sepal length, sepal width, petal length, and petal width.

To obtain a PCA embedding $x_i \in \mathbb{R}^p \mapsto z_i \in \mathbb{R}^d$ of dimension $d \leqslant p$, project onto the first $d$ right singular vectors, extending the basis if $d>r$:

The corresponding reconstruction is

Thus $\tilde x_i = \mathop{\mathrm{Proj}}_{\tilde T}(x_i)$, where $\tilde T\mathrel{:=}\hat m+ \mathop{\mathrm{Span}}_{k=1}^d(v_k)$ is the affine reconstruction space.

Figure 13.3 shows an example of PCA for 2-D and 3-D visualization.

Optimality analysis.

We prove that PCA minimizes the $\ell^2$ reconstruction error among linear dimensionality-reduction methods.

An optimal affine reconstruction space can be chosen to contain the empirical mean: for any fixed direction space, translating it to the mean cannot increase squared reconstruction error. After centering, we therefore assume $\hat m=0$, hence $X=\tilde X$.

We recall that $X=\sqrt{n}U\mathop{\mathrm{diag}}(\sigma)V^\top$ with observations in rows, and $\hat C=X^\top X/n=V\mathop{\mathrm{diag}}(\sigma_i^2)V^\top$. We assume $1\leqslant k\leqslant p$ and extend $V$ to an orthogonal basis of $\mathbb{R}^p$ when necessary.

For the centered data, compare linear compression maps into dimension $k$ followed by reconstruction:

Although $f(\cdot,S)$ and $f(R,\cdot)$ are convex separately, $f$ is not jointly convex. Alternating minimization in $R$ and $S$ therefore need not reach a global minimizer.

Nevertheless, a global minimizer has a simple explicit form.

The compressor $S^\top$ and decompressor $R=S$ give precisely the PCA maps in (13.3) and (13.4).

We first prove that one can restrict attention to orthogonal projection matrices.

The next lemma bounds the objective by a linear program. We then show that PCA attains this bound, proving optimality.

13.1.2 Clustering and $k$-means¶

Clustering assigns a label $y_i \in \{1,\ldots,k\}$ to each observation $x_i$. Observations with the same label form a cluster.

$k$-means

One approach is to seek compact clusters by minimizing a measure of within-cluster dispersion. In Euclidean space, the $k$-means objective measures distances from observations to their centroids $c=(c_\ell)_{\ell=1}^k$, with $c_\ell \in \mathbb{R}^p$. Related formulations are possible in general metric spaces, although the centroid computations may be more difficult. The optimization problem is

The $k$-means algorithm alternates between updating the labels and updating the centroids, a form of block coordinate minimization.

Initialize the centroids $c$, for example by choosing well-separated observations. We discuss a more systematic initialization below.

For fixed centroids $c$, minimize $y \mapsto \mathcal{E}(y,c)$ by assigning each observation to its nearest centroid:

For fixed labels $y$, minimize $c \mapsto \mathcal{E}(y,c)$ by setting each centroid to the mean of its cluster:

If a cluster becomes empty, its centroid does not contribute to the objective for the current labels. One may leave it in place or reseed it at a data point before the next assignment step. The centroid formula applies only to nonempty clusters.

Each step decreases or preserves $\mathcal{E}$. With a fixed tie-breaking rule that retains existing assignments at ties, and without indefinite reseeding, a changed assignment strictly decreases the objective. There are finitely many label assignments, so the algorithm terminates after finitely many such changes.

Because the objective is nonconvex, $k$-means need not find globally optimal clusters.

A reseeding heuristic moves centroids $c_\ell$ from redundant or low-contribution clusters to regions with high distortion. This can help the iterations escape a poor local configuration.

Figure 13.8 shows an example of $k$-means iterations on the Iris dataset.

$k$-means++

The result of $k$-means depends strongly on initialization. Well-separated initial centers tend to cover the data more effectively than tightly grouped centers.

The randomized $k$-means++ initialization has an approximation guarantee even before Lloyd iterations are applied. In practice, Lloyd iterations usually improve the resulting centers. First choose $c_1$ uniformly among the sample points. After choosing $c_1,\ldots,c_\ell$, select $c_{\ell+1}$ from the sample according to a probability $\pi^{(\ell)}$ on $\{1,\ldots,n\}$ proportional to the squared distance from the nearest selected center

Points far from the previously seeded centers are more likely to be chosen. If every $d_i=0$, the current centers already achieve zero distortion and the procedure can stop.

The following theorem, due to David Arthur and Sergei Vassilvitskii, bounds the expected seeding cost by a logarithmic factor times the optimal cost. Computing a global optimum is NP-hard in general.

Lloyd algorithm and continuous densities.

The $k$-means iterations are also known as Lloyd’s algorithm and are used in vector quantization for compression. The same alternating procedure applies to a probability measure $\mu$ on $\mathbb{R}^p$ with finite second moment.

The energy to minimize becomes

where $(\mathcal{V}_\ell)_\ell$ is a partition of the domain.

Step (13.7) is replaced by the computation of a Voronoi cell

For distinct centroids, the Voronoi cells are polyhedra bounded by perpendicular bisectors. Resolve distance ties consistently to obtain a partition; this matters if $\mu$ assigns positive mass to a boundary. In low dimensions, computational geometry provides efficient algorithms for constructing the cells.

Step (13.8) is then replaced by

The centroid update applies to cells of positive $\mu$-mass; zero-mass cells require a separate convention, as in the discrete algorithm. For a uniform planar density, hexagonal Voronoi patterns are characteristic of high-resolution quantization away from the boundary. They need not be exact finite-$k$ minimizers on a bounded domain.

Figure 13.9 displays two examples of Lloyd iterations on 2-D densities on a square domain.

13.2 Empirical Risk Minimization¶

We first describe empirical risk minimization, a framework shared by regression and classification. For classification, predictions may represent class probabilities rather than labels.

To make fitting $y_i \approx f(x_i)$ computationally tractable and obtain useful predictions, we restrict the admissible functions to a class of controlled complexity. The class may become richer as the sample size $n$ grows, but its complexity must remain compatible with the available data.

13.2.1 Empirical Risk¶

Let $\mathcal{F}_n$ be a class of functions that may depend on the sample size. A common learning procedure is empirical risk minimization (ERM)

The loss function $L: \mathcal{Y}^2 \rightarrow \mathbb{R}^+$ measures prediction error for the chosen task. Classification and regression typically require different loss functions.

We sometimes write $\hat f_n$ to emphasize the estimator’s dependence on the sample size $n$.

13.2.2 Prediction and Consistency¶

For the analysis, assume the training pairs $(x_i,y_i)$ are i.i.d. with distribution $\pi$ on $\mathcal{X}\times \mathcal{Y}$. The corresponding population optimization problem is

The desired consistency property is $\hat f_n \rightarrow \bar f$ as $n \rightarrow +\infty$. One way to express this is through expected prediction error over the training samples $(x_i,y_i)_i$:

Here the expectation includes both the test input $\mathbf{x}$, drawn from the marginal $\pi_\mathcal{X}$, and the $n$ i.i.d. training pairs $(x_i,y_i)\sim\pi$ that determine $\hat f_n$.

The comparison loss $\tilde L$ measures discrepancies between predictions in $\mathcal{Y}$; for example, we may take $\tilde L=L$. One can also study convergence in probability, i.e.

These properties define consistency in expectation and in probability, respectively.

A convergence rate gives an explicit upper bound on the decay of $E_n$ or $E_{\varepsilon,n}$.

For $\tilde L(y,y')=|y-y'|^r$, convergence in expectation implies convergence in probability by Markov’s inequality:

13.2.3 Parametric Approaches and Regularization¶

Instead of specifying $\mathcal{F}_n$ through a hard constraint, we can favor simple or regular functions through a penalty.

A typical way to achieve this is by using a parametric model $y \approx f(x,\beta)$ where $\beta\in \mathcal{B}$ parametrizes the function $f(\cdot,\beta) : \mathcal{X}\rightarrow \mathcal{Y}$. The empirical risk minimization procedure (13.9) now becomes

where $J$ is a regularizer. For example, $J=|\!| \cdot |\!|_2^2$ controls parameter magnitude, while $J=|\!| \cdot |\!|_1$ promotes sparse coefficients and, in a linear feature model, feature selection. Here $\lambda_n>0$ is a regularization parameter, and its asymptotic scaling must balance approximation and estimation; consistency often requires $\lambda_n\to0$ at a controlled rate.

The population counterpart of (13.10) defines the parameter $\bar \beta$. A limiting estimator as $n \rightarrow +\infty$ then takes the form $\bar f = f(\cdot,\bar \beta)$, with $\bar \beta$ satisfying

Prediction vs. estimation risks.

We may also ask how accurately $\hat\beta$ estimates $\bar\beta$, as measured by the parameter error $|\!| \hat\beta-\bar\beta |\!|$ for a chosen norm $|\!| \cdot |\!|$. Parameter estimation is generally more demanding than prediction: different parameters can produce similar predictions, especially when the model is poorly identified.

13.2.4 Validation, Test Sets, and Cross-validation¶

The errors $E_n$ and $E_{\varepsilon,n}$ cannot be evaluated directly because the population optimum $\bar f$ is unknown.

To tune a parameter such as the regularization strength $\lambda$, we estimate the prediction risk $\mathbb{E}(L(\hat f(\mathbf{x}),\mathbf{y}))$ on data held out from training.

Use a second sample $(\bar x_j,\bar y_j)_{j=1}^{\bar n}$, called a validation set. In the statistical model, these observations are i.i.d. with distribution $\pi$ and independent of the training sample. The validation risk is

which converges to $\mathbb{E}(L(\hat f(\mathbf{x}),\mathbf{y}))$ for large $\bar n$.

Minimizing $R_{\bar n}$ over hyperparameters (such as $\lambda_n$) is holdout validation. Cross-validation repeats training and validation across folds. A separate test set, untouched by hyperparameter selection, estimates the performance of the final selected procedure.

13.3 Supervised Learning: Regression¶

In supervised learning, the training data consist of pairs $(x_i,y_i)\in\mathcal{X}\times\mathcal{Y}$, with $\mathcal{X}=\mathbb{R}^p$ here for simplicity. We seek a function $f:\mathcal{X}\rightarrow\mathcal{Y}$ that captures the relationship $y_i\approx f(x_i)$ and predicts an output $f(x)$ for a new input $x$.

When $\mathcal{Y}$ is finite and discrete, the task is supervised classification, studied in Section 13.4. Binary classification uses $\mathcal{Y}=\{0,1\}$; for example, in a medical application, $y_i=0$ may indicate a healthy subject and $y_i=1$ a subject with the condition of interest.

When $\mathcal{Y}$ is continuous, typically $\mathcal{Y}=\mathbb{R}$, the task is regression.

13.3.1 Linear Regression¶

For regression, we specialize empirical risk minimization to $\mathcal{Y}=\mathbb{R}$ with the quadratic loss $L(y,y')=\frac{1}{2}|y-y'|^2$.

Nonlinear regression can be formulated using a dictionary of features, such as polynomials. This lifts the data into a higher-dimensional space and leads to the kernel methods of Section 13.5.

Least squares and conditional expectation.

If $f$ ranges over measurable functions and $\mathbb{E}(\mathbf{y}^2)<\infty$, a population minimizer $\bar f$ in (13.10) is the conditional mean of the response given the input. It is defined up to equality almost everywhere under the input distribution.

Suppose for simplicity that $\pi$ has density $\frac{\mathrm{d}\pi}{\mathrm{d}x \mathrm{d}y}$ with respect to a product measure $\mathrm{d}x \mathrm{d}y$, such as Lebesgue measure. Then

where $(\mathbf{x},\mathbf{y})$ has law $\pi$ and the formula applies where the denominator is positive.

If $\mathcal{X}$ and $\mathcal{Y}$ are discrete, let $\pi_{x,y}$ denote the probability of $(\mathbf{x}=x,\mathbf{y}=y)$. Then

with no prescribed value when the marginal of $\pi$ on $\mathcal{X}$ vanishes at $x$.

Replacing $\pi$ by the empirical distribution averages responses at each observed input but leaves predictions unspecified elsewhere. Generalization therefore requires assumptions on regularity or model complexity.

Penalized linear models.

A linear model constrains $f$ to have the form $f(x,\beta)=\langle x,\,\beta\rangle$ with parameters $\beta\in \mathcal{B}=\mathbb{R}^p$. Affine predictors are included through the identity $\langle x,\,\beta\rangle+\beta_0=\langle (x,1),\,(\beta,\beta_0)\rangle$, by appending a constant coordinate to $x$ to form $(x,1)$. We therefore use the linear notation below.

Under the square loss, the regularized ERM (13.11) is conveniently rewritten as

where we introduced the empirical second-moment matrix (equal to the covariance (13.1) for centered data) and cross-moment vector

As $n\to+\infty$, suitable moment assumptions on $\pi$ give the following almost sure limits by the law of large numbers:

For appropriately chosen $\lambda_n \rightarrow 0$, the estimator can converge as $n \rightarrow +\infty$ to the population parameter

Problem (13.14) is equivalent to the regularized resolution of inverse problems (9.9), with $X/\sqrt n$ in place of $\Phi$ and $y/\sqrt n$ as the observation vector.

Random design introduces sampling error into the operator as well: $\hat C$ estimates $C$. The typical scale is $1/\sqrt{n}$, in the sense that

assuming $\mathbb{E}(\mathbf{y}^4)<+\infty$ and $\mathbb{E}(|\!| \mathbf{x} |\!|^4)<+\infty$, which give finite second moments for $\mathbf{x}\mathbf{x}^*$ and $\mathbf{x}\mathbf{y}$. Using a linear predictor does not require a correctly specified model $\mathbf{y}=\langle \mathbf{x},\,\beta\rangle + w$ with noise $w$ independent of $\mathbf{x}$. The aim is to estimate the best linear predictor $\bar\beta$.

Extensions of Theorems 9.4, 11.9, and 11.11 can account for covariance-estimation error and yield prediction bounds of the form

and estimation rates of the form

under suitable source conditions involving $C$ and $u$.

Since the noise level is roughly $n^{-\frac{1}{2}}$, the ideal cases are when $\kappa=\kappa'=1$, which is the so-called linear rate regime.

Support-recovery guarantees can also be derived by extending Theorem 11.15. We now focus on quadratic regularization, a standard regression method. Sparse penalties such as $\ell^1$ are useful when feature selection or a sparse model is appropriate; their statistical benefit depends on the data and assumptions.

In imaging inverse problems, sparse regularization expresses prior structure of the object being reconstructed.

Ridge regression (quadratic penalization).

For $J=|\!| \cdot |\!|^2/2$, the estimator (13.14) is obtained in closed form as

This is often called ridge regression in the literature.

The Woodbury identity gives the equivalent expression

When $n\gg p$, formula (13.16) requires solving the smaller, $p\times p$ system. When $p\gg n$, formula (13.17) uses an $n\times n$ system and is preferable; it also extends to infinite-dimensional feature spaces.

Convergence to the minimum-norm population solution $\bar\beta=C^+u$ requires control of both the regularization bias and sampling error. The condition $\lambda_n\to0$ alone is insufficient when small eigenvalues amplify noise. Here is one elementary bound that also applies in a Hilbert feature space.

In finite dimensions, a source representation exists for every $\bar\beta\in\mathop{\mathrm{Im}}(C)$, but its radius $\rho$ depends on the spectrum and on the target. In infinite dimensions, it is an additional regularity assumption. The bound displays no explicit feature dimension, although $\kappa$, $B$, and $\rho$ may depend on it. Standard ridge regularization has bias saturation beyond source exponent $\gamma=1$; stronger source assumptions alone do not justify extending this estimate to $\gamma=2$.

13.4 Supervised Learning: Classification¶

For classification, the labels are discrete: $y_i \in \mathcal{Y}= \{1,\ldots,k\}$.

We study two standard methods: nearest neighbors and logistic classification.

Both methods provide useful baselines for assessing more elaborate models.

Nearest-neighbor methods also apply to regression.

13.4.1 Nearest-Neighbor Classification¶

The $R$-nearest-neighbor method ($R$-NN) classifies an input using the labels of its $R$ nearest training observations. Increasing $R$ reduces sensitivity to individual noisy labels, but excessive smoothing can obscure class boundaries. For consistency, a standard asymptotic choice has $R=R_n\to\infty$ and $R_n/n\to0$.

Nearest-neighbor classification is a useful baseline, especially when the feature dimension $p$ is small.

The prediction $\hat f(x)\in\mathcal{Y}$ is the most frequent label among those $R$ observations, with a fixed rule for breaking ties. The method is nonparametric: its representation depends on the training sample rather than on a fixed-dimensional parameter vector.

First compute the Euclidean distance from the query $x$ to each training observation $x_i$.

Sorting the distances generates an indexing $\sigma$ (a permutation of $\{1,\ldots,n\}$) such that

For a given $R$, one can compute the “local” histogram of classes around $x$

The predicted class at $x$ is an index attaining the largest histogram value:

Choose $R$ by validation or cross-validation. For classification, the validation risk (13.13) typically uses the 0–1 loss, giving the fraction of misclassified observations

where $\delta(0)=0$ and $\delta(s)=1$ if $s \neq 0$.

The method also applies to features in a general metric space in place of $\mathbb{R}^p$.

Fast nearest-neighbor search can avoid explicitly sorting every distance.

Figure 13.14 shows the predicted class regions $\left\{ x \;;\; \hat f(x)=\ell \right\}$, for $\ell=1,\ldots,k$, in a two-dimensional projection of the Iris dataset.

Increasing $R$ leads to smoother class boundaries.

13.4.2 Binary Logistic Classification¶

Logistic classification is a widely used baseline for binary and multiclass prediction. It outputs a probability for each class, providing a richer prediction than a class label alone. The accuracy of these probabilities must still be assessed on held-out data.

The standard name for this model is logistic regression, even when the prediction task is classification.

Support vector machines provide another common approach. Their hinge loss is nonsmooth, and their scores are not probabilities without an additional calibration step. The choice between these methods depends on the task and evaluation criterion.

For the binary formulas below, use labels $y_i \in \mathcal{Y}= \{-1,1\}$. In logistic classification, the prediction $f(\cdot,\beta) \in [0,1]$ is a probability for the positive class rather than a class label. We retain the notation $f$ for this probability predictor.

Approximate risk minimization.

A linear score $\langle x,\,\beta\rangle$ defines the classifier $\mathop{\mathrm{sign}}(\langle x,\,\beta\rangle)$, with either class assigned when the score is zero. The empirical 0–1 objective counts classification errors, leading to the minimization problem

where $\ell_0=1_{[0,+\infty)}$ counts zero margins as errors. Away from ties, misclassification corresponds to $\langle x_i,\,\beta\rangle$ and $y_i$ having different signs, so that in this case $\ell_0(-y_i \langle x_i,\,\beta\rangle)=1$ (and 0 otherwise for correct classification).

The loss $\ell_0$ is nonconvex, and solving (13.20) globally is NP-hard in general. Convex upper bounds on $\ell_0$ provide tractable surrogate objectives.

Two common surrogates are

These are, respectively, the nonsmooth hinge loss used in support vector machines and the smooth logistic loss. The $1/\log(2)$ is just a constant which makes $\ell_0 \leqslant\ell$.

AdaBoost uses the exponential loss $\ell(u)=e^u$.

Least squares uses $\ell(u)=(1+u)^2$. Its growth for large negative margins makes it a poor pointwise approximation of $\ell_0$, although it can still give useful classifiers.

Unregularized hinge-loss minimization is equivalent to a linear program with nonnegative slack variables:

Logistic loss probabilistic interpretation.

Logistic regression applies a sigmoid to the linear score from Section 13.3.1. The resulting probability predictor is nonlinear. The predicted probability of class +1 is

which is often called the “logit” model.

Despite its simplicity, a linear classifier can be effective when the feature dimension $p$ is large.

With $s=\langle x,\,\beta\rangle$, the predicted probability of class -1 is $1-f(x,\beta)=\theta(-s)$. The identity $\theta(-s)=1-\theta(s)$ makes the two probabilities sum to one. It does not require balanced class frequencies; an intercept can account for unequal prior frequencies.

For a nonzero parameter, $\beta/|\!| \beta |\!|$ determines the normal direction of the separating hyperplane, while $1/|\!| \beta |\!|$ controls the width of the probability transition. As $|\!| \beta |\!| \rightarrow +\infty$ along a fixed direction, the transition approaches a sharp decision boundary.

The predictor $f(x,\beta)$ is also a single-layer perceptron with a logistic (sigmoid) activation; see Chapter 17.

For conditionally independent labels $y_i\in\{-1,+1\}$, write $\bar y_i=(y_i+1)/2\in\{0,1\}$, $s_i=\langle x_i,\,\beta\rangle$, and $p_i=\theta(s_i)$. The likelihood is

Thus the negative log-likelihood is

After division by the sample size, this gives the empirical-risk objective (13.11) with logistic loss:

where the logistic loss reads

Problem (13.22) is a smooth convex minimization. If $X$ has full column rank, $E$ is strictly convex, so it has at most one minimizer. A finite minimizer need not exist for separable data. A positive quadratic penalty guarantees existence and uniqueness.

Figure 13.16 compares the 0–1, logistic, and hinge losses.

Gradient descent method.

Write the objective as

Its gradient is

where $\odot$ is the pointwise multiplication operator, i.e. .* in Matlab.

Once $\beta^{(\ell=0)} \in \mathbb{R}^{p}$ is initialized (for instance at $0_{p}$), one step of gradient descent (14.13) reads

Figure 13.17 illustrates the method on samples from a Gaussian mixture, whose overlap is controlled by an offset $\omega$.

Overlaying the data on the probability map highlights the separating hyperplane $\left\{ x \;;\; \langle \beta,\,x\rangle=0 \right\}$.

13.4.3 Multiclass Logistic Classification¶

For multiclass logistic classification, introduce one weight vector for each of the $k$ classes. Collect the vectors $\beta= ( \beta_\ell )_{\ell=1}^k$ as columns of $\beta\in \mathbb{R}^{p \times k}$.

For an input $x \in \mathbb{R}^p$, convert the class scores into probabilities using the softmax model:

The vector $f(x,\beta)\in[0,1]^k$ gives the probabilities that $x$ belongs to each class, and satisfies $\sum_\ell f(x,\beta)_\ell=1$.

Estimate $\beta$ by maximizing the log-likelihood over the training sample:

where $y_i \in \mathcal{Y}= \{1,\ldots,k\}$ is the label of observation $x_i$.

This is conveniently rewritten as

where $D \in \{0,1\}^{n \times k}$ is the one-hot class matrix

and LSE is the log-sum-exp operator

For $k=2$, identify the first class with label +1 and the second with label -1. The model then agrees with binary logistic classification in Section 13.4.2, with parameter $\beta_1-\beta_2$. A single weight vector therefore suffices in that case.

Direct evaluation of log-sum-exp can overflow when $S_{i,\ell}$ is large. A stable evaluation subtracts the largest entry in each row before exponentiation, using the identity

for a rowwise constant $a$. This shift is especially useful for well-separated classes, whose parameter vectors $\beta_\ell$ may grow large.

The gradient of the sum of the rowwise log-sum-exp values is the softmax operator

As with log-sum-exp, subtract the largest entry in each row before evaluating the exponentials.

Once the matrix $D$ is computed, the gradient of $\mathcal{E}$ is computed as

The objective $\mathcal{E}$ can then be minimized by gradient descent.

To illustrate the method, we use a dataset of $n$ images of size $p = 8 \times 8$, representing digits from 0 to 9 (so there are $k=10$ classes).

Figure 13.18 displays a few representative examples as well as 2-D and 3-D PCA projections.

Figure 13.19 visualizes the fitted class probabilities $f(x,\hat\beta)$ by assigning colors on a regular image grid.

13.5 Kernel Methods¶

Linear models may fail to capture the geometry needed for a regression or classification task. Moreover, inputs such as strings and graphs need not belong to a vector space.

Kernel methods introduce a feature map into a richer space and apply linear prediction there. This produces nonlinear regression and classification rules while retaining linear-system solvers and convex optimization. The kernel trick expresses the computations through similarities between the $n$ observations, without constructing the feature vectors explicitly. Such methods are nonparametric when the number of fitted coefficients grows with $n$, allowing the model to become more expressive as data accumulate.

Algorithms that use the observations $x_i$ only through inner products can often be kernelized. Examples include linear regression, nearest-neighbor methods, SVMs, logistic classification, and PCA. We introduce the construction through ridge regression and then extend it to other losses and distance-based methods.

13.5.1 Feature Map and Kernels¶

Let $\varphi: \mathcal{X}\rightarrow \mathcal{H}$ be a feature map into a real Hilbert space $\mathcal{H}$, which may be infinite dimensional.

The input space $\mathcal{X}$ need not be a vector space. When $\mathcal{X}=\mathbb{R}^p$, the choice $\varphi(x)=x$ recovers ordinary linear methods.

For regression, we consider predictions which are linear functions over this lifted space, i.e. of the form $x \mapsto \langle \varphi(x),\,\beta\rangle_\mathcal{H}$ for some weight vector $\beta\in \mathcal{H}$. For binary classification, one uses $x \mapsto \mathop{\mathrm{sign}}(\langle \varphi(x),\,\beta\rangle_\mathcal{H})$.

For scalar inputs ($p=1$), the feature map $\varphi(x) = (1,x,x^2,\ldots,x^k)\in\mathbb{R}^{k+1}$ gives polynomial regression. Multivariate monomials extend this construction to higher-dimensional inputs.

Let $\Phi = ( \varphi(x_i)^* )_{i=1}^n$ be the feature matrix, whose rows are $\varphi(x_i)$. For $\mathcal{H}=\mathbb{R}^{\bar p}$, it is an ordinary matrix $\Phi \in \mathbb{R}^{n \times \bar p}$. In an infinite-dimensional feature space, the same notation represents a linear operator.

For $\lambda>0$, consider the regularized empirical-risk problem

Start with regression under the squared loss $\ell(y,z)=\frac{1}{2}(y-z)^2$:

which is the ridge regression problem studied in Section 13.3.1.

Setting the gradient to zero gives

Direct computation may be impractical or impossible when $\mathcal{H}$ is infinite dimensional. The following Woodbury identity converts the problem into a finite system.

This is the feature-space counterpart of (13.17), with $\Phi$ replacing $X$ and $\lambda$ replacing $n\lambda_n$.

The identity gives the alternative representation

where the empirical kernel matrix is

The coefficients $c^\star \in \mathbb{R}^n$ depend on the inputs only through $\kappa$, without explicit computation of $\varphi(x)$. If $\kappa$ is inexpensive to evaluate, forming $K$ takes $O(n^2)$ kernel evaluations. A dense direct solve costs $O(n^3)$ operations; an iterative solve costs $O(n^2)$ per matrix-vector product, with the iteration count determined by conditioning and accuracy.

The fitted predictor at a new input $x \in \mathcal{X}$ also uses only kernel evaluations: