1 Notes on this Version and Current Status↩︎

**This is a preprint of a book in preparation by the authors.**

**We anticipate the focus and content of this manuscript not to change except for the correction of typos and other errors and other minor changes such as including additional homework problems. It can already by used as textbook for a graduate course in Mathematics of Data Science; it has been used as such by the authors at their home institutions.**

We welcome suggestions and comments, and would like to learn about any possible errors and typos.

**Please contact the authors atbandeira@math.ethz.ch,
strohmer@math.ucdavis.edu, oramits@math.princeton.edu. In fact, this version was improved using many suggestions sent from, and typos found by, many students and readers. Thank you!**

Thank you,

Afonso, Thomas, and Amit.

2 Introduction↩︎

2.1 Origins↩︎

The idea of this book emerged from the classes on “mathematics of data science” that the authors have been teaching at their home institutions Princeton University (AS), UC Davis (TS), and MIT, NYU, and ETH (ASB)¹. The lecture notes eventually became the basis to this book. While there are specialized books and many advanced research papers on each topic covered in this book, we believe it fills in the gap of being a single source that provides an accessible entry point at an advanced undergraduate, or graduate, level to mathematics of data science. The book can serve readers as a segue to those more advanced and specialized topis. In that light, one of the main purposes of this book is to be used as a textbook at both the (advanced) undergraduate and graduate levels. This book also differs from most texts in the way that it interweaves mathematical theory with some of its applications in data science.

2.2 The data science revolution↩︎

The data science revolution is one of the most profound shifts in modern history, reshaping industries, economies, and societies with unprecedented speed and scale. At its core, the data science revolution is driven by the exponential growth in data, advances in computing power, and the development of sophisticated analytical techniques.

Data science and its cousin artificial intelligence have become a cornerstone of modern decision-making and innovation, penetrating almost all aspects of our life. In healthcare, for example, data science is driving personalized medicine and predictive analytics. By analyzing patient data, healthcare providers can tailor treatments to individual needs, predict potential health issues before they arise, and improve patient outcomes. In finance, data science is revolutionizing risk management, fraud detection, and investment strategies. Algorithms analyze market trends, assess creditworthiness, and identify unusual transactions to mitigate risks and enhance decision-making.

Retailers leverage data science to understand consumer behavior, optimize supply chains, and personalize marketing efforts. The transportation industry benefits from data science through advancements in autonomous vehicles, route optimization, and predictive maintenance. Companies use data to enhance logistics, reduce costs, and improve safety. The energy sector uses data science to optimize resource allocation, predict equipment failures, and enhance energy efficiency. Data-driven insights help in managing renewable energy sources and reducing environmental impact.

From the statistical techniques that form the basis of data analysis to the complex algorithms that power machine learning and artificial intelligence, mathematics is at the heart of data science. Understanding the connection between data science and mathematics is crucial for appreciating the depth and breadth of this revolution.

2.3 Mathematics as a foundation of data science↩︎

Data science is an inherently interdisciplinary field that merges techniques from mathematics, statistics, computer science, signal processing, and domain-specific knowledge to extract insights from data. However, at its core, data science is fundamentally about understanding and leveraging patterns in data. This understanding is built on a complex and nuanced foundation of mathematical principles that enable the development of models and algorithms capable of making sense of complex datasets. This book aims to unravel many of these mathematical principles and provide a thorough grounding in the essentials of data science from a rigorous, mathematical perspective.

Mathematics provides the language and framework for formalizing the processes involved in data analysis, from data collection and preprocessing to modeling and evaluation. It offers the tools necessary for designing effective algorithms, making precise predictions, visualizing complex datasets, quantifying uncertainty, and interpreting the results of data-driven experiments. As the field of data science keeps evolving at a rapid speed, new algorithms and techniques are constantly being developed. A strong mathematical foundation enables practitioners to quickly understand and adopt these innovations because they can grasp the underlying principles rather than just following black-box implementations. Moreover, with a solid grasp of mathematics, learning new methods becomes easier.

Also, when algorithms fail to produce correct results in real world applications, we would like to know why they failed. Is it because of mistakes in the experimental setup, did we gather the wrong kind of data, do we have insufficient data, is it because of corrupted measurements, calibration errors, or incorrect modeling assumptions, or perhaps due to fundamentally unrealistic expectations, or is it due to a deficiency of the algorithm itself? If it is the latter, can it be fixed by a better initialization, a more careful tuning of the parameters, or by choosing a different algorithm? Or is a more fundamental modification required, such as developing a different model, including additional prior information, taking more measurements, or a better compensation of calibration errors? How reliable are our algorithms to small changes in the input data and can we explain why an algorithm arrives at a certain prediction?

As we are dealing with larger datasets and we aim for faster and faster throughput, it becomes increasingly important to address the aforementioned challenges in a systematic and principled manner. Thus, a rigorous and thorough study of computational algorithms both from a theoretical and numerical viewpoint is not a luxury, but emerges as an imperative ingredient towards effective data-driven discovery.

While a thorough mathematical framework can never completely avoid some guessing and trial-and-error, it can dramatically reduce the time and resources we spend on trial-and-error. Good mathematical theory can help in predicting if a problem is solvable at all in feasible time. Moreover, fundamental mathematical analysis often leads to new discoveries such as faster or more robust algorithms, and it may reveal hitherto unknown connections to other, already well studied problems.

At the heart of data science lies a rich tapestry of mathematical concepts. Key areas include statistics and probability, linear algebra, optimization, numerical algorithms, analysis, and graph theory, but this list is by no means exhaustive. This book emphasizes the geometric viewpoint of data sets. Specifically, it is a very powerful way to think about a dataset as a point cloud in a high dimensional space for a wide range of analysis tasks such as clustering, classification, dimension reduction, prediction, denoising, outlier detection, and more. In this geometric viewpoint, linear models of data are interpreted as linear subspaces in a Euclidean space, while non-linear models can be viewed as manifolds embedded in the ambient space. Equally powerful is the representation of data points and their interconnections and similarities as a graph and learning using graph algorithms and analysis. In particular, the graph representation provides a powerful way to extend classical methods for function approximation and interpolation to the setting of complex and high dimensional setting of modern data sets. The computational machinery provided by classical harmonic analysis on Euclidean spaces such as Fourier analysis can be extended in this way to more complex geometries and networks.

As data science continues to evolve, new mathematical techniques and approaches will emerge, expanding the possibilities for analyzing and interpreting data. This book provides a solid foundation upon which readers can build their expertise, adapt to new developments, and contribute to the advancement of the field.

Data science also opens a new playground for mathematics, motivating new questions and new areas of mathematical inquiry. For this reason, another goal of this book is to motivate mathematicians (and, in particular, mathematics students) to engage with these topics and to create new mathematics motivated by questions arising from data and computation.

2.4 Whom is this book for?↩︎

Many of the concepts, algorithms, and methods presented in this book do not only belong in the toolbox of a well trained data scientist, but are also found in the standard repertoire of researchers and practitioners of machine learning and AI.

The purpose of this book is twofold. On the one hand. the book is designed to guide readers through many of the mathematical foundations of data science in a structured and accessible manner. It assumes a basic familiarity with mathematical concepts and aims to build upon this foundation to explore more advanced topics relevant to data science. On the other hand, the selection of data science topics for this book is also guided by our desire to demonstrate how data science leads to deep and exciting mathematical problems, and how (sometimes seemingly unrelated) mathematical concepts can provide key insight into data science tasks. This mutual enrichment is akin to how physics and mathematics have had a vibrant cross-fertilization over many centuries.

While the tools and methods of data science and machine learning that are currently en vogue tend to change quickly, a solid grasp of the underlying mathematical principles provides an enduring foundation that allows one to adapt to, and even shape, the innovations of future data science. We sincerely hope that this book can make a contribution to this endeavor.

What prerequisites should a reader have? An undergraduate level knowledge of linear algebra and probability is assumed.

2.5 Book organization and some omissions↩︎

What is covered in this book? This book contains the material that one might want to teach in either a year-long course or a one-semester course on the mathematical foundations of data science.

The book is separated into sixteen chapters, and an effort was made to intertwine mathematical techniques and motivating applications. It starts with a description of the curious phenomena of high dimensional geometry, taking the opportunity to introduce and recall some probability theory in Chapter 3. It proceeds to discussing linear dimension reduction techniques (such as PCA and truncated SVD) and doing a recap of several linear algebra concepts in Chapter 4. Linear regression and least squares are the central topics in Chapter 5.

We move on to graph theory (including spectral graph theory, and spectral clustering in Chapter 6) and non-linear dimension reduction (including diffusion maps in Chapter 7), continuing to highlight how central the spectrum of matrices formed from data are to understanding the data geometry. Chapter 8 focuses on dimension reduction via randomized methods, including a brief introduction to randomized linear algebra. In Chapter 9 we dive into the tools of optimization that play a key role in data science, such as convex optimization and gradient descent. We explore several classical classification methods in Chapter 10 and present a gentle introduction to Deep Learning in Chapter 11.

In Chapter 12 we analyze the large sample limit of diffusion maps and other non-linear dimensionality reduction methods. The concept and power of convex relaxations is illustrated via the community detection problem in Chapter 13. Chapters 14 and 15 are devoted to probability, the first to Gaussian analysis and concentration of measure, and the second to matrix concentration inequalities. Finally, Chapters 16 and 17 deal with sparse data models, the first one is dedicated to Compressive Sensing, while the second one covers low-rank matrix recovery.

What is not covered in this book? This book does not cover everything related to mathematics of data science. Indeed, there are several topics not included to keep the book at a manageable length. Perhaps the most notable omissions are relevant topics from statistics such as statistical inference and parameter estimation (e.g., Bayesian inference, etc.). There are already excellent textbooks dedicated to these topics, which are also being typically taught in separate courses. Such courses in statistics are of course recommended to those interested in data science, and can be taken either before, after, or in parallel to a course based on this book. The book also does not cover more specialized (but definitely very important) topics like data privacy, which is deeply connected to the mathematics of data science. The book is also at fault to not discuss optimal transport which has gained a lot of attention in mathematical data science. And we do not dive too deep into learning theory (such as VC dimension), as there are already plenty of resources on this topic. Notwithstanding these deliberate omissions, we believe that this book covers most of the foundational mathematical pillars of data science.

2.6 Some related resources↩︎

At https://video.ethz.ch/lectures/d-math/2025/autumn/401-4944-20L you can access videos of lectures given by one of the authors on parts of this book (or see https://people.math.ethz.ch/~abandeira/videos.html for more videos on related topics). For the readers that find themselves looking for more problems or exercises: you might enjoy this note [@JigSawStatCompGaps2025] which presents many exciting new results in theoretical computer science and high dimensional statistics via mathematical exercises. For open problems related to many topics of this book, you can visit https://randomstrasse101.math.ethz.ch/ or [@RandomstrasseProblems2024] (or [@Afonso_10L42P] for a less updated list).

Acknowledgment↩︎

We want to thank our students, colleagues, and readers of earlier drafts of this book, who provided valuable feedback, alerted us to errors, and pointed out typos. Please keep doing so! We are grateful to Junda (Albie) Sheng and Stefan Broecker for generating some of the figures and to Yuan Ni for providing code for some numerical simulations. We are very thankful for the invaluable help that teaching assistants provided while we teach our courses on this topic, including the designing of countless problems and exercises on these topics. A particular thank you to Pedro Abdalla, Kevin Lucca, Anastasia Kireeva, Antoine Maillard, Petar Nizíc-Nikolac, and Almut Rödder, who have greatly helped with the teaching of classes on Mathematics of Data Science at ETH; many exercises in the book originated as homework problems in these classes [@ExercisesInMDS2025]. Amit Singer acknowledges support from the Simons Foundation Math+X Investigator Award, AFOSR FA9550-20-1-0266, AFOSR FA9550-23-1-0249, NSF DMS 2009753, NSF DMS 2510039, and NIH/NIGMS R01GM136780-01. Thomas Strohmer acknowledges support from NSF DMS-2208356, NIH R01HL16351, P41EB032840, and DE-SC0023490.

3 Curses, Blessings, and Surprises in High Dimensions↩︎

Most of the material in this book will naturally “live” in the context of high dimensions, either motivated by the analysis of high dimensional data or by associating with networks a space whose dimension is the number of network nodes, to name a few. This chapter is about High Dimensions. It discusses the curse of dimensionality, but also many of its blessings. The first is caused by the exponential increase in volume associated with adding extra dimensions to Euclidean space. The latter is a manifestation of an intriguing phenomenon called the concentration of measure. This concentration phenomenon will give rise to many surprising facts about high dimensional geometry that we will discuss. Since several of the results discussed in this chapter require basic tools from probability, we will also review some fundamental probabilistic concepts.

3.0.0.1 The curse of dimensionality

The curse of dimensionality refers to the fact that many algorithmic approaches to problems in \(\mathbb{R}^d\) become exponentially more difficult as the dimension \(d\) grows. The expression “curse of dimensionality” was coined by Richard Bellman to describe the problem caused by the exponential increase in volume associated with adding extra dimensions to Euclidean space [@bellman1957].

For instance, if we want to sample the unit interval such that the distance between adjacent points is at most 0.01, then 100 evenly-spaced sample points suffice; an equivalent sampling of a five-dimensional unit hypercube with a grid with a spacing of 0.01 between adjacent points would require \(10^{10}\) sample points. Thus, a modest increase in dimensions results in a dramatic increase in required data points to cover the space at the same density.

Intimately connected to the curse of dimensionality is the problem of overfitting and underfitting. Here, overfitting refers to the issue that an algorithm may show good performance on the training data, but poor generalization to other data. Underfitting in turn, corresponds to poor performance on the training data (and poor generalization to other data). This problem manifests itself in many machine learning algorithms.

We will discuss a toy example from image classification in more detail to illustrate the underlying issues. Assume we want to classify images into two groups, cars and bicycles, say. From the vast number of images depicting cars or bicycles, we are only able to obtain a small number of training images, say five images of cars and five images of bicycles. We want to train a simple linear classifier based on these ten labeled training images to correctly classify the remaining unlabeled car/bicycle images. We start with a simple feature, e.g.the amount of red pixels in each image. However, this is unlikely to give a linear separation of the training data. We add more features and eventually the training images become linearly separable. This might suggest that increasing the number of features until perfect classification of the training data is achieved, is a sound strategy. However, as we linearly increase the dimension of the feature space, the density of our training data decreases exponentially with the feature dimension.

In other words, to maintain a comparable density of our training data, we would need to increase the size of the dataset exponentially – the curse of dimensionality. Thus, we risk producing a model that could be very good at predicting the target class on the training set, but it may fail miserably when faced with new data. This means that our model does not generalize from the training data to the test data, leading to large generalization error. The generalization error is an important concept in data science and machine learning, since it measures how accurately an algorithm can predict outcomes for new, unseen data. As such, it indicates the data science model’s ability to perform well beyond the training data, and is closely related to the concept of overfitting.

3.1 Geometry of spheres and cubes in high dimension↩︎

When we peel an orange, then after having removed the rind we are still left with the majority of the orange. Suppose now we peel a \(d\)-dimensional orange for large \(d\), then after removing the orange peel we would be left with essentially nothing. The reason for this – from a healthy nutrition viewpoint discouraging – fact is that for a \(d\)-dimensional unit ball almost all of its volume is concentrated near the boundary sphere. This is just one of many surprising phenomena in high dimensions. Many of these surprises are actually a manifestation of some form of concentration of measure that we will analyze in more detail in the next section (and then these surprises are no longer so startling).

When introducing data analysis concepts, we typically use only a few dimensions in order to facilitate visualization. However, our intuition about space, which is naturally based on two and three dimensions, can often be misleading in high dimensions. Many properties of even very basic objects become counterintuitive in higher dimensions. Understanding these paradoxical properties is essential in data analysis as it allows us to avoid pitfalls in the design of algorithms and statistical methods for high-dimensional data. It is therefore instructive to analyze the shape and properties of some basic geometric forms that we understand very well in dimensions two and three, in high dimensions.

To that end, we will look at some of the properties of the sphere and the cube as the dimension increases. The \(d\)-dimensional hyperball of radius \(R\) is defined by \[B^d(R) =\{x \in \mathbb{R}^d : x_1^2 + \dots +x_d^2 \le R^2\},\] the \(d\)-dimensional hypersphere (or \(d\)-sphere) of radius \(R\) is given by \[S^{d-1}(R) =\{x \in \mathbb{R}^d : x_1^2 + \dots +x_d^2 = R^2\},\] and the \(d\)-dimensional hypercube with side length \(2R\) is the subset of \(\mathbb{R}^d\) defined as the \(d\)-fold product of intervals \([-R, R]\): \[C^d(R) = \underbrace{ [-R,R] \times \dots \times [-R,R]}_{\text{d times}}.\] If there is no danger of confusion, we may write \(B^d\) for the unit ball \(B^d(1)\), \(S^{d-1}\) for the unit sphere \(S^{d-1}(1)\), and \(C^d\) for the unit hypercube \(C^d(\frac{1}{2})\).

3.1.1 High dimensional geometry and high dimensional probability↩︎

In the rest of this section we will establish a variety of (at first) surprising properties of high dimensional geometry. To illustrate the connections between geometry and probability we will derive some of the results first using analytical tools (computing volumes and integrals) and then using tools from Probability, namely concentration inequalities (which the reader might agree gives more elegant arguments). The goal here is not to get the best bounds, but to illustrate these properties and showcase the duality between geometry and probability.

3.1.1.1 Volume of the hyperball

Recall that for even \(d\), we have \(\Gamma(d/2) = (d/2 - 1)!\). To understand the asymptotic behavior \(\mathop{\mathrm{Vol}}(B^d(R))\) we can use the celebrated Stirling’s Formula, \[\begin{align} \Gamma(n) \approx \sqrt{\frac{2\pi}{n}}\left(\frac{n}{e}\right)^n \end{align}\] we obtain an approximation for the volume of the unit \(d\)-ball for large \(d\) \[\begin{align} \label{approxvolume} \mathop{\mathrm{Vol}}(B^d) \approx \frac{1}{\sqrt{d\pi}} \Big( \frac{2\pi e}{d} \Big)^{\frac{d}{2}}, \end{align}\tag{3}\] where \(\approx\) means that the quotient goes to \(1\) as the relevant parameter, in this case \(d\), goes to \(\infty\).

Since the denominator in the parenthesis of equation 3 goes to infinity, the volume of the unit \(d\)-sphere goes rapidly to \(0\) as the dimension \(d\) increases to infinity, see also Figure 1. Notice that this would still be the case for any fixed radius \(R\).

Thus, unit spheres in high dimensions have almost no volume—compare this to the unit cube, which has volume \(1\) in any dimension. For \(B^d(R)\) to have volume equal to 1, its radius \(R\) must be approximately (asymptotically) equal to \(\sqrt{\frac{d}{2\pi e}}\).

3.1.1.2 Concentration of the volume of a ball near its equator

If we take an orange and cut it into slices, then the slices near the center are larger since the sphere is wider there. This effect increases dramatically (exponentially with the dimension) with increasing dimension. Assume we want to cut off a slab around the “equator²” of the \(d\)-unit ball such that 99% of its volume is contained inside the slab. In two dimensions the width of the slab has to be almost 2, so that 99% of the volume are captured by the slab. But as the dimension increases the width of the slab gets rapidly smaller. Indeed, in high dimensions only a very thin slab is required, since nearly all the volume of the unit ball lies a very small distance away from the equator. The following theorem makes the considerations above precise.

3.1.2 Concentration of the volume of a ball on shells↩︎

We consider two concentric balls \(B^d(1)\) and \(B^d(1-\varepsilon)\). Using equation ?? , the ratio of their volumes is \[\begin{align} \frac{\mathop{\mathrm{Vol}}\big(B^d(1-\varepsilon)\big)}{\mathop{\mathrm{Vol}}\big(B^d(1)\big)} = (1-\varepsilon)^d. \end{align}\] Clearly, for every \(\varepsilon\) this ratio tends to zero as \(d \to \infty\). This implies that the spherical shell given by the region between \(B^d(1)\) and \(B^d(1-\varepsilon)\) will contain most of the volume of \(B^d(1)\) for large enough \(d\) even if \(\varepsilon\) is very small. How quickly does the volume concentrate at the surface of \(B^d(1)\)? We choose \(\varepsilon\) as a function of \(d\), e.g. \(\varepsilon= \frac{t}{d}\), then \[\frac{\mathop{\mathrm{Vol}}(B^d(1-\varepsilon))}{\mathop{\mathrm{Vol}}(B^d(1))} = \Big(1-\frac{t}{d}\Big)^d \to e^{-t}.\] Thus, almost all the volume of \(B^d(R)\) is contained in an annulus of width \(R/d\).

Therefore, if we peel a \(d\)-dimensional orange—and even if we peel it very carefully so that we remove only a very thin layer of its peel—we will have removed most of the orange and are left with almost nothing.

3.1.2.1 Geometry of the hypercube

We have seen that most of the volume of the hypersphere is concentrated near its surface. A similar result also holds for the hypercube, and it is indeed a tendency for high-dimensional geometric objects. Yet, the hypercube exhibits an even more interesting volume concentration behavior, which we will establish below.

We start with a basic observation.

The above proposition, although mathematically trivial, hints already at a somewhat counterintuitive behavior of the cube in high dimensions. Its corners seem to get “stretched out” more and more, while the rest of the cube must “shrink” to keep the volume constant. This property becomes even more striking when we compare the cube with the sphere as the dimension increases.

In two dimensions (Figure 2), the unit square is completely contained in the unit sphere. The distance from the center to a vertex (radius of the circumscribed sphere) is \(\frac{\sqrt{2}}{2}\) and the apothem (radius of the inscribed sphere) is \(\frac{1}{2}\). In four dimensions (Figure 3), the distance from the center to a vertex is \(1\), so the vertices of the cube touch the surface of the sphere. However, the apothem is still \(\frac{1}{2}\). The result, when projected in two dimensions no longer appears convex, however all hypercubes are convex. This is part of the strangeness of higher dimensions - hypercubes are both convex and “pointy.” In dimensions greater than \(4\) the distance from the center to a vertex is \(\frac{\sqrt{d}}{2} > 1\), and thus the vertices of the hypercube extend far outside the sphere, cf. Figure 4.

The considerations above suggest the following observation:

3.2 Basic concepts from probability↩︎

We briefly review some fundamental concepts from probability theory, which are helpful or necessary to understand the blessings of dimensionality and some of the surprises encountered in high dimensions. More advanced probabilistic concepts will be presented in Chapter 14. We assume that the reader is familiar with elementary probability as is covered in introductory probability courses (see, for example [@durrett2019probability; @ross2014introduction]).

The two most basic concepts in probability associated with a random variable \(X\) are expectation (or mean) and variance, denoted by \[\mathbb{E}[X] \qquad \text{and} \qquad \mathop{\mathrm{Var}}(X) := \mathbb{E}[ (X- \mathbb{E}[X])^2 ],\] respectively.

A useful identity is the Law of Total Variance. This fundamental theorem in probability theory decomposes the variance of a random variable into two distinct components: the variance that can be explained by another variable and the variance that remains as “noise.” Let \(X\) and \(Y\) be random variables on the same probability space, and assume that the variance of \(X\) is finite. The Law of Total Variance states \[\begin{align} \label{eq:totalvariance} \mathop{\mathrm{Var}}(X) = \mathbb{E}[\mathop{\mathrm{Var}}(X \mid Y)] + \mathop{\mathrm{Var}}(\mathbb{E}[X \mid Y]). \end{align}\tag{6}\] The two terms have a natural interpretation:

• \(\mathbb{E}[\mathop{\mathrm{Var}}(X \mid Y)]\) is the unexplained variance — the variability in \(X\) that persists even after knowing \(Y\).

• \(\mathop{\mathrm{Var}}(\mathbb{E}[X \mid Y])\) is the is the variance of the conditional mean — the variability in \(X\) explained by \(Y\).

An important tool to describe probability distributions is the moment generating function (MGF) of \(X\), defined by \[\label{eq:MomentGeneratingFunction} M_X(t) = \mathbb{E}[e^{tX}], \quad t\in \mathbb{R},\tag{7}\] the choice of nomenclature can be easily justified by expanding \(M_X(t)\) in a series. The \(p\)-th moment of \(X\) is defined by \(\mathbb{E}[X^p]\) for \(p >0\) and the \(p\)-th absolute moment is \(\mathbb{E}[|X|^p]\).

We can introduce \(L^p\)-norms of random variables by taking the \(p\)-th root of moments, i.e., \[\|X\|_{L^p} := \big(\mathbb{E}[|X|^p]\big)^{\frac{1}{p}}, \qquad p \in [0,\infty],\] with the usual extension to \(p=\infty\) by setting \[\|X\|_{\infty} := \operatorname{ess}\sup |X|.\]

Let \((\Omega,\Sigma,\mathbb{P})\) be a probability space, where \(\Sigma\) denotes a \(\sigma\)-algebra on the sample space \(\Omega\) and \(\mathbb{P}\) is a probability measure on \((\Omega,\Sigma).\) For fixed \(p\) the vector space \(L^p(\Omega,\Sigma,\mathbb{P})\) consists of all random variables \(X\) on \(\Omega\) with finite \(L^p\)-norm, i.e., \[L^p(\Omega,\Sigma,\mathbb{P}) = \{X : \|X\|_{L^p} < \infty\}.\] The random variable can be viewed as a (measurable) map \(X:\Omega\to\mathbb{R}\) (or \(\mathbb{C}\)). We will usually not mention the underlying probability space. For example, we will often simply write \(L^p\) for \(L^p(\Omega,\Sigma,\mathbb{P})\).

The case \(p=2\) deserves special attention since \(L^2\) is a Hilbert space with inner product and norm \[\langle X, Y \rangle_{L^2} = \mathbb{E}[XY], \qquad \|X\|_{L^2} = \big(\mathbb{E}[X^2] \big)^{\frac{1}{2}},\] respectively. Note that the standard deviation \(\sigma(X):=\sqrt{\mathop{\mathrm{Var}}(X)}\) of \(X\) can be written as \[\sigma(X) = \| X - \mathbb{E}[X] \|_{L^2}.\] The covariance of the random variables \(X\) and \(Y\) is \[\label{covariance} \operatorname{Cov}(X,Y) = \mathbb{E}[ (X - \mathbb{E}[X]) (Y - \mathbb{E}[Y])] = \langle X - \mathbb{E}[X], Y - \mathbb{E}[Y] \rangle_{L^2} .\tag{8}\]

We recall a few classical inequalities for random variables. Hölder’s inequality states that for random variables \(X\) and \(Y\) on a common probability space and \(p,q \ge 1\) with \(1/p +1/q =1\), there holds \[\label{hoelder} | \mathbb{E}[XY] | \le \|X\|_{L^p} \|Y\|_{L^q}.\tag{9}\] The special case \(p=q=2\) is the Cauchy-Schwarz inequality \[\label{cauchy-schwarz} | \mathbb{E}[XY] | \le \sqrt{ \mathbb{E}[|X|^2] \mathbb{E}[|Y|^2]}.\tag{10}\]

Jensen’s inequality states that for any random variable \(X\) and a convex function \(\phi: \mathbb{R}\to \mathbb{R}\), we have \[\label{Jensen} \phi(\mathbb{E}[X]) \le \mathbb{E}[\phi(X)].\tag{11}\]

Since \(\phi(x) = x^{q/p}\) is a convex function for \(q\geq p \geq 0\), it follows immediately from Jensen’s inequality that \[\| X\|_{L^p} \le \|X\|_{L^q} \qquad \text{for 0 \le p \le q < \infty.}\]

Minkovskii’s inequality states that for any \(p \in [0,\infty]\) and any random variables \(X, Y\), we have \[\label{minkovski} \|X + Y\|_{L^p} \le \|X \|_{L^p} + \|Y\|_{L^p},\tag{12}\] which can be viewed as the triangle inequality.

The cumulative distribution function of \(X\) is defined by \[F_X (t) =\mathbb{P}(X \le t), \quad t \in \mathbb{R}.\] We have \(\mathbb{P}\{X > t\} = 1 -F_X(t),\) where the function \(t \mapsto \mathbb{P}\{ |X| \ge t \}\) is called the tail of \(X\). The following lemma establishes a close connection between expectation and tails.

Given an event \(E\) with non-zero probability,\(\mathbb{P}(\cdot |E)\) denotes conditional probability, furthermore for a random variable \(X\) we use \(\mathbb{E}[X|E]\) to denote the conditional expectation.

3.2.1 Tail bounds↩︎

Markov’s inequality is a fundamental tool to bound the tail of a random variable in terms of its expectation.

We provide two versions of the same proof, one using the language of conditional expectations.

An important consequence of Markov’s inequality is Chebyshev’s inequality.

Chebyshev’s inequality, which follows by applying Markov’s inequality to the non-negative random variable \(Y = (X - \mathbb{E}[X])^2\), is a form of concentration inequality, as it guarantees that \(X\) must be close to its mean \(\mu\) whenever the variance of \(X\) is small. Both, Markov’s and Chebyshev’s inequality are sharp, i.e., in general they cannot be improved.

Markov’s inequality only requires the existence of the first moment. We can say a bit more if in addition the random variable \(X\) has a moment generating function in a neighborhood around zero, that is, there is a constant \(b>0\) such that \(\mathbb{E}[ e^{\lambda (X-\mu)}]\) exists for all \(\lambda \in [0,b]\). In this case we can apply Markov’s inequality to the random variable \(Y = e^{\lambda (X-\mu)}\) and obtain the generic Chernoff bound \[\label{basicchernoff} \mathbb{P}\{ X -\mu \ge t \} = \mathbb{P}\{ e^{\lambda (X-\mu)} \ge e^{\lambda t} \} \le \frac{ \mathbb{E}[e^{\lambda (X-\mu)}] }{ e^{\lambda t} }.\tag{14}\] In particular, optimizing over \(\lambda\) in order to obtain the tightest bound in 14 gives \[\log \mathbb{P}\{ X -\mu \ge t \} \le - \sup_{ \lambda \in [0,b]} \{ \lambda t - \log \mathbb{E}[e^{\lambda (X-\mu)} ] \}.\]

Gaussian random variables are among the most important random variables. A Gaussian random variable \(X\) with mean \(\mu\) and standard deviation \(\sigma\) has a probability density function given by \[\label{def:gaussian} \psi(t) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left( - \frac{(t-\mu)^2}{2\sigma^2} \right).\tag{15}\] We write \(X \sim {\mathcal{N}}(\mu,\sigma^2)\). We call a Gaussian random variable \(X\) with \(\mathbb{E}[X] = 0\) and \(\mathbb{E}[X^2] = 1\) a standard Gaussian or standard normal (random variable). In this case we have the following tail bound.

3.2.2 Sub-gaussian random variables↩︎

Note that \(\sigma^2\) is not necessarily the variance of \(X\). If \(X\) satisfies the above definition, we also say that \(X\) is sub-Gaussian with parameter \(\sigma\), or \(X\) is \((\mu,\sigma)\) sub-Gaussian in case we want to emphasize \(\mu\) as well. Clearly, owing to the symmetry in the definition, \(-X\) is sub-Gaussian if and only if \(X\) is sub-Gaussian. Obviously, any Gaussian random variable with variance \(\sigma^2\) is sub-Gaussian with parameter \(\sigma\). We refer to [@vershynin2018] for other, equivalent, definitions of sub-Gaussian random variables.

Combining the moment condition in Definition 3.1 with calculations similar to those that lead us to the Gaussian tail bounds in 3.5, yields the following concentration inequality for sub-Gaussian random variables.

An important example of non-Gaussian, but sub-Gaussian random variables are Rademacher random variables. A Rademacher random variable \(\varepsilon\) takes on the values \(\pm 1\) with equal probability and is sub-Gaussian with parameter \(\sigma\). Indeed, any bounded random variable is sub-Gaussian (see Remark 3.1).

While many important random variables have a sub-Gaussian distribution, this class does not include several frequently occurring distributions with heavier tails. A classical example is the chi-squared distribution, which we will discuss at the end of this chapter.

3.2.3 Sub-exponential random variables↩︎

Relaxing slightly the condition on the moment-generating function in Definition 3.1 leads to the class of sub-exponential random variables.

Clearly, a sub-Gaussian random variable is sub-exponential (set \(\nu = \sigma\) and \(b=0\), where \(1/b\) is interpreted as \(+\infty\)). However, the converse is not true. Take for example \(X \sim {\mathcal{N}}(0,1)\) and consider the random variable \(Z=X^2\). For \(\lambda < \frac{1}{2}\) it holds that \[\label{subexp2} \mathbb{E}[ e^{\lambda(Z-1)} ] = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{\lambda (x^2-1)} e^{-x^2/2} dx = \frac{ e^{-\lambda} }{\sqrt{1-2\lambda}}.\tag{16}\] However, for \(\lambda \ge \frac{1}{2}\) the moment-generating function does not exist, which implies that \(X^2\) is not sub-Gaussian. But \(X^2\) is sub-exponential. Indeed, a brief computation shows that \[\frac{ e^{-\lambda} }{\sqrt{1-2\lambda}} \le e^{2\lambda^2} = e^{4\lambda^2/2}, \qquad \text{for all |\lambda| \le 1/4,}\] which in turn implies that \(X^2\) is sub-exponential with parameters \((\nu,b)=(2,4)\).

Following a similar procedure that yielded sub-Gaussian tail bounds produces concentration inequalities for sub-exponential random variables. However, in this case we see two different types of concentration emerging, depending on the value of \(t\). In particular, we see a sub-exponential tail behavior.

Both the sub-Gaussian property and the sub-exponential property is preserved under summation for independent random variables, and the associated parameters transform in a simple manner (see Remark 3.1).⁴

3.3 Concentration of measure, a “blessing of dimensionality”↩︎

Suppose we wish to predict the outcome of an event of interest. One natural approach would be to compute the expected value of the object. However, how can we tell how good the expected value is to the actual outcome of the event? Without further information of how well the actual outcome concentrates around its expectation, the expected value is of little use. We would like to have an estimate for the probability that the actual outcome deviates from its expectation by a certain amount. This is exactly the role that concentration inequalities play in probability and statistics.

The concentration of measure phenomenon plays a central role in the asymptotic geometry of Banach spaces and in high dimensional probability; see the work of Vitali Milman and Ledoux’s monograph [@milman1986asymptotic; @ledoux2001concentration].

The celebrated law of large numbers of classical probability theory is the most well known form of concentration of measure; it states that sums of independent random variables are, under very mild conditions, close to their expectation with a large probability. We will see various quantitative versions of such concentration inequalities throughout this course. Some deal with sums of scalar random variables, others with sums of random vectors or sums of random matrices. Such concentration inequalities are instances of what is sometimes called Blessings of dimensionality (cf. [@donoho2000high]). This expression refers to the fact that certain random fluctuations can be well controlled in high dimensions, while it would be very complicated to make such predictive statements in moderate dimensions.

3.3.1 Large deviation inequalities↩︎

Concentration and large deviations inequalities are among the most useful tools when understanding the performance of some algorithms. We start with two of the most fundamental results in probability. We refer to Sections 1.7 and 2.4 in [@durrett2019probability] for the proofs and variations.

The celebrated central limit theorem tells us that the limiting distribution of a sum of i.i.d.random variables is always Gaussian. The best known version is probably due to Lindeberg-Lévy.

The strong law of large numbers and the central limit theorem give us qualitative statements about the behavior of a sum of i.i.d.random variables. In many applications it is desirable to be able to quantify how such a sum deviates around its mean. This is where concentration inequalities come into play.

The intuitive idea is that if we have a sum of independent random variables \[X = X_1 + \cdots + X_n,\] where \(X_i\) are i.i.d.centered random variables, then while the value of \(X\) can be of order⁵ \(\mathcal{O}(n)\) it will very likely be of order \(\mathcal{O}(\sqrt{n})\) (note that this is the order of its standard deviation). The inequalities that follow are ways of very precisely controlling the probability of \(X\) being larger (or smaller) than \(\mathcal{O}(\sqrt{n})\). While we could use, for example, Chebyshev’s inequality for this, in the inequalities that follow the probabilities will be exponentially small, rather than just quadratically small, which will be crucial in many applications to come. Moreover, unlike the classical central limit theorem, the concentration inequalities below are non-asymptotic in the sense that they hold for all fixed \(n\) and not just for \(n \to \infty\) (but the larger the \(n\), the stronger the inequalities become).

The inequality implies that fluctuations larger than \(\mathcal{O}\left(\sqrt{n}\right)\) have small probability. For example, if \(a_i=a\) for all \(i\), setting \(t=a\sqrt{2n\log n}\) yields that the probability is at most \(\frac{2}{n}\).

A natural way to attempt to capture this behaviour is by noting that the variance of \(\sum_{i=1}^n r_i\) depends on \(p\) as \(\mathrm{Var}(r_i) = p\). The inequality that follows, Bernstein’s inequality, uses the variance of the summands to improve over Hoeffding’s inequality.

The way this is going to be achieved is by strengthening the proof above, more specifically in step 17 we will use the bound on the variance to get a better estimate on \(\mathbb{E}[e^{\lambda X_i}]\) essentially by realizing that if \(X_i\) is centered, \(\mathbb{E}X_i^2 = \sigma^2\), and \(|X_i|\leq a\) then, for \(k\geq 2\), \(\mathbb{E}X_i^k \leq \mathbb{E}|X_i|^k \leq a^{k-2}\mathbb{E}|X_i|^{2} \leq a^{k-2} \sigma^2= \left( \frac{\sigma^2}{a^2} \right)a^k\).

Proof.

As before, we will prove \[{\mathbb{P}}\left\{ \sum_{i=1}^n X_i \geq t \right\} \leq \exp\left(-\frac{t^2}{2n\sigma^2 + \frac{2}{3}at}\right),\] and then union bound with the same result for \(- \sum_{i=1}^n X_i\), to prove the Theorem.

For any \(\lambda >0\) we have \[\begin{align} {\mathbb{P}}\left\{ \sum_{i=1}^n X_i \geq t \right\} &=& {\mathbb{P}}\{e^{\lambda \sum X_i} \geq e^{\lambda t} \} \\ &\leq & \frac{\mathbb{E}[e^{\lambda \sum X_i}]}{e^{\lambda t}} \\ &=& e^{-\lambda t} \prod_{i=1}^n \mathbb{E}[e^{\lambda X_i}] \end{align}\]

The following calculations reveal the source of the improvement over Hoeffding’s inequality. \[\begin{align} \mathbb{E}[e^{\lambda X_i}] &=& \mathbb{E}\left[1 + \lambda X_i + \sum_{m=2}^\infty \frac{\lambda ^m X_i^m}{m!} \right] \\ &\leq & 1 + \sum_{m=2}^\infty \frac{\lambda ^m a^{m-2} \sigma^2}{m!} \\ &=& 1 + \frac{\sigma^2}{a^2} \sum_{m=2}^\infty \frac{(\lambda a)^m}{m!} \\ &=& 1+ \frac{\sigma^2}{a^2} \left(e^{\lambda a} - 1 - \lambda a \right) \end{align}\] Therefore, \[{\mathbb{P}}\left\{ \sum_{i=1}^n X_i \geq t \right\} \leq e^{-\lambda t} \left[1+ \frac{\sigma^2}{a^2} \left(e^{\lambda a} - 1 - \lambda a \right) \right]^n\]

We will use a few simple inequalities (that can be easily proved with calculus) such as⁷ \(1+x \leq e^x,\; \text{for all } x\in\mathbb{R}\).

This means that, \[1 + \frac{\sigma^2}{a^2}\left(e^{\lambda a} - 1 - \lambda a \right) \leq e^{\frac{\sigma^2}{a^2}(e^{\lambda a} - 1 - \lambda a)},\]

which readily implies \[{\mathbb{P}}\left\{ \sum_{i=1}^n X_i \geq t \right\} \leq e^{-\lambda t} e^{\frac{n\sigma^2}{a^2}(e^{\lambda a} - 1 - \lambda a)}.\]

As before, we try to find the value of \(\lambda >0\) that minimizes \[\min_{\lambda}\left\{ -\lambda t + \frac{n\sigma^2}{a^2}(e^{\lambda a}-1-\lambda a) \right\}\] Differentiation gives \[-t + \frac{n\sigma^2}{a^2} (ae^{\lambda a} - a) = 0\] which implies that the optimal choice of \(\lambda\) is given by \[\lambda ^* = \frac{1}{a}\log\left(1 + \frac{at}{n\sigma^2}\right)\]

If we set \[u = \frac{at}{n\sigma^2},\] then \(\lambda^* = \frac{1}{a}\log(1 + u)\).

Now, the value of the minimum is given by \[-\lambda^* t + \frac{n\sigma^2}{a^2} (e^{\lambda^*a} -1 -\lambda^*a) = -\frac{n\sigma^2}{a^2}\left[(1+u)\log(1+u) - u \right].\]

This means that \[\begin{align} {\mathbb{P}}\left\{ \sum_{i=1}^n X_i \geq t \right\} &\leq& \exp\left(-\frac{n\sigma^2}{a^2}\left\{(1+u)\log(1+u)-u\right\} \right) \end{align}\]

The rest of the proof follows by noting that, for every \(u>0\), \[\left(1+u\right)\log(1+u) - u \geq \frac{u}{\frac{2}{u}+\frac{2}{3}},\] which implies:

\[\begin{align} {\mathbb{P}}\left\{ \sum_{i=1}^n X_i \geq t \right\} &\leq& \exp\left(-\frac{n\sigma^2}{a^2}\frac{u}{\frac{2}{u}+\frac{2}{3}} \right)\\ &=& \exp\left(-\frac{t^2}{2n\sigma^2 + \frac{2}{3}at} \right). \end{align}\]

\(\Box\)

We refer to [@vershynin2018] for several useful variations of Bernstein’s inequality.

3.3.2 The geometry of the hypercube revisited↩︎

Equipped with the probabilistic tools from the previous sections, we will give an alternative proof for Fact 3.2.

The proof of this statement will be based on a probabilistic argument, thereby illustrating (again) the nice and fruitful connection between geometry and probability in high dimension. Pick a point at random in \(C^d\), which corresponds to \(\left[-\frac{1}{2},\frac{1}{2}\right]^d\). We want to upperbound the probability that the point is also in the sphere or radius \(R\).

Let \(x = (x_1, \dots, x_d) \in \mathbb{R}^d\) and each \(x_i \in \left[-\frac{1}{2},\frac{1}{2}\right]\) is chosen uniformly at random. The event that \(x\) also lies in the sphere of radius \(R\) corresponds to the inequality \[\|x\|_2 = \sqrt{\sum_{i=1}^d x_i^2} \leq R.\] Let \(z_i = x_i^2\) and note that \[\begin{align} \mathbb{E}[z_i] = \int_{-1/2}^{1/2} t^2 \, dt = 2\frac{(1/2)^3}{3}=\frac{1}{12}, \end{align}\] which implies \(\mathbb{E}[\|x\|_2^2] = \frac{d}{12}\). Since the random variables \(z_i - \mathbb{E}[z_i]=z_i-\frac{1}{12}\) are independent, centered, and bounded (\(|z_i - \mathbb{E}[z_i]|\leq \max\left\{\frac{1}{12},\frac{1}{4}-\frac{1}{12}\right\}=\frac{2}{12}\)), we can use Hoeffding’s inequality (Theorem 3.10). This gives \[\begin{align} \mathbb{P}(\|x\|_2^2 \leq R^2) &= \mathbb{P}\left(\sum_{i=1}^d x_i^2 \leq R^2\right) = \mathbb{P}\left(\sum_{i=1}^d (z_i - \mathbb{E}[z_i]) \leq R^2 - \frac{d}{12}\right) \\ &\leq \exp\left[-\frac{\left(R^2-\frac{d}{12}\right)^2}{2d\left( \frac{2}{12} \right)^2}\right] = \exp\left[-\frac{\left(12R^2-d\right)^2}{8d}\right]. \end{align}\] Taking \(R=\frac{\sqrt{d}}{4}\) gives \(\mathbb{P}\left(\|x\|_2 \leq \frac{\sqrt{d}}{4}\right)\leq \exp\left[-c_0d\right]\) for a universal constant \(c_0\) (which can be computed in this case to be \(c_0=\frac{\left(1-\frac{12}{16}\right)^2}{8}=\frac{1}{128}\), but that’s not important).

Since we now have gained a better understanding of the properties of the cube in high dimensions, we can use this knowledge to our advantage. For instance, this “pointiness” of the hypercube (in form of the \(\ell_1\)-ball) turns out to very useful in the areas of compressive sensing and sparse recovery, see Chapter 16.

3.3.3 Tail bounds for the chi-squared distribution↩︎

If \(X_1,\dots,X_n\) are independent, standard normal random variables, then the sum of their squares, \(Z = \sum_{k=1}^n X_k^2\) is distributed according to the chi-squared distribution with \(n\) degrees of freedom. We denote this by \(Z \sim \chi^2(n)\). Its probability density function is \[\phi(t) = \begin{cases} \frac{t^{\frac{n}{2}-1} e^{-\frac{n}{2}}}{{ 2^{\frac{n}{2}} \Gamma(\frac{n}{2})}}, & t > 0. \\ 0, & \text{else.} \end{cases}\] Since the \(X_k^2, k=1,\dots,n\) are subexponential with parameters \((2,4)\) and independent, \(Z = \sum_{k=1}^n X_k^2\) is subexponential with parameters \((2\sqrt{n},4)\) (analogously to Remark 3.1). Therefore, using ?? , we obtain the \(\chi^2\) tail bound \[\label{chi-tail} \mathbb{P}\left( \left| \frac{1}{n} \sum_{k=1}^n X_k^2 - 1 \right| \ge t \right) \le \begin{cases} 2 e^{-nt^2/8} & \text{for t \in (0,1).} \\ 2e^{-nt/8} & \text{if t \ge 1.} \\ \end{cases}\tag{20}\] A very useful variation of this bound, commonly known as “Lemma 1 in Laurent and Massart [@Laurent_Massart_TailChiSquare]” is the following Theorem.

Note that if \(a_k=1\), for all \(k\), then \(Z\) is a \(\chi_2\) random variable with \(n\) degrees of freedom, so this theorem immediately gives a deviation inequality for \(\chi_2\) random variables, readily comparable to 20 .

3.3.4 How to generate random points on a sphere↩︎

How can we sample a point uniformly at random from \(S^{d-1}\)? The first approach that may come to mind is the following method to generate random points on a unit circle. Independently generate each coordinate uniformly at random from the interval \([-1,1]\). This yields points that are distributed uniformly at random in a square that contains the unit ball. We could now project all points onto the unit ball. However, the resulting distribution will not be uniform since more points fall on a line from the origin to a vertex of the square, than fall on a line from the origin to the midpoint of an edge due to the difference in length of the diagonal of the square to its side length.

To remedy this problem, we could discard all points outside the unit ball and project the remaining points onto the sphere. However, if we generalize this technique to higher dimensions, the analysis in the previous section has shown that the ratio of the volume of \(S^{d-1}(1)\) to the volume of \(C^d(1)\) decreases rapidly. This makes this process not practical, since almost all the generated points will be discarded in this process and we end up with essentially no points inside (and thus, after projection, on) the sphere.

Instead we can proceed as follows: Recall that the multivariate Gaussian distribution is rotation invariant about the origin. This means that we can simply construct a vector in \(\mathbb{R}^d\) whose entries are independently drawn from a univariate Gaussian distribution, and then normalize it to lie on the sphere. This gives a distribution of points that is uniform over the sphere.

Having a method of generating points uniformly at random on \(S^{d-1}\) at our disposal, we can now give a probabilistic proof that points on \(S^{d-1}\) concentrate near its equator.

Without loss of generality, we pick the first canonical basis vector \(e_1\) to represents the “north pole”, and the intersection of the sphere with the plane \(x_1=0\) forms our equator. We create a random vector by sampling \((z_1,\dots,z_d) \sim {\mathcal{N}}(0,I_d)\) and normalize the vector to get \(x = (x_1,\dots, x_d) = \frac{1}{\sqrt{\sum_{k=1}^d z_k^2}} (z_1,\dots,z_d)\). Because \(x\) is on the sphere, it holds that \(\|x\|_2=1\). We are interested in understanding the random variable \(x_1^2\). Notice that \(1=\sum_{k=1}^d x_k^2 = \mathbb{E}\left[\sum_{k=1}^d x_k^2\right] =d\,\mathbb{E}x_1^2\) since all coordinates are identically distributed (but not independent), this means that \(\mathbb{E}x_k^2=\frac{1}{d}\) for all coordinates \(x_k\). Our goal is to upperbound the probability that \(|x_1|=\frac{|z_1|}{\sqrt{\sum_{k=1}^d z_k^2}}\) is large, we will do that by union bounding the events that \(z_1\) is large (in absolute value) and that \(\sum_{k=1}^d z_k^2\) is small (which can be done with 20 ). We will not attempt to derive the strongest possible bound, but rather to showcase an useful argument. Let \(\rho>0\) and \(0<\varepsilon<1\) be parameters that we will set later (perhaps depending on \(d\)). We have

\[\begin{align} \mathbb{P}\left( |x_1| \geq \sqrt{\frac{\rho}{d}} \right) & \leq \mathbb{P}\left( |z_1| \geq \sqrt{(1-\varepsilon)\rho} \,\text{ or }\, \sum_{k=1}^d z_k^2 \leq (1-\varepsilon)d \right) \\ & \leq \mathbb{P}\left( |z_1| \geq \sqrt{(1-\varepsilon)\rho} \right) + \mathbb{P}\left( \frac{1}{d}\sum_{k=1}^d z_k^2 \leq 1 - \varepsilon\right) \\ & \leq 2\exp\left[ -\frac{(1-\varepsilon)\rho}{2} \right] + 2\exp\left[ -\frac{d\varepsilon^2}{2} \right], \end{align}\] where in the last inequality we used the Gaussian tail bound (Proposition 3.5) in the first term and the \(\chi^2\) bound 20 in the second.

For \(d\) large enough, taking \(\varepsilon\) sufficiently small we have \[\mathbb{P}\left( |x_1| \geq \sqrt{\frac{\rho}{d}} \right) \leq \exp\left[ -\frac{\rho}{3} \right],\] where the \(3\) could be replaced by any real number larger than \(2\).

A more general version of this phenomenon is that \(k\) coordinates of a uniformly random vector in \(\mathbb{S}^{d-1}\) have \(\ell_2\) norm squared that concentrates well around \(\frac{k}{d}\). We omit the the proof (it can be found in [@Dasgupta_Gupta_JLsimpleproof], and the ideas are similar to the ones used to obtain other inequalities in this Chapter).

We will dive deeper into the topic of probability in high dimensions in Chapters 14 and 15. There are several excellent texts about different aspects of high dimensional probability, two notable examples are [@vanHandel_LectureNotesProb_14] and [@vershynin2018].

Exercises↩︎

4 Singular Value Decomposition and Principal Component Analysis↩︎

Data is most often represented as a matrix, even network data and graphs are often naturally represented by their adjacency matrix. For this reason Linear Algebra is one of the key tools in data analysis. Perhaps more surprising is the fact that spectral properties of matrices representing data play a crucial role in data analysis. We illustrate this importance with a discussion of the Singular Value Decomposition (SVD) and Principal Component Analysis (PCA) that are often used for data compression, denoising, and dimension reduction. Tools from random matrix theory are then utilized to better understand the performance of SVD and PCA in the high dimensional regime.

4.1 Singular Value Decomposition↩︎

The SVD is like a mathematical Swiss Army knife, with its many uses across different areas including linear algebra, operator theory, functional analysis, optimization, numerical analysis, statistics, and information theory. The SVD is also one of the most useful tools for analyzing data. We recommend [@HJ90] and [@Golub_MatrixComputations] as base references in linear algebra and for the SVD.

Given a matrix \(A\in\mathbb{R}^{m\times n}\), the SVD of \(A\) is defined as \[\label{eq:SVD:PCAnotes} A = U \Sigma V^\top ,\tag{21}\] where \(U\in O(m)\), \(V\in O(n)\) are orthogonal matrices (meaning that \(U^\top U = UU^\top =\operatorname{I}_{m\times m}\) and \(V^\top V = VV^\top = \operatorname{I}_{n\times n}\)) and \(\Sigma\in\mathbb{R}^{m\times n}\) is a matrix with non-negative entries on its diagonal and otherwise zero entries.

The columns of \(U\) and \(V\), denoted by \(u_k\) and \(v_k\), are referred to as the left and right singular vectors of \(A\), respectively and the diagonal elements of \(\Sigma\) are the singular values of \(A\), denoted by \(\sigma_1 \geq \sigma_2 \geq \sigma_r > 0\). In particular, \(\operatorname{rank}(A)=r\), that is, the number of non-zero singular values \(r\) is the rank of \(A\). Through the SVD, any matrix can be written as a sum of rank-1 matrices \[\label{eq:svdrank1sum} A = \sum_{k=1}^r \sigma_k u_k v_k^\top .\tag{22}\]

Furthermore, it is straightforward to show that the left singular vectors \(u_1,\dots,u_m\) are eigenvectors of the matrix \(A A^\top\), the right singular vectors \(v_1,\dots,v_n\) are eigenvectors of the matrix \(A^\top A\), and the squared singular values \(\sigma_1,\dots,\sigma_r\) are the non-zero eigenvalues of the matrices \(A A^\top\) and \(A^\top A\), respectively.

4.1.1 Matrix norms↩︎

A very powerful modeling tool in data science is given by low rank matrices. In fact, we will devote the whole Chapter 17 to this topic. As already suggested in the expansion 22 the SVD will play an important role in this, being used to provide low rank approximation of data matrices.

In order to be able to talk about low rank approximations of matrices, we need a notion of distance between matrices. Just like with vectors, the distance between matrices can be measured using a suitable norm of the difference. One popular norm is the Frobenius norm, or the Hilbert-Schmidt norm, defined as \[\|A\|_F = \sqrt{\sum_{i,j} a_{ij}^2},\] where \(a_{ij}\) are the entries of the matrix \(A\). This norm is simply the Euclidean norm of a vector of length \(mn\) comprised of the matrix elements. The Frobenius norm can also be expressed in terms of the singular values. To see this, first express the Frobenius norm in terms of the trace of \(A^\top A\) as \[\label{eq:SVD:frob} \|A\|_F^2 = \sum_{i,j} a_{ij}^2 = \operatorname{Tr}(A^\top A),\tag{23}\] where we recall that the trace of a square matrix \(A\) is defined as \[\operatorname{Tr}(A) = \sum_{i} a_{ii}.\] A particularly important property of the trace is that for any \(A\) of size \(m \times n\) and \(B\) of size \(n\times m\) \[\label{eq:SVD:trace} \operatorname{Tr}(AB) = \operatorname{Tr}(BA).\tag{24}\] Note that this implies that, e.g., \(\operatorname{Tr}(ABC) = \operatorname{Tr}(CAB)\), but it does not imply that, e.g., \(\operatorname{Tr}(ABC) = \operatorname{Tr}(ACB)\) which is not true in general. Now, plugging the SVD (21 ) into (23 ) gives \[\label{eq:SVD:MF} \|A\|_F^2 = \operatorname{Tr}(A^\top A) = \operatorname{Tr}(V \Sigma^\top U^\top U \Sigma V^\top ) = \operatorname{Tr}(\Sigma^\top \Sigma) = \sum_{k=1}^r \sigma_k^2,\tag{25}\] where we used the orthogonality of \(U\) and \(V\) and the trace property (24 ). We conclude that the Frobenius norm of a matrix equals the Euclidean norm of its vector of singular values.

A different way to define the size of a matrix is by viewing it as an operator and measuring by how much it can dilate vectors. For example, the operator 2-norm is defined as \[\|A\|_2 = \sup_{\|x\|=1}\|A x\|.\] This operator norm can also be succinctly expressed in terms of the singular values. Indeed, for any \(x \in \mathbb{R}^n\) \[Ax = \sum_{k=1}^r \sigma_k u_k (v_k^\top x).\] Using the orthogonality of the left singular vectors \(u_k\) we get \[\|Ax \|^2 = \sum_{k=1}^r \sigma_k^2 \langle v_k , x \rangle^2 \leq \sigma_1^2 \sum_{k=1}^r \langle v_k, x \rangle^2 \leq \sigma_1^2 \sum_{k=1}^n \langle v_k, x \rangle^2 = \sigma_1^2 \|x\|^2,\] where the last equality is due to the orthogonality of the right singular vectors \(v_k\). Moreover, we get equality by choosing \(x=v_1\). We conclude that the 2-norm is simply the largest singular value \[\|A\|_2 = \sigma_1.\]

Before proceeding we briefly recall Hölder’s inequality

4.1.2 Existence of the SVD↩︎

The matrix 2-norm plays an important role in proving that any matrix has an SVD.

4.1.3 Low rank matrix approximation↩︎

A very important property of the SVD is that it provides the best low rank approximation of a matrix, when the approximation error is measured in terms of either the operator norm or the Frobenius norm. Specifically, for any \(0 \leq s \leq r\) consider the rank-\(s\) matrix \(A_s = \sum_{k=1}^s \sigma_k u_k v_k^\top\). Then, among all matrices of rank \(s\), \(A_s\) best approximates \(A\) in terms of the operator norm (or Frobenius norm) error. When measuring error using the operator norm, the approximation error equals \(\sigma_{s+1}\), which is the largest singular value among the remaining \(r-s\) smallest singular values \[\|A - A_s \|_2 = \inf_{B \in \mathbb{R}^{m\times n}, \operatorname{rank}(B) \leq s} \|A - B \|_2 = \sigma_{s+1}.\] A similar result holds for the Frobenius norm, with the approximation error given in terms of the remaining \(r-s\) smallest singular values as \[\|A - A_s \|_F = \inf_{B \in \mathbb{R}^{m\times n}, \operatorname{rank}(B) \leq s} \|A - B \|_F = \sqrt{\sum_{k=s+1}^r \sigma_k^2}.\] In fact, \(A_s\) is the best low rank approximation for any univariate matrix norm satisfying \(\|U A V\| = \| A \|\) for any \(U \in O(m), V\in O(n)\), that is, norms that are invariant to multiplication by orthogonal matrices [@mirsky1960symmetric].

We prove the low rank approximation property for the operator norm case⁸.

For the Frobenius norm case, the low rank approximation property follows from von Neumann’s trace inequality [@vonNeumann1937some]:

The reader is referred to [@mirsky1975trace] for a proof of Theorem 4.4. We are now ready to prove that the SVD provides the best low rank approximation in the Frobenius norm.

The low rank approximation property has a wide ranging implication on data compression. The storage size of an \(m\times n\) data matrix is \(mn\). If that matrix is of rank \(r\), then storage size reduces from \(mn\) to \((n+m+1)r\) (for storing \(r\) left and right singular vectors and values). For \(r \ll \min\{n,m\}\) this reduction can be quite dramatic. For example, if \(r=10\) and \(n=m=10^6\), then storage reduces from \(10^{12}\) entries to just \(2\cdot 10^7\). But even if the matrix is not precisely of rank \(r\), but only approximately, in the sense that \(\sigma_{r+1} \ll \sigma_1\), then we are guaranteed by the above approximation results to incur only a small approximation error due to compression using the top \(r\) singular vectors and values. In many cases, the singular values of large data matrices decay very quickly, motivating this type of low rank approximation which oftentimes is the only way to handle massive data sets that otherwise cannot be stored and/or manipulated efficiently. Remarkably, even treating an image as a matrix of pixel intensity values and compressing it this way yields good image compression and de-noising algorithms (as it mitigates the noise corresponding to singular values that are truncated). This is illustrated using a photo of size \(1224 \times 1632\) pixels taken at Prospect Garden in Princeton University on April 28, 2015 by Amit Singer.

The original photo, see Figure 5 (a), is considered as three matrices (one for each RGB color channel). The main features of the image appear already in the rank-10 approximation, Figure 5 (b), although many fine details are clearly missing and there are also some visible artifacts (such as striping, smearing, etc.). The rank-100 approximation, Figure 5 (c), has exceptional resemblance with the original photo. Of course there are still differences between the two images (such as artifacts on the curbstone) but these are less visible to the naked eye on first inspection. Note that compression using \(r=10\) and \(r=100\) require less than \(1\%\) and \(10\%\), respectively, storage of the original image. The decay of the singular values of the three matrices is illustrated in Figure 6. We use a logarithmic scale due to the rapid decay of singular values. Notice that the singular values of the three matrices are very similar.

4.1.4 Spectral Decomposition↩︎

If \(M\in\mathbb{R}^{n\times n}\) is symmetric then it admits a spectral decomposition \[M = V\Lambda V^\top ,\] where \(V\in O(n)\) is a matrix whose columns \(v_k\) are the eigenvectors of \(M\) and \(\Lambda\) is a diagonal matrix whose diagonal elements \(\lambda_k\) are the eigenvalues of \(M\). Similarly, we can write \[M = \sum_{k=1}^n\lambda_k v_k v_k^\top .\]

When all of the eigenvalues of \(M\) are non-negative we say that \(M\) is positive semidefinite and write \(M\succeq 0\). In that case we can write \[M = \left( V \Lambda^{1/2} \right) \left( V \Lambda^{1/2} \right)^\top .\] A decomposition of \(M\) of the form \(M=UU^\top\) (such as the one above) is called a Cholesky decomposition (oftentimes, in a Cholesky decomposition, the matrix \(U\) is required to be triangular, here we do not make this requirement).

For symmetric matrices, the operator 2-norm is also known as the spectral norm, given by \[\|M\| = \max_k \left|\lambda_k(M)\right|.\]

4.1.5 Quadratic forms↩︎

In both this and following chapters, we will be interested in solving problems of the type \[\max_{\substack{V\in \mathbb{R}^{n\times d} \\ V^\top V = \operatorname{I}_{d\times d}}} \operatorname{Tr}\left(V^\top M V \right),\] where \(M\) is a symmetric \(n\times n\) matrix.

Note that this is equivalent to \[\label{eq:0:evalue:vard} \max_{\substack{v_1,\dots,v_d \in \mathbb{R}^{n} \\ v_i^\top v_j = \delta_{ij}}} \sum_{k=1}^d v_k^\top M v_k,\tag{26}\] where \(\delta\) is the Kronecker delta (\(\delta_{ij}=1\) for \(i=j\) and \(\delta_{ij}=0\) otherwise).

When \(d=1\) this reduces to the more familiar \[\label{eq:0:evalue:var} \max_{\substack{v\in \mathbb{R}^n \\ \|v\|_2 = 1}} v^\top M v.\tag{27}\]

It is easy to see (for example, using the spectral decomposition of \(M\)) that 27 is maximized by the leading eigenvector of \(M\) and \[\max_{\substack{v\in \mathbb{R}^n \\ \|v\|_2 = 1}} v^\top M v = \lambda_{\max} (M).\]

Furthermore 26 is maximized by taking \(v_1,\dots,v_d\) to be the \(d\) leading eigenvectors of \(M\) and its value is simply the sum of the \(d\) largest eigenvalues of \(M\). This follows, for example, from a Theorem of Fan (see page 3 of [@Fan_variational_evalue_UAU]). A fortunate consequence is that the solution to 26 can be computed sequentially: we can first solve for \(d=1\), computing \(v_1\), then update the solution for \(d=2\) by simply computing \(v_2\) without changing the first vector.

4.1.6 SVD and the Moore-Penrose Pseudoinverse↩︎

The SVD is intimately connected to the Moore-Penrose pseudoinverse of a matrix [@ben2003generalized]. The Moore-Penrose pseudoinverse (named after mathematicians E. H. Moore and Roger Penrose) of \(A\), denoted as \(A^\dagger\), is a generalization of the matrix inverse to matrices that may not be square, or if square, may not have full rank. For an invertible matrix, the pseudoinverse is identical to the standard inverse. However, for a rectangular or rank-deficient square matrix, the pseudoinverse provides a unique solution that satisfies a set of four specific conditions known as the Penrose conditions, that uniquely define the pseudoinverse \(A^\dagger\) for any matrix \(A\). These conditions are a generalization of the properties of a standard matrix inverse. A matrix \(A^\dagger\) is the pseudoinverse of \(A\) if and only if it satisfies all four of the following equations: \[\begin{align} \label{eq:pinvproperties} \begin{aligned} AA^\dagger A & = A, \\ A^\dagger AA^\dagger & = A^\dagger, \\ (AA^{\dagger})^{T} & = AA^\dagger, \\ (A^{\dagger} A)^\top & = A^\dagger A. \end{aligned} \end{align}\tag{28}\] These conditions ensure that the pseudoinverse acts “as much as possible” as an inverse for matrices that are not invertible, including rectangular or rank-deficient square matrices. The above four conditions are equivalent to the matrices \(A A^\dagger\) and \(A^\dagger A\) being orthogonal projectors onto the range of \(A\) and the range of \(A^\top\) respectively.

The primary purpose of the pseudoinverse is to find “best-fit” solution to a system of linear equations \(Ax=b\) when a unique solution does not exist. In cases where the system has multiple solutions, the pseudoinverse provides the one with the minimum Euclidean norm. For systems with no exact solution (an overdetermined system), it provides the least-squares solution, which minimizes the error \(\|Ax - b\|\). This makes the pseudoinverse a tremendously useful tool in fields like numerical mathematics, statistics, machine learning, and signal processing. We will explore this in more detail for example in Chapter 5.5.

The SVD provides a straightforward and computationally stable way to calculate \(A^\dagger\). Given \(A = U \Sigma V^\top\), we can represent \(\Sigma\) in block form as \[\Sigma = \begin{bmatrix} \Sigma_r & \mathbf{0}_{r \times (n-r)} \\ \mathbf{0}_{(m-r) \times r} & \mathbf{0}_{(m-r) \times (n-r)} \end{bmatrix},\] where \(\Sigma_r = \text{diag}(\sigma_1, \sigma_2, \dots, \sigma_r)\) is an \(r \times r\) diagonal matrix containing the non-zero singular values. The pseudoinverse of \(A\) is given by the formula. \[\label{pinv} A^\dagger = V \Sigma^\dagger U^\top ,\tag{29}\] where the matrix \(\Sigma^\dagger\) is given by \[\Sigma^\dagger = \begin{bmatrix} \Sigma_r^{-1} & \mathbf{0}_{r \times (m-r)} \\ \mathbf{0}_{(n-r) \times r} & \mathbf{0}_{(n-r) \times (m-r)} \end{bmatrix},\] with \(\Sigma_r^{-1} = \text{diag}\left(\frac{1}{\sigma_1}, \frac{1}{\sigma_2}, \dots, \frac{1}{\sigma_r}\right)\).

This method is particularly powerful because it works for any matrix—square or rectangular, full rank or rank-deficient—and is numerically robust, avoiding the issues with ill-conditioned matrices that can plague other inversion methods. Indeed, in the practical computation of \(\Sigma^\dagger\) one usually chooses a threshold \(\tau\), such that all singular values less than \(\tau\) are treated as zero.

The SVD and the pseudoinverse play a central role in linear regression and least squares problems, see Chapter 5.

4.2 Principal Component Analysis and dimension reduction↩︎

When faced with a high dimensional dataset, a natural approach is to attempt to reduce its dimension, either by projecting it to a lower dimensional space or by finding a better representation for the data using a small number of meaningful features. Beyond data compression and visualization, dimension reduction facilitates downstream analysis such as clustering and regression that perform significantly better in lower dimensions. We will explore a few different ways of reducing the dimension, both linearly and non-linearly.

We start with the classical Principal Component Analysis (PCA). PCA continues to be one of the most effective and simplest tools for exploratory data analysis. It dates back to a 1901 paper by Karl Pearson [@Pearson_PCA_1901] and was rediscovered several times by different research communities where it is also known under different names such as the proper orthogonal decomposition (POD) and the Karhunen–Loève transform (KLT) among others.

Suppose we have \(n\) data points \(x_1,\dots,x_n\) in \(\mathbb{R}^p\), and we are interested in (linearly) projecting the data to \(d< p\) dimensions. This is particularly useful if, say, one wants to visualize the data in two or three dimensions (\(d=2, 3\)). There are a couple of seemingly different criteria we can use to choose this projection:

1. Find the \(d\)-dimensional affine subspace for which the projections of \(x_1,\dots,x_n\) on it best approximate the original points \(x_1,\dots,x_n\).

2. Find the \(d\)-dimensional projection of \(x_1,\dots,x_n\) that preserves as much variance of the data as possible.

As we will see below, these two approaches are equivalent and they correspond to Principal Component Analysis.

Before proceeding, we recall a couple of simple statistical quantities associated with \(x_1,\dots,x_n\), that will reappear below.

Given \(x_1,\dots,x_n\) we define its sample mean as \[\mu_n = \frac{1}{n}\sum_{k=1}^n x_k,\] and its sample covariance as \[\Sigma_n = \frac{1}{n-1}\sum_{k=1}^n \left(x_k-\mu_n\right)\left(x_k-\mu_n\right)^\top .\]

4.2.0.1 PCA as the best \(d\)-dimensional affine fit

We start with the first interpretation of PCA and then show that it is equivalent to the second. We are trying to approximate each \(x_k\) by \[\label{eq:PCA:01} x_k \approx \mu + \sum_{i=1}^d \left(\beta_k\right)_iv_i,\tag{30}\] where \(v_1,\dots,v_d\) is an orthonormal basis for the \(d\)-dimensional subspace, \(\mu \in \mathbb{R}^p\) represents the translation, and \(\beta_k \in \mathbb{R}^d\) corresponds to the coefficients of \(x_k\). Without loss of generality we take \[\label{eq:bk95sum950} \sum_{k=1}^n \beta_k = 0,\tag{31}\] as any joint translation of \(\beta_k\) can be absorbed into \(\mu\).

If we represent the subspace by \(V = [v_1\cdots v_d]\in\mathbb{R}^{p\times d}\) then we can rewrite 30 as \[\label{nxfcjldz} x_k \approx \mu + V\beta_k,\tag{32}\] where \(V^\top V=\operatorname{I}_{d\times d}\), because the vectors \(v_i\) are orthonormal.

We will measure goodness of fit in terms of least squares and attempt to solve \[\label{eq:PCA:goodnessoffit0} \min_{\substack{\mu,\;V,\;\beta_k \\ V^\top V=\operatorname{I}}} \sum_{k=1}^n \left\| x_k - \left( \mu + V\beta_k \right) \right\|_2^2\tag{33}\]

We start by optimizing for \(\mu\). It is easy to see that the first order condition for \(\mu\) corresponds to \[\nabla_\mu \sum_{k=1}^n \left\| x_k - \left( \mu + V\beta_k \right) \right\|_2^2 = 0 \Longleftrightarrow \sum_{k=1}^n \left( x_k - \left( \mu + V\beta_k \right) \right) = 0.\] Thus, the optimal value \(\mu^\ast\) of \(\mu\) satisfies \[\left( \sum_{k=1}^n x_k \right) - n\mu^\ast - V\left( \sum_{k=1}^n \beta_k \right)=0.\] Since we assumed in 31 that \(\sum_{k=1}^n \beta_k = 0\), we have that the optimal \(\mu\) is given by \[\mu^\ast = \frac{1}{n} \sum_{k=1}^n x_k = \mu_n,\] the sample mean.

We can then proceed to finding the solution for 33 by solving \[\label{eq:PCA:goodnessoffit1} \min_{\substack{V,\;\beta_k \\ V^\top V=\operatorname{I}}} \sum_{k=1}^n \left\| x_k - \mu_n - V\beta_k \right\|_2^2.\tag{34}\]

Let us proceed by optimizing for \(\beta_k\). The problem almost fully decouples in each \(k\), the only constraint coupling them being 31 . We will ignore this constraint, solve the decoupled problems, and verify that it is automatically satisfied. Hence we focus on, for each \(k\), \[\label{eq:PCA:goodnessoffit12} \min_{\beta_k} \left\| x_k - \mu_n - V\beta_k \right\|_2^2 = \min_{\beta_k} \Big\| x_k - \mu_n - \sum_{i=1}^d \left(\beta_k \right)_i v_i \Big\|_2^2.\tag{35}\] Since \(v_1,\dots,v_d\) are orthonormal, it is easy to see that the solution is given by \(\left(\beta_k^\ast \right)_i = v_i^\top \left( x_k - \mu_n\right)\) which can be succinctly written as \(\beta_k^* = V^\top \left(x_k - \mu_n\right)\), which satisfies 31 . Thus, 34 is equivalent to \[\label{eq:PCA:goodnessoffit13} \min_{V^\top V=\operatorname{I}} \sum_{k=1}^n \left\| \left( x_k - \mu_n \right) - VV^\top \left(x_k - \mu_n\right) \right\|_2^2.\tag{36}\]

Note that \[\begin{align} \left\| \left( x_k - \mu_n \right) - VV^\top \left(x_k - \mu_n\right) \right\|_2^2 & = & \left( x_k - \mu_n \right)^\top \left( x_k - \mu_n \right) \\ & & - 2\left( x_k - \mu_n \right)^\top VV^\top \left(x_k - \mu_n\right) \\ & & +\left(x_k - \mu_n\right)^\top V\left(V^\top V\right)V^\top \left(x_k - \mu_n\right) \\ & = & \left( x_k - \mu_n \right)^\top \left( x_k - \mu_n \right) \\ & & - \left( x_k - \mu_n \right)^\top VV^\top \left(x_k - \mu_n\right). \end{align}\]

Since \(\left( x_k - \mu_n \right)^\top \left( x_k - \mu_n \right)\) does not depend on \(V\), minimizing 36 is equivalent to \[\label{eq:PCA:goodnessoffit1395max} \max_{V^\top V=\operatorname{I}} \sum_{k=1}^n \left( x_k - \mu_n \right)^\top VV^\top \left(x_k - \mu_n\right).\tag{37}\]

A few algebraic manipulations using properties of the trace yields: \[\begin{align} \sum_{k=1}^n \left( x_k - \mu_n \right)^\top VV^\top \left(x_k - \mu_n\right) & = &\sum_{k=1}^n \operatorname{Tr}\left[ \left( x_k - \mu_n \right)^\top VV^\top \left(x_k - \mu_n\right) \right] \\ & = & \sum_{k=1}^n\operatorname{Tr}\left[ V^\top \left(x_k - \mu_n\right)\left( x_k - \mu_n \right)^\top V \right] \\ & = & \operatorname{Tr}\left[ V^\top \sum_{k=1}^n\left(x_k - \mu_n\right)\left( x_k - \mu_n \right)^\top V \right] \\ & = & (n-1)\operatorname{Tr}\left[ V^\top \Sigma_n V \right]. \end{align}\] This means that the solution to 37 is given by the solution of \[\label{eq:PCA:goodnessoffit1395min} \max_{V^\top V=\operatorname{I}} \operatorname{Tr}\left[ V^\top \Sigma_n V \right].\tag{38}\] As we saw above (recall 26 ) the solution is given by \(V=\left[v_1,\cdots,v_d\right]\) where \(v_1,\dots,v_d\) correspond to the \(d\) leading eigenvectors of \(\Sigma_n\).

4.2.0.2 PCA as the \(d\)-dimensional projection that preserves the most variance

We now show that the alternative interpretation of PCA, of finding the \(d\)-dimensional projection of \(x_1,\dots,x_n\) that preserves the most variance, also arrives to the optimization problem 38 . We aim to find an orthonormal basis \(v_1,\dots,v_d\) (organized as \(V=\left[v_1,\dots,v_d\right]\) with \(V^\top V = \operatorname{I}_{d\times d}\)) of a \(d\)-dimensional space such that the projection of \(x_1,\dots,x_n\) onto this subspace has the most variance. Equivalently we can ask for the points \[\left\{ \left[\begin{array}{c} v_1^\top x_k \\ \vdots \\ v_d^\top x_k\end{array}\right] \right\}_{k=1}^n,\] to have as much variance as possible. Hence, we are interested in solving \[\label{eq:PCA:maxvariance95min} \max_{V^\top V=\operatorname{I}} \sum_{k=1}^n\left\| V^\top x_k - \frac{1}{n} \sum_{r=1}^n V^\top x_r \right\|^2.\tag{39}\] Note that \[\sum_{k=1}^n\left\| V^\top x_k- \frac{1}{n} \sum_{r=1}^n V^\top x_r \right\|^2 = \sum_{k=1}^n\left\| V^\top \left(x_k-\mu_n\right) \right\|^2 = (n-1) \operatorname{Tr}\left(V^\top \Sigma_n V\right),\] showing that 39 is equivalent to 38 and that the two interpretations of PCA are indeed equivalent.

4.2.0.3 Finding the Principal Components

When given a dataset \(x_1,\dots,x_n\in\mathbb{R}^p\), in order to compute the Principal Components one needs to compute the leading eigenvectors of \[\Sigma_n = \frac{1}{n-1}\sum_{k=1}^n \left(x_k-\mu_n\right)\left(x_k-\mu_n\right)^\top .\] A naïve way of doing this is to construct \(\Sigma_n\) (which takes \(\mathcal{O}(np^2)\) work) and then finding its spectral decomposition (which takes \(\mathcal{O}(p^3)\) work). This means that the computational complexity of this procedure is \(\mathcal{O}\left( \max \left\{ np^2, p^3 \right\} \right)\) (e.g.see [@HJ90] or [@Golub_MatrixComputations]).

An alternative is to use the Singular Value Decomposition 21 . Let \(X = \left[x_1 \cdots x_n \right]\), and recall that \[\Sigma_n = \frac{1}{n-1} \left( X - \mu_n\mathbf{1}^\top \right) \left( X - \mu_n\mathbf{1}^\top \right)^\top .\]

Let us take the SVD of \(\frac{1}{\sqrt{n-1}}\left(X - \mu_n\mathbf{1}^\top \right) = U_L D U_R^\top\) with \(U_L \in O(p)\), \(D\) diagonal, and \(U_R^\top U_R = \operatorname{I}\). Then, \[\Sigma_n = \frac{1}{n-1} \left( X - \mu_n\mathbf{1}^\top \right) \left( X - \mu_n\mathbf{1}^\top \right)^\top = U_L D U_R^\top U_R D U_L^\top = U_L D^2 U_L^\top ,\] meaning that \(U_L\) corresponds to the eigenvectors of \(\Sigma_n\). Computing the SVD of \(X - \mu_n\mathbf{1}^\top\) takes \(\mathcal{O}(\min \{n^2p,p^2n\})\) work but if one is interested in simply computing the top \(d\) eigenvectors then this computational cost reduces to \(\mathcal{O}(dnp)\). This can be further improved with randomized algorithms. In Chapter 8.2 we will discuss randomized algorithms that compute an approximate solution in \(\mathcal{O}\left( pn\log d +(p+n)d^2 \right)\) time, see also [@HMT11; @Rokhlin_Szlam_Tygert_randPCA; @Musco_Musco_BlockLanczos].

Numerical stability is another important reason why computing the principal components using the SVD is preferable. Since the eigenvalues of \(\Sigma_n\) are proportional to the squares of the singular values of \(X-\mu_n \mathbf{1}^\top\), problems arise when the ratio of singular values exceeds \(10^8\), causing the ratio of the corresponding eigenvalues of \(\Sigma_n\) to be larger than \(10^{16}\). In this case, the smaller eigenvalue would be rounded to zero (due to machine precision), which is certainly not desirable.

4.2.0.4 Example: The eigenspace of handwritten digits (MNIST)

Every aspiring data scientist will at some point experiment with the MNIST dataset, a collection of handwritten digits created by researchers Burges, Cortes, and LeCun from data by the National Institute of Standards and Technology (NIST) [@lecun1998mnist]. The most common version of MNIST has 60,000 training images and 10,000 test images, each image containing \(28 \times 28\) grayscale pixels. MNIST has become a highly influential, frequently used benchmark in the data science/machine learning community⁹.

We will apply PCA to the MNIST dataset to analyze whether a high-dimensional space (784 dimensions in this case) can be compressed into a low-dimensional latent space that is still able to capture the “essence” of these handwritten digits. To process the data via PCA, we first flatten each image into a vector \(x\in\mathbb{R}^{784}\). The resulting data matrix \(X\) has a size of \(60000\times784\). In this matrix, each row is a single handwritten digit (an observation), while each column represents the intensity of a specific pixel location (a feature).

Before performing PCA, we compute the mean digit \(\bar{x}\), which represents the “average” handwriting across the entire dataset. By subtracting this mean from every row \(X_{\text{centered}} =X - \bar{x}\), we shift the coordinate system so that the origin lies at the center of the data cloud. This ensures that the first principal component truly represents the direction of maximum variance. When we compute the SVD of \(X_{\text{centered}}\), the right singular vectors (the columns of V) are the Principal Components. We can reshape these \(784 \times 1\) vectors back into \(28 \times 28\) images, sometimes called eigen-digits, to visualize what the algorithm has “learned.” The first twelve eigen-digits are depicted in Figure 7.

4.2.0.5 Which \(d\) should we pick?

Given a dataset, if the objective is to visualize it then picking \(d=2\) or \(d=3\) might make the most sense. However, PCA is useful for many other purposes, for example:

1. Denoising: often times the data belongs to a lower dimensional space but is corrupted by high dimensional noise. In such cases, PCA helps reduce the noise while keeping the signal.

2. Downstream analysis: One may be interested in running an algorithm (clustering, regression, etc.) that would be too computationally expensive or too statistically insignificant to run in high dimensions. Dimension reduction using PCA may help there.

In these applications (and many others) it is not clear how to pick \(d\). A fairly popular heuristic is to try to choose the cut-off at a component that has significantly more variance than the one immediately after. Since the total variance is \(\operatorname{Tr}(\Sigma_n) = \sum_{k=1}^p \lambda_k\), the proportion of variance in the \(i\)’th component is nothing but \(\frac{\lambda_i}{\operatorname{Tr}(\Sigma_n)}\). A plot of the values of the ordered eigenvalues, also known as a scree plot, helps to identify a reasonable choice of \(d\). It is common to then try to identify an “elbow” on the scree plot to choose the cut-off, as illustrated in Figure 8 (a). Figure 8 (b) shows the scree plot for the MNIST example, where the elbow is less well defined, as often happens in practice. In this case, the point of maximum curvature of the graph can serve as a surrogate for the elbow.

Figure 9 shows the choice of different values of \(k=d\) for the PCA approximation of MNIST, illustrated via some typical example digits.

In the next section we will look into Random Matrix Theory to better understand the behavior of the eigenvalues of \(\Sigma_n\) and gain insight into choosing cut-off values.

4.3 PCA in high dimensions and Marčenko-Pastur law↩︎

Let us assume that the data points \(x_1,\dots,x_n\in\mathbb{R}^p\) are independent draws of a zero mean Gaussian random variable \(g\sim \mathcal{N}(0,\Sigma)\) with some covariance matrix \(\Sigma \in \mathbb{R}^{p\times p}\). In this case, when we use PCA we are hoping to find a low dimensional structure in the distribution, which should correspond to the large eigenvalues of \(\Sigma\) (and their corresponding eigenvectors). For that reason, and since PCA depends on the spectral properties of \(\Sigma_n\), we would like to understand whether the spectral properties of the sample covariance matrix \(\Sigma_n\) (eigenvalues and eigenvectors) are close to those of \(\Sigma\), also known as the population covariance.

Since \(\mathbb{E}\Sigma_n = \Sigma\), if \(p\) is fixed and \(n\to \infty\) the law of large numbers guarantees that indeed \(\Sigma_n\to \Sigma\). However, in many modern applications it is not uncommon to have \(p\) of the same order as \(n\) (or, sometimes, even larger). For example, if our dataset is composed of images then \(n\) is the number of images and \(p\) the number of pixels per image; it is conceivable that the number of pixels and the number of images are comparable. Unfortunately, in that case, it is no longer clear that \(\Sigma_n\to \Sigma\). Dealing with this type of difficulties is an important goal of high dimensional statistics.

For simplicity we will try to understand the spectral properties of \[S_n = \frac{1}{n} X X^\top ,\] where \(x_1,\dots,x_n\) are the columns of \(X\). Since \(x\sim \mathcal{N}(0,\Sigma)\) we know that \(\mu_n \to 0\) (and, clearly, \(\frac{n}{n-1}\to 1\)), hence the spectral properties of \(S_n\) will be essentially the same as \(\Sigma_n\).¹⁰

Let us start by looking into a simple example, \(\Sigma = \operatorname{I}\). In that case, the distribution has no low dimensional structure, as the distribution is rotation invariant. Figure 10 depicts a normalized histogram (left panel) and a scree plot (right panel) of the eigenvalues of a realization of \(S_n\) (when \(\Sigma = \operatorname{I}\)) for \(p = 500\) and \(n=1000\). The red line is the eigenvalue distribution predicted by the Marčenko-Pastur distribution 40 , that we will discuss below.

As one can see, there are many eigenvalues considerably larger than \(1\), as well as many eigenvalues significantly smaller than \(1\). Notice that, if given this profile of eigenvalues of \(\Sigma_n\) one could potentially be led to believe that the data has low dimensional structure, when in truth the distribution it was drawn from is isotropic.

Understanding the distribution of eigenvalues of random matrices is at the core of Random Matrix Theory (there are many good books on Random Matrix Theory, e.g. [@bai2010spectral; @Tao_topicsRMT; @Anderson_Guionnet_Zeitouni_IntroRandomMatrices]). This particular limiting distribution was first established in 1967 by Marčenko and Pastur [@Marchenko_Pastur_1967] and is now referred to as the Marčenko-Pastur distribution. They showed that, if \(p\) and \(n\) both grow indefinitely to \(\infty\) with their ratio fixed \(p/n = \gamma \leq 1\), the sample distribution of the eigenvalues of \(S_n\) (like the histogram above), in the limit, will be

\[\label{eq:MPdist95PCA1} dF_\gamma(\lambda) = \frac{1}{2\pi} \frac{\sqrt{\left( \gamma_+ - \lambda \right)\left( \lambda - \gamma_- \right)}}{\gamma \lambda} 1_{\left[ \gamma_-, \gamma_+ \right]}(\lambda) d\lambda,\tag{40}\] with support \(\left[ \gamma_-, \gamma_+ \right]\), where \(\gamma_- = (1-\sqrt{\gamma})^2\), \(\gamma_+= (1+\sqrt{\gamma})^2\), and \(\gamma = p/n\). This is plotted as the red line in the figure above.

4.3.1 Spike models and the BBP phase transition↩︎

What if there actually is some (linear) low dimensional structure in the data? When can we expect to capture it with PCA? A particularly simple, yet relevant, example to analyze is when the covariance matrix \(\Sigma\) is a rank-1 perturbation of the identity matrix and is of the form \(\Sigma = I + \beta uu^\top\), where \(u\) is a unit norm vector and \(\beta > 0\). This is a particular case of a spike model.

One way to think about this spike model is that each data point \(x\) consists of a signal part \(\sqrt{\beta}g_0 u\) where \(g_0\) is a scalar standard Gaussian \(\mathcal{N}(0,1)\) (i.e. a normally distributed multiple of a fixed vector \(\sqrt{\beta}u\)) and a noise part \(g\sim\mathcal{N}(0,I)\) (independent of \(g_0\)). Then \(x = g + \sqrt{\beta}g_0 u\) is a Gaussian random variable \[\label{eq:spikymodel} x\sim \mathcal{N}(0,I + \beta uu^\top ).\tag{41}\] Whereas the signal part \(\sqrt{\beta}g_0 u\) resides on a central line in the direction of \(u\), the noise part is high dimensional and isotropic. We therefore refer to \(\beta\) as the signal-to-noise ratio (SNR). Indeed, \(\beta\) is the ratio of the signal variance (in the \(u\)-direction) to the noise variance (in each direction).

A natural question is whether this rank-\(1\) perturbation can be seen in \(S_n = \frac{1}{n} XX^\top\), where the columns of \(X\) are now given by 41 . In other words, is it possible to detect the direction of the line \(u\) from noisy measurements in high dimension? Let us build some intuition with an example. Figure 11 shows the histogram of the eigenvalues of a random realization of \(S_n\) for \(p=500\), \(n=1000\), \(u\) is the first element of the canonical basis \(u=e_1\), and \(\beta = 1.5\):

The histogram suggests that there is an eigenvalue of \(S_n\) that “pops out” of the support of the Marčenko-Pastur distribution (below we will estimate the location of this eigenvalue, and that estimate corresponds to the red “x”). It is worth noting that the largest eigenvalues of \(\Sigma\) is simply \(1+\beta = 2.5\) while the largest eigenvalue of \(S_n\) appears considerably larger than that.

Let us try now the same experiment with \(\beta = 0.5\), see Figure 12 for the result. It appears that, for \(\beta = 0.5\), the histogram of the eigenvalues is indistinguishable from when \(\Sigma = \operatorname{I}\). In particular, no eigenvalue is separated from the Marčenko-Pastur distribution.

This motivates the following question:

Question: For which values of \(\gamma\) and \(\beta\) do we expect to see an eigenvalue of \(S_n\) popping out of the support of the Marčenko-Pastur distribution, and what is the limiting value that we expect it to take?

As we will see below, there is a critical value of \(\beta\), denoted \(\beta_c\), below which we do not expect to see a change in the distribution of eigenvalues and above which we expect one of the eigenvalues to pop outside of the support. This phenomenon is known as the BBP phase transition (after Baik, Ben Arous, and Péché [@BPP_MC_BPP_2005]). There are many very nice papers about this and similar phenomena, including [@Paul_Marchenko_Pastur_BPP_TechReport; @Johnston_MC_BPP; @BPP_MC_BPP_2005; @Paul_Marchenko_Pastur_BPP_ActualPaper; @Baik_Silverstein_BPP_05; @Karoui_BPP_05; @BenaychGeorges_Nadakuditi_PCA; @BenaychGeorges_Nadakuditi_SVD].¹¹

In what follows we will find the critical value \(\beta_c\) and estimate the location of the largest eigenvalue of \(S_n\) for any \(\beta\). While the argument we will use can be made precise (and is borrowed from [@Paul_Marchenko_Pastur_BPP_TechReport]) we will be ignoring a few details for the sake of exposition. In other words, the argument below can be transformed into a rigorous proof, but it is not one at the present form.

We want to understand the behavior of the leading eigenvalue of the sample covariance matrix \[S_n = \frac{1}{n}\sum_{i=1}^n x_ix_i^\top .\] Since \(x \sim \mathcal{N}(0,I+\beta uu^\top )\) we can write \(x = (I + \beta uu^\top )^{1/2} z\) where \(z \sim \mathcal{N}(0,I)\) is an isotropic Gaussian. Then, \[S_n = \frac{1}{n} \sum_{i=1}^n (I + \beta uu^\top )^{1/2} z_i z_i^\top (I + \beta uu^\top )^{1/2} = (I + \beta uu^\top )^{1/2} Z_n (I + \beta uu^\top )^{1/2},\] where \(Z_n = \frac{1}{n} \sum_{i=1}^n z_i z_i^\top\) is the sample covariance matrix of independent isotropic Gaussians. The matrices \(S_n = (I + \beta uu^\top )^{1/2} Z_n (I + \beta uu^\top )^{1/2}\) and \(Z_n (I + \beta uu^\top )\) are related by a similarity transformation, and therefore have exactly the same eigenvalues. Hence, it suffices to find the leading eigenvalue of the matrix \(Z_n(I + \beta uu^\top )\), which is a rank-1 perturbation of \(Z_n\) (indeed, \(Z_n(I + \beta uu^\top ) = Z_n + \beta Z_n uu^\top\)). We already know that the eigenvalues of \(Z_n\) follow the Marčenko-Pastur distribution, so we are left to understand the effect of a rank-1 perturbation on its eigenvalues.

To find the leading eigenvalue \(\lambda\) of \(Z_n(I + \beta uu^\top )\), let \(v\) be the corresponding eigenvector, that is, \[Z_n(I + \beta uu^\top ) v = \lambda v.\] Subtract \(Z_n v\) from both sides to get \[\beta Z_n uu^\top v = (\lambda I - Z_n)v.\] Assuming \(\lambda\) is not an eigenvalue of \(Z_n\), we can multiply by \((\lambda I - Z_n)^{-1}\) to get¹² \[\beta (\lambda I - Z_n)^{-1} Z_n uu^\top v = v.\] Our assumption also implies that \(u^\top v \neq 0\), for otherwise \(v=0\). Multiplying by \(u^\top\) gives \[\beta u^\top (\lambda I - Z_n)^{-1} Z_n u (u^\top v) = u^\top v.\] Dividing by \(\beta u^\top v\) (which is not 0 as explained above) yields \[\label{eq:SVD:spike1} u^\top (\lambda I - Z_n)^{-1} Z_n u = \frac{1}{\beta}.\tag{42}\] Suppose \(w_1,\ldots,w_p\) are orthonormal eigenvectors of \(Z_n\) (with corresponding eigenvalues \(\lambda_1,\ldots, \lambda_p\)), and expand \(u\) in that basis: \[u = \sum_{i=1}^p \alpha_i w_i.\] Plugging this expansion in (42 ) gives \[\sum_{i=1}^p \frac{\lambda_i}{\lambda-\lambda_i} \alpha_i^2 = \frac{1}{\beta}\label{eq:SVD:spike2}\tag{43}\] For large \(p\), each \(\alpha_i^2\) concentrates around its mean value \(\mathbb{E}[\alpha_i^2]=\frac{1}{p}\) (again, this statement can be made rigorous), and 43 becomes \[\lim_{p\to \infty} \frac{1}{p}\sum_{i=1}^p \frac{\lambda_i}{\lambda-\lambda_i} = \frac{1}{\beta}\] Since the eigenvalues \(\lambda_1,\ldots,\lambda_p\) follow the Marčenko-Pastur distribution, the limit on the left hand side can be replaced by the integral \[\int_{\gamma_{-}}^{\gamma_{+}} \frac{t }{\lambda-t}\,d F_\gamma(t) = \frac{1}{\beta}\] Using an integral table (or an integral software), we find that \[\frac{1}{\beta} = \int_{\gamma_{-}}^{\gamma_{+}} \frac{t }{\lambda-t}\,d F_\gamma(t) =\frac{1}{4 \gamma} \left[2 \lambda -(\gamma_{-}+\gamma_{+})- 2 \sqrt{(\lambda-\gamma_{-})(\lambda-\gamma_{+})}\right]. \label{eq:SVD:Integral}\tag{44}\]

For \(\lambda=\gamma_{+}\), that is, when the top eigenvalue touches the right edge of the Marčenko-Pastur distribution, 44 becomes \(\frac{1}{4 \gamma} (\gamma_{+}-\gamma_{-})\). This is the critical point that one gets the pop out of the top eigenvalue from the bulk of the Marčenko-Pastur distribution. To calculate the critical value \(\beta_c\), we recall that \(\gamma_{-}=(1-\sqrt{\gamma})^2\) and \(\gamma_{+}=(1+\sqrt{\gamma})^2\), hence \[\frac{1}{\beta_c} = \frac{1}{4 \gamma} \left((1+\sqrt{\gamma})^2-(1-\sqrt{\gamma})^2\right).\] Therefore, the critical SNR is \[\label{eq:SVD:critical} \beta_{c}= \sqrt{\gamma}=\sqrt{\frac{p}{n}}.\tag{45}\] When \(\beta > \sqrt{\frac{p}{n}}\) one can observe the pop out of the top eigenvalue from the bulk.

Eq. (45 ) illustrates the interplay of the SNR \(\beta\), the number of samples \(n\), and the dimension \(p\). Low SNR, small sample size, and high dimensionality are all obstacles for detecting linear structure in noisy high dimensional data.

More generally, inverting the relationship between \(\beta\) and \(\lambda\) given by 44 (which simply amounts to solving a quadratic), we find that the largest eigenvalue \(\lambda\) of the sample covariance matrix \(S_n\) has the limiting value \[\label{eq:SVD:pop} \lambda \to \left\{\begin{array}{ccc} (\beta + 1)\left(1 + \frac{\gamma}{\beta} \right) & & \text{for }\;\; \beta \geq \sqrt{\gamma}, \\ & & \\ (1+\sqrt{\gamma})^2 & & \text{for }\;\; \beta < \sqrt{\gamma}. \end{array} \right.\tag{46}\] In the finite sample case \(\lambda\) will be fluctuating around that value.

Notice that the critical SNR value, \(\beta_{c} = \sqrt{\gamma}\) is buried deep inside the support of the Marčenko-Pastur distribution, because \(\sqrt{\gamma} < \gamma_{+} = (1+\sqrt{\gamma})^2\). In other words, the SNR does not have to be greater than the operator norm of the noise matrix in order for it to pop out. We see that the noise effectively pushes the eigenvalue to the right (indeed, \(\lambda > \beta\)).

The asymptotic squared correlation \(|\langle u, v\rangle |^2\) between the top eigenvector \(v\) of the sample covariance matrix and true signal vector \(u\) can be calculated in a similar fashion. The limiting correlation value turns out to be \[\label{eq:SVD:corr} |\langle v, u \rangle|^2 \to \left\{\begin{array}{ccc} \frac{1-\frac{\gamma}{\beta^2}}{1+\frac{\gamma}{\beta^2}} & & \text{for }\;\; \beta \geq \sqrt{\gamma} \\ & & \\ 0 & & \text{for }\;\; \beta < \sqrt{\gamma} \end{array} \right.\tag{47}\] Notice that the correlation value tends to 1 as \(\beta \to \infty\), but is strictly less than 1 for any finite SNR.

4.3.2 Wigner matrices↩︎

Another very important random matrix model is the Wigner matrix (and it will make appearances in Chapters 14, 15 and 13). Given an integer \(n\), a standard Gaussian Wigner matrix \(W\in\mathbb{R}^{n\times n}\) is a symmetric matrix with independent \(\mathcal{N}(0,1)\) off-diagonal entries (except for the fact that \(W_{ij}=W_{ji}\)) and jointly independent \(\mathcal{N}(0,2)\) diagonal entries. In the limit, the eigenvalues of \(\frac{1}{\sqrt{n}}W\) are distributed according to the so-called semi-circular law \[\label{eq:WignerSC} dSC(x) = \frac{1}{2\pi}\sqrt{4-x^2}1_{[-2,2]}(x)dx,\tag{48}\] and there is also a BBP like transition for this matrix ensemble [@Feral_Peche_BPPWigner]. More precisely, if \(v\) is a unit-norm vector in \(\mathbb{R}^n\) and \(\xi\geq 0\) then the largest eigenvalue of \(\frac{1}{\sqrt{n}}W+\xi vv^\top\) satisfies

• If \(\xi \leq 1\) then \[\lambda_{\max}\left(\frac{1}{\sqrt{n}}W+\xi vv^\top \right) \to 2,\]

• and if \(\xi > 1\) then \[\label{eq:1spike95Wigner:PCA1} \lambda_{\max}\left(\frac{1}{\sqrt{n}}W+\xi vv^\top \right) \to \xi + \frac{1}{\xi}.\tag{49}\]

The typical correlation, with \(v\), of the leading eigenvector \(v_{\max}\) of \(\frac{1}{\sqrt{n}}W+\xi vv^\top\) is also known:

• If \(\xi \leq 1\) then \[\left| \left\langle v_{\max},v \right\rangle \right|^2 \to 0,\]

• and if \(\xi > 1\) then \[\left| \left\langle v_{\max},v \right\rangle \right|^2 \to 1-\frac{1}{\xi^2}.\]

From a statistical viewpoint, a central question is to understand for different distributions of matrices, when it is possible to detect and estimate a spike in a random matrix [@AlexAmeliaAfonsoAnkur_PCA]. When the underlying random matrix corresponds to the adjacency matrix of a random graph and the spike to a bias on the distribution of the graph edges, corresponding to structural properties of the graph, the estimates above are able to predict important phase transitions in community detection in networks, as we will see in Chapter 13.

4.3.3 Rank and covariance estimation↩︎

The spike model and random matrix theory thus offer a principled way for determining the number of principal components, or equivalently of the rank of the hidden linear structure: simply count the number of eigenvalues to the right of the Marčenko-Pastur distribution. In practice, this approach for rank estimation is often too simplistic for several reasons. First, for actual datasets, \(n\) and \(p\) are finite, and one needs to take into account non-asymptotic corrections and finite sample fluctuations [@kritchman2008determining; @kritchman2009non]. Second, the noise may be heteroskedastic (that is, noise variance is different in different directions). Moreover, the noise statistics could also be unknown and it can be non-Gaussian [@liu2018pca]. In some situations it might be possible to estimate the noise statistics directly from the data and to homogenize the noise (a procedure sometimes known as “whitening”) [@leeb2018optimal]. These situations call for careful analysis, and many open problems remain in the field.

Another popular method for rank estimation is using permutation methods. In permutation methods, each column of the data matrix is randomly permutated, so that the low-rank linear structure in the data is destroyed through scrambling, while only the noise is preserved. The process can be repeated multiple times, and the statistics of the singular values of the scrambled data matrices are then used to determine the rank. In particular, only singular values of the original (unscrambled) data matrix that are larger than the largest singular value of the scrambled matrices (taking fluctuations into account of course) are considered as corresponding to signal and are counted towards the rank. The mathematical analysis of permutation methods is another active field of research [@dobriban2017permutation; @dobriban2018deterministic].

In some applications, the objective is to estimate the low rank covariance matrix of the clean signal \(\Sigma\) from the noisy measurements. We saw that in the spike model, the eigenvalues of the sample covariance matrix are inflated due to noise. It is therefore required to shrink the computed eigenvalues of \(S_n\) in order to obtain a better estimate of the eigenvalues of \(\Sigma\). That is, if \[S_n = \sum_{i=1}^p \lambda_i v_i v_i^\top\] is the spectral decomposition of \(S_n\), then we seek an estimator of \(\Sigma\), denoted \(\hat{\Sigma}\) of the form \[\hat{\Sigma} = \sum_{i=1}^p \eta(\lambda_i) v_i v_i^\top .\] The scalar nonlinearity \(\eta : \mathbb{R}^+ \to \mathbb{R}^+\) is known as the shrinkage function. An obvious shrinkage procedure is to estimate \(\beta = \eta(\lambda)\) from the computed \(\lambda\) by inverting (46 ) (and setting \(\beta=0\) for \(\lambda < \gamma_{+}\)). It turns out that this particular shrinker is optimal in terms of the operator norm loss. However, for other loss functions (such as the Frobenius norm loss), the optimal shrinkage function takes a different form [@donoho2018optimal]. The reason why the shrinker depends on the loss function is that the eigenvectors of \(S_n\) are not perfectly correlated with those of \(\Sigma\) but rather make some non-trivial angle, as in (47 ). In other words, the eigenvectors are noisy, and it may require more aggressive shrinkage to account for that error in the eigenvector. It can be shown that the eigenvector \(v\) of the sample covariance is uniformly distributed in a cone around \(u\) whose opening angle is given by (47 ). While we can improve the estimation of the eigenvalue via shrinkage, it is however unclear how to improve the estimation of the eigenvector (without any a priori knowledge about it). Finally, we remark that eigenvalue shrinkage also plays an important role in denoising.

Exercises↩︎

5 Linear Regression and Regularization↩︎

It is dawn at Mavericks, California. The swell forecast reads 10 meters at 18 seconds—a massive storm system has sent energy across thousands of miles of Pacific Ocean. For big wave surfers, these are exciting conditions! But they also come with all the dangers of a wipeout.

The hold-down—the time a surfer is trapped beneath the surface after a wipeout—is one of the most critical variables in a big-wave wipeout. Professional big wave surfers can hold their breath for four to five minutes under controlled conditions, but after the physical exertion of paddling into a 12-meter wave, that window shrinks dramatically. Add the disorientation of being tumbled by thousands of tons of water, and you have maybe roughly 20-30 seconds before panic sets in, and perhaps 60-90 seconds before blackout. This motivates a practical and (for some readers) urgent question: Given ocean conditions and wave characteristics, can we predict how long is a surfer likely to be held underwater?

Imagine we have collected data from 50 big wave sessions at various locations (Mavericks, Nazaré, Jaws, Teahupo’o,...) using underwater cameras and GPS trackers on professional surfers. For each wipeout, we have recorded hold-down time (seconds) \(y\), wave height (meter, measured face) \(x_1\), wave period (seconds between waves) \(x_2\), water depth at impact zone (meter) \(x_3\), and wind speed (knots) \(x_4\). Tomorrow, a massive northwest swell is forecast with wave height 15m, period 20s, depth 12m, and wind 18 kts. How long should we expect the hold-down to be?

While the true physics is nonlinear (involving fluid dynamics, buoyancy, and turbulent flow), a linear approximation may capture the first-order effects within the range of observed conditions. We thus postulate that hold-down time is approximately a linear function of the wave characteristics \[y = \theta_0 + \theta_1 x_1 + \theta_2 x_2 + \theta_3 x_3 + \theta_4 x_4 + \varepsilon,\] where \(\theta_0\) is the baseline hold-down time, \(\theta_1\) is the increase in hold-down per meter of wave height, \(\theta_2\) is the effect of wave period, \(\theta_3\) is the effect of water depth (deeper water means shorter hold-downs), \(\theta_4\) is the effect of wind, and \(\varepsilon\) represents random variation (wave shape, surfer position, etc.)

With our 50 observations, we can write the model in matrix form \[y = X \theta+ \varepsilon\] where \(y \in \mathbb{R}^{50}\) is the vector of observed hold-down times, \(X \in \mathbb{R}^{50 \times 5}\) contains the wave characteristics (with a column of 1’s for the intercept), and \(\theta\in \mathbb{R}^5\) are the unknown coefficients we wish to estimate. Once we have found \(\hat{\theta}\), we can predict tomorrow’s expected hold-down.

But how do we find the “best” estimate \(\hat{\theta}\)? Is our estimator optimal in some sense? How accurate are our predictions? Should we include all variables, or is simpler better? Linear regression and least squares will provide answers to these questions.

This chapter studies linear regression and least squares from four complementary viewpoints. After introducing the basic idea of linear regression, we first show that least squares is an orthogonal projection problem. Second, we show that it arises naturally from Gaussian noise and conditional-mean prediction. Third, we analyze its bias, variance, and prediction risk. Finally, we show how ridge regression stabilizes least squares by shrinking low-information directions, and we reinterpret this shrinkage as Wiener filtering.

5.1 Linear regression↩︎

Linear regression is the foundational supervised learning model in data science, statistics, and applied mathematics [@hastieelements; @LVanderberghe_SBoyd_book; @Golub_MatrixComputations]. While it is often motivated by statistical inference, its mathematical core is a geometric problem: finding the “best” projection of a target vector onto the subspace spanned by the predictors.

Let us formulate the linear regression problem in greater generality. Consider a dataset consisting of \(n\) observations and \(p\) features \(x_1,\dots,x_p\), represented by the matrix \(X \in \mathbb{R}^{n \times p}\) (often called the design matrix), and a vector of responses \(y \in \mathbb{R}^n\). Now we introduce a rather strong assumption on the underlying data model: we assume a linear relationship \(y = X\theta+ \varepsilon\), where, as in the previous example, \(\theta\in \mathbb{R}^p\) is the vector of unknown coefficients and \(\varepsilon\) represents noise. The goal is now to find the best estimate \(\hat{\theta}\) of \(\theta\).

Linear regression serves two primary purposes:

1. Prediction: Given a new predictor \(x_{\text{new}}\), we wish to predict the corresponding response \(y_{\text{new}}\). Our postulated linear model suggests the prediction \(\hat{y}_{\text{new}} = x_{\text{new}}^\top \hat{\theta}\) where \(\hat{\theta}\) is an estimate of \(\theta\).

2. Inference: We wish to understand the relationship between predictors and response, quantified by \(\theta\), such as: Which predictors are important? How does changing \(x_j\) affect \(y\)? What is the uncertainty in our estimates?

But how can we find an estimate \(\hat{\theta}\) in the first place? Because of noise, we do not expect to recover \(\theta\) exactly from finite data, nor do we expect \(y\)to lie exactly in the column space of \(X\). So, which error function (or loss function in the terminology of data science, statistics, and machine learning) shall we use to quantify a good (or even optimal) fit?

We propose to compute the \(\hat{\theta}\) that minimizes the sum of squared residuals \[\label{eq:l2reg} f(\theta) = \|y - X\theta\|^2.\tag{50}\]

This proposal immediately raises the question, why did we use the least-squares error in 50 and not some other measure of fit? After all, we could consider many alternatives, such as for instance the sum of absolute errors \(\sum_{i=1}^n |y_i -x_i^\top\theta|\), or the maximum absolute error \(\max_{i=1,\ldots,n} |y_i - x_i^\top\theta|\).

The choice of the squared loss 50 is not arbitrary. It emerges from multiple complementary perspectives: geometric interpretation, computational tractability, and statistical optimality. We examine each justification in turn¹³.

It is helpful to keep in mind the fundamental differences between linear regression and PCA (discussed in Chapter 4.2). While linear regression and PCA both involve fitting a linear structure to data, they have fundamentally different goals, and assumptions. Linear regression is a supervised method. The goal is to model the relationship between a set of predictor variables \(X\) and a response variable \(y\). It assumes the \(X\) values are known perfectly¹⁴ and all noise is in \(y\). In contrast, PCA is an unsupervised method that does not involve a response variable. Instead, it seeks a low-dimensional linear structure that captures as much variation in the given data \(X\) as possible. In a nutshell, regression asks: “How well can \(X\) explain \(y\)?” while PCA asks: “What is the most efficient way to represent \(X\)?”

5.2 Geometry, existence, and uniqueness of least squares↩︎

The least squares loss has some appealing geometric interpretations that are not shared by other loss functions, which we collect below (and the verification of which is left as exercises).

Let \(S \subset \mathbb{R}^n\) be a linear subspace. For every \(y \in \mathbb{R}^n\), there exists a unique decomposition \[\label{th:orthoproj} y = s + r, \quad \text{such that s \in S, and r \in S^\perp.}\tag{51}\] Furthemore, \(s\) is the unique minimizer of \[\min_{u \in S} \|y - u\|^2.\]

Let \(\mathop{\mathrm{Col}}(X)\) denote the column space of \(X\) and \(P_{\mathop{\mathrm{Col}}(X)}\) the orthogonal projection onto \(\mathop{\mathrm{Col}}(X)\). Applying the subspace decomposition above to the case \(S = \mathop{\mathrm{Col}}(X)\) implies that for any \(y \in \mathbb{R}^n\), there exists \(\hat{\theta}\in \mathbb{R}^p\) such that \[X\hat{\theta}= P_{\mathop{\mathrm{Col}}(X)}y.\] Moreover, we have \[\hat{\theta}\in \arg\min_{\theta\in \mathbb{R}^p} \|y - X\theta\|^2.\] In the hold-down example, the fitted vector \(X \hat{\theta}\) is the orthogonal projection of the observed hold-down times onto the subspace spanned by the measured wave features.

Let \(X \in \mathbb{R}^{n \times p}\) and \(y \in \mathbb{R}^n\). A vector \(\hat{\theta}\in \mathbb{R}^p\) minimizes \[\|y - X\theta\|^2\] if and only if \[X^\top (y - X\hat{\theta}) = 0.\] Equivalently, \[\label{eq:normalequations} X^\top X \hat{\theta}= X^\top y,\tag{52}\] which are also referred to as the normal equations.

Furthermore, if \(X\) has full column rank (\(\mathrm{rank}(X)=p\)), then \(X^T X\) is invertible and the normal equations have as unique solution the least squares estimator given by \[\label{eq:ols95standard} \hat{\theta}_{\text{LS}} = (X^\top X)^{-1}X^\top y.\tag{53}\]

However, this approach has both practical and theoretical limitations: (i) Computing \(X^\top X\) squares the condition number, amplifying rounding errors. This is not a big deal when the condition number of \(X\) is small (say, \(\kappa(X) = \mathcal{O}(1)\)), but becomes a major concern when \(X\) is not well conditioned. (ii) When \(\mathop{\mathrm{rank}}(X) < p\), the matrix \(X^\top X\) is singular and cannot be inverted; (iii) When \(p > n\), the system is underdetermined.

The SVD provides an elegant solution that avoids forming \(X^\top X\) explicitly and handles all cases (full rank, rank deficient, underdetermined) in a unified framework¹⁵.

The SVD solution has an appealing geometric interpretation: We project \(y\) onto each left singular vector \(u_i^\top y\), then scale by \(1/\sigma_i\), weight the result by the corresponding right singular vector \(v_i\), and sum over all positive singular values. This explicitly constructs the solution as a linear combination of the right singular vectors, with coefficients determined by the singular values and projections of \(y\).

For dense full-rank least-squares problems, QR factorization is usually the standard method because it is numerically stable and significantly cheaper than the SVD [@Golub_MatrixComputations]. The SVD is more robust, especially when the matrix is ill-conditioned or rank-deficient, because it reveals the numerical rank and naturally produces minimum-norm solutions.

5.3 Statistical properties of least squares↩︎

Before exploring the benefits of using the least squares error function from a statistical viewpoint, we will first impose the following standard assumptions on our linear model \(y= X\theta+\varepsilon\), which we will make throughout this section unless otherwise noted:

1. (A1) Exogeneity: \(\mathbb{E}[\varepsilon| X] = 0\).

2. (A2) Homoscedasticity: \(\mathop{\mathrm{Var}}(\varepsilon| X) = \sigma^2 I\).

3. (A3) Full rank: \(\mathop{\mathrm{rank}}(X) = p\) (assumes \(n \geq p\)).

It is instructive to explain the role of each of these assumptions.

Exogeneity is the most important statistical assumption. It stipulates that the errors are uncorrelated with the predictors. Exogeneity would mean, for instance, that the unmodeled part of the hold-down time in our surfing example is not systematically correlated with the recorded wave height, period, depth, or wind speed. It also means there is no “omitted variable bias”, i.e., we have not left out important confounders. Thus, exogeneity is required for unbiasedness and consistency of the least squares solution, as well as for a causal interpretation of coefficients¹⁶.

Homoscedasticity means the errors have constant variance, i.e., the “variability” of the errors around the regression line is the same for all values of the predictors¹⁷. Without homoscedasticity, the least squares estimator may still be unbiased but can be inefficient and produce misleading standard errors.

We already encountered the full rank assumption in the previous section, and have seen that it ensures invertibility of \(X^T X\) and uniqueness of the solution. While non-invertibility is not a “deal-breaker”, enforcing uniqueness of the solution would require additional assumptions and regularization (for instance we can replace the inverse with the Moore-Penrose pseudo-inverse, as previously shown), which we will discuss in more detail in Chapter 5.4.

The most compelling statistical justification of the least squares error function 50 arguably comes from maximum likelihood theory.

Let us recall the concept of the likelihood function. Assume we have a set of i.i.d.observations \(\mathcal{D} = \{x_1, x_2, \dots, x_n\}\) drawn from a probability density (or mass) function \(f(x | \theta)\). The joint probability of the entire dataset is the product of the individual probabilities \(\mathbb{P}(\mathcal{D} | \theta) = \prod_{i=1}^n f(x_i | \theta).\) When we treat the data \(\mathcal{D}\) as fixed and the parameter \(\theta\) as a variable, this function is called the likelihood function, denoted as \(\mathcal{L}(\theta)\), \[\mathcal{L}(\theta; x_1, \dots, x_n) = \prod_{i=1}^n f(x_i | \theta).\] The Maximum Likelihood Estimator (MLE) \(\hat{\theta}_{MLE}\) is defined as the argument that maximizes the likelihood function, i.e., \[\hat{\theta}_{MLE} := \operatornamewithlimits{argmax}_{\theta\in \Theta} \mathcal{L}(\theta),\] where \(\Theta\) denotes a parameter space.

Under suitable regularity conditions (identifiability, correct model specification, and mild technical assumptions), one can show (using e.g.the Uniform Law of Large Numbers and the Argmax Theorem) that the maximum likelihood estimator is consistent [@wald1949note; @van2000asymptotic], i.e. \[\hat{\theta}_n \xrightarrow{P} \theta_0 \quad \text{as } n \to \infty,\] where \(\theta_0\) is the true parameter and \(\hat{\theta}_n\) is the MLE based on \(n\) samples.

Notice that to actually compute the MLE numerically, one usually considers the log-likelihood function \(\ell(\theta) = \log \mathcal{L}(\theta)\). Since the natural logarithm is a monotonically increasing function, the value of \(\theta\) that maximizes \(\mathcal{L}(\theta)\) also maximizes the log-likelihood \[\ell(\theta) = \sum_{i=1}^n \log f(x_i | \theta),\] which has numerically much more stable properties. The score is the gradient of the log-likelihood function with respect to \(\theta\). To find the MLE, we typically solve the score equation by setting the gradient of the log-likelihood to zero, \[\nabla_\theta\ell(\theta) = 0.\]

Except for a few cases (like least squares, see also Theorem 5.2 below) the MLE typically lacks a closed-form solution. While many standard models (like logistic regression, see Chapter 10.2) remain computationally tractable due to convexity, more complex models result in non-convex optimization problems where finding the global MLE is often computationally intractable.

When we add a distributional assumption on the errors, we can show that the least squares estimator is actually the maximum likelihood estimator.

Note that the MLE 57 is biased; the unbiased estimator uses \(n-p\) instead of \(n\) in the denominator (see also Exercise 4.21).

Theorem 5.3 tells us that if we believe \(\mathbb{E}[Y | x] = x^\top\theta\), then minimizing the squared error is the natural way to estimate \(\theta\); we are finding the linear function that best approximates the conditional mean. We leave it to the reader to verify that using instead the sum of absolute errors estimates the conditional median.

A natural question is whether the least squares estimator is optimal among all unbiased linear estimators?

5.4 Bias, variance, and risk↩︎

A fundamental tension in statistical estimation is the bias-variance tradeoff: reducing variance often increases bias, and vice versa. In fact, the bias–variance tradeoff is a fundamental identity of squared-error prediction and randomness in estimation, and not specific to linear models. It arises from basic properties of expectations and variance.

Also recall that \(\mathop{\mathrm{Bias}}(\hat{\theta}) = \mathbb{E}[\hat{\theta}] - \theta\) and \(\operatorname{Cov}(\hat{\theta}) = \mathbb{E}\left[ (\hat{\theta} - \mathbb{E}[\hat{\theta}])(\hat{\theta} - \mathbb{E}[\hat{\theta}])^T \right]\).

We point out that Theorem 5.6 does allow for heteroscedasticity and for correlated errors, Moreover, it does not impose any distributional assumptions beyond zero-mean and finite second moments. Thus, the decomposition holds for any estimator under very mild conditions. This generality makes it a fundamental result in statistical learning theory.

While Theorem 5.6 concerns estimating the parameter vector \(\theta\), in practice we often care about predicting new observations. The bias-variance decomposition for prediction includes an additional irreducible error term.

The decomposition 59 in Corollary 5.1 reveals three sources of prediction error:

1. Bias: Systematic error from using a biased estimator or misspecified model. For least squares, this is zero since \(\mathbb{E}[\hat{\theta}_{\text{LS}}] = \theta\). For ridge regression (discussed below), this is nonzero but can be outweighed by variance reduction.

2. Variance: Random error from the particular training sample. Predictions vary across different training datasets. This term is larger when \(x_{\text{new}}\) is far from the training data (extrapolation), the design matrix \(X\) is ill-conditioned, or the sample size \(n\) is small relative to dimension \(p\).

3. Irreducible error: Fundamental noise in the response variable that cannot be eliminated by any estimator. This represents the intrinsic randomness in \(y\) that cannot be explained by \(x\), even with infinite data.

The bias-variance tradeoff refers to the tension between the first two terms. This brings up the question whether we can accept some bias to achieve substantial variance reduction to achieve lower total prediction error? This is indeed possible and can be accomplished via proper regularization.

For example, for the case of linear regression, when the matrix \(X\) is ill-conditioned, the least squares solution remains unbiased but can have extremely high variance. Historically, the focus was on unbiasedness, because datasets often had rather small values of \(p\) compared to \(n\), making ill-conditionedness of \(X\) a much less common scenario. The ill-conditionedness of \(X\) happens typically more often when the number of variables \(p\) is comparable to the sample size \(n\) (here, we still assume \(p \le n\))¹⁸. In this regime, minimizing prediction error rather than preserving unbiasedness becomes central, and regularization, such as ridge regression, improves performance by trading bias for variance reduction. We will explore this topic in Chapter 5.5.

Corollary 5.1 also forms the basis for constructing prediction intervals, see e.g. [@seber2003linear; @hastieelements] for details. Specializing to the case of a linear model \(y = X\theta + \varepsilon\), if the model is correctly specified and the estimator is unbiased, then under Gaussian errors a \((1-\alpha)\) prediction interval for \(y_{\text{new}}\) is \[\label{eq:predictioninterval} \hat{y}_{\text{new}} \pm t_{1-\alpha/2, n-p} \cdot \sqrt{\hat{\sigma}^2 + \hat{\sigma}^2 \, x_{\text{new}}^\top(X^\top X)^{-1}x_{\text{new}}},\tag{61}\] where \(\hat{\sigma}^2 = \frac{1}{n-p}\|y - X\hat{\theta}\|^2\) and \(t_{1-\alpha/2, n-p}\) is the critical value of the \(t\)-distribution. Note that the first term under the square root in 61 accounts for the irreducible error and the second term for parameter uncertainty. We emphasize that the probability is in the procedure, not in a single realization. Namely, if we repeatedly collect a training dataset, fit a model (with the same method), construct a prediction interval, and observe a new data point, then on average, \((1-\alpha) \times 100\%\) of these intervals will contain the new observation. We emphasize that for a biased estimator (such as ridge regression or lasso) constructing prediction intervals is more involved [@efron1994introduction; @hastieelements].

Generalization gap. In statistical learning theory, the bias-variance tradeoff tells us how the components of risk behave as model complexity changes. Let us recall the notion of risk, which represents the average cost incurred by a predictive model when applied to the entire population of possible data points. Formally, we define risk as follows:

Using the definition of expectation, the risk can be written as the integral over the joint distribution \[R(f) = \int_{\mathcal{X} \times \mathcal{Y}} L(y, f(x)) \, dP(x, y).\]

In data science, the distribution \(P\) is almost always unknown. This leads to the distinction that is central to learning theory:

• Theoretical Risk (\(R(f)\)): This is the risk we want to minimize but cannot calculate in practice (also known as True Risk or Population Risk).

• Empirical Risk (\(\hat{R}_n(f)\)): This is what we actually minimize using a training set of size \(n\); it is given by \[\hat{R}_n(f) = \frac{1}{n} \sum_{i=1}^n L(y_i, f(x_i)).\]

The generalization gap is the difference between these two, \(|R(f) - \hat{R}_n(f)|\). An important part of theory in data science and machine learning is dedicated to bounding this gap to ensure that a model that performs well on training data will also perform well “in the wild”.

Sometimes it is stated in the literature that the generalization gap is the difference between test error and training error. However, the test error is only an estimator of the true risk. Hence, while the test error minus training error is a standard proxy, the true generalization gap is defined against the population risk, and conflating the two can obscure the effects of test set noise, data reuse, and distribution shift.

Model mismatch, underfitting, overfitting: In reality, the prediction error decomposition 59 is almost always wrong, since the true data-generating process does rarely perfectly fit the chosen model. For instance in case of linear regression, when the underlying data is generated by a process that is more complex than a linear model (as is the case with the surfing/hold-down time example used in the beginning of this chapter), we encounter model mismatch, also known as specification error or structural bias.

Model mismatch is closely related to the issue of underfitting/overfitting, but they are not identical. Model mismatch is a design-time issue: the chosen model class does not contain a function capable of representing the true relationship in the data. Even with infinite data, the best possible predictor within the model class cannot recover the true signal.

Underfitting, in contrast, is a run-time phenomenon: the fitted model fails to capture the structure of the data. This may occur because the model class is too restrictive (model mismatch), for example fitting a line to a quadratic relationship. But it can also arise when the estimation procedure overly constrains the parameters, for example through strong regularization or insufficient training. Underfitting typically corresponds to large approximation bias. Overfitting occurs when the model class is too flexible, and thus the estimator fits noise in finite samples. Examples include using a high-degree polynomial or linear regression when \(p > n\) without regularization (see Chapter 5.5). Overfitting manifests as large variance.

5.5 Regularization and ridge regression↩︎

Assume we want to predict the optimal drug dosage for a patient using the expression levels of thousands of genes. In such problems the number of features (genes) can be far larger than the number of patients in the study, leading to a high-dimensional regression problem with \(p \gg n\).

This problem falls into the area of high-dimensional linear regression (assuming one considers the underlying model to be linear). More generally, high-dimensional linear regression refers broadly to regimes where \(p\) is comparable to or larger than \(n\). This includes the cases \(p > n\), \(p \gg n\), and \(p/n \to \gamma > 0\).

Standard least squares works well when the number of predictors \(p\) is small relative to the sample size \(n\) and the design matrix \(X\) is well-conditioned. However, when predictors are highly correlated or \(p\) approaches \(n\), least squares becomes unstable, since small changes in the data can lead to large changes in coefficient estimates. In the case \(p > n\), \(X^\top X\) is singular, and the minimizer is not unique.

Ridge regression¹⁹ addresses these problems by adding a penalty term that shrinks coefficients toward zero, trading a small amount of bias for a potentially substantial variance reduction [@hoerl1970ridge; @hastieelements]. In numerical analysis and optimization ridge regression is also known as Tikhonov regularization (named for Andrey Tikhonov [@tikhonov1963solution]).

The ridge objective balances two competing goals:

• Fit the data: Minimize \(\|y - X\theta\|^2\),

• Keep coefficients small: Minimize \(\|\theta\|^2\) (regularization).

The regularization parameter \(\lambda\) controls the trade-off between these two goals.

Unlike many regularized methods, ridge regression has an explicit solution.

When \(X\) has full column rank, the least squares solution is given by \(\hat{\theta}_{\text{LS}} = (X^\top X)^{-1}X^\top y\). Ridge regression modifies this by adding \(\lambda I_p\) to \(X^\top X\), which stabilizes the inversion when eigenvalues of \(X^\top X\) are small and ensures invertibility even when \(X^\top X\) is singular.

The SVD provides deep insight into how ridge regression works.

It is instructive to compare 66 with the least squares solution (when \(r = p\)): \[\label{eq:LS95svd} \hat{\theta}_{\text{LS}} = \sum_{i=1}^p \frac{1}{\sigma_i}(u_i^\top y)v_i.\tag{67}\] In the surfing exmaple, if wave height and wave period are highly correlated in the dataset, ridge regression prevents the fitted coefficients from becoming unstable by shrinking poorly determined combinations of features.

So, how does ridge regression fare with respect to the bias-variance tradeoff? The following result is straightforward (the proof is left to the reader as Exercise 5.4).

Unlike in Theorem 5.4, we do not require that \(X\) has full rank (assumption (A3)), since being able to handle the rank-deficient case is exactly a key motivation for introducing ridge regression in the first place.

Choosing the regularization parameter: Proposition 5.9 suggests that by choosing an appropriate value for the regularization parameter \(\lambda\) we can achieve a good balance between overfitting and underfitting of the solution. See Figure 14 for an illustration of this bias-variance tradeoff for the prediction error and the interplay between regularization and over/underfitting. The optimal choice for \(\lambda\) via the minimum of the MSE, indicated in Figure 14, is unfortunately not computable in practice, since it would require access to the true solution \(\theta\). This begs the question how to find \(\lambda\) in practice.

Cross-validation, a widely used technique in data science and machine learning for hyperparameter selection, is a standard way to determine a good choice for \(\lambda\). We describe the cross-validation approach in the context of ridge regression, but it should be understood that this technique is more broadly applicable.

If we know that the noise is Gaussian and we have a good estimate for the noise variance \(\sigma^2\), instead of the rather expensive cross validation we can use Stein’s unbiased risk estimate (SURE). SURE provides an analytic estimate of prediction risk under Gaussian noise [@stein1981estimation] and can therefore be used to select hyperparameters such as the ridge penalty or the lasso regularization parameter (see 68 ) by minimizing the SURE value. While SURE often has lower variance than cross-validation and uses the full dataset, it relies on stronger assumptions—besides assuming Gaussian noise and knowledge of the noise variance, it imposes a certain smoothness condition (weak differentiability) on the estimator. Cross-validation, though computationally heavier, is thus more broadly applicable.

The work of Donoho and Johnstone [@donoho1995adapting] was instrumental in moving SURE from a theoretical result in statistics into a cornerstone of signal and image processing, in particular for data denoising (see also [@blu2007sure]).

Stein’s SURE has its roots in his earlier work [@stein1956inadmissibility; @james1961estimation]. In fact, while the Gauss-Markov theorem states that the least squares estimator is the best linear unbiased estimator, Charles Stein proved in a seminal paper [@stein1956inadmissibility] the startling result that in dimensions \(p\ge 3\), the least squares estimator is “inadmissible”. This result, also known as Stein’s Paradox, means there exists a biased “shrinkage” estimator—the James-Stein estimator [@james1961estimation]—that achieves a lower MSE for any true parameter value by pulling individual estimates toward a common average or the origin [@gruber2017improving].

Linear regression and sparsity: In many applications in high-dimensional settings—such as in financial econometrics and risk management, natural language processing, recommender systems, neuroimaging, or cheminformatics—it is often believed that only a small subset of variables truly influence the process we try to model. Consequently, it is desirable to construct estimators that set many coefficients exactly to zero. Such sparse models can improve predictive performance by reducing variance, and they also lead to models that are much easier to interpret.

This motivation leads to the lasso (least absolute shrinkage and selection operator) [@Tib96; @hastieelements], which replaces the \(\ell_2\)-penalty with an \(\ell_1\)-penalty, namely we consider the regularized regression problem \[\label{eq:lasso1} \min_{\theta} \, \, \|y - X\theta\|^2 + \lambda \|\theta\|_1,\tag{68}\] for some \(\lambda > 0\). Unlike ridge regression, the \(\ell_1\)-penalty encourages coefficients to be exactly zero, thereby performing variable selection as part of the estimation procedure. We will discuss the reason why the \(\ell_1\)-penalty promotes sparsity, associated theoretical guarantees, and numerical considerations extensively in Chapter 16 in connection with compressive sensing.

Outliers: An important consideration in linear regression (or in any data science task for that matter) is how different error functions handle outliers, i.e., data points that deviate significantly from the pattern exhibited by the majority of observations. The least squares solution is sensitive to outliers, since a single observation with large error \((y_i - x_i^\top\theta)\) contributes quadratically to the objective, potentially dominating the fit. In contrast, the solution of the \(\ell_1\) objective (also called \(\ell_1\)-regression or least absolute deviations regression) given by \[\min_{\theta} \, \, \|y - X\theta\|_1,\] is more robust, since outliers contribute only linearly. However, one downside of the \(\ell_1\) loss function is its lack of smoothness, since the function \(f(x)=|x|\) is not differentiable at the origin. Also, \(\ell_1\) estimators are typically less statistically efficient compared to \(\ell_2\) estimators when the noise is Gaussian.

A good compromise is provided by the Huber method [@huber1992robust]. The Huber loss is defined as \[L_\delta(r_i) = \begin{cases} \frac{1}{2}r_i^2 & \text{if } |r_i| \leq \delta, \\ \delta |r_i| - \frac{1}{2}\delta^2 & \text{if } |r_i| > \delta, \end{cases}\] where \(r_i = y_i - x_i^\top \theta\) is the \(i\)-th residual and \(\delta > 0\) is a tuning parameter (see also Figure 15 (b)). Because the derivative is continuous everywhere (including the transition points at \(\pm \delta\)), the Huber loss is \(C^1\)-continuous. This ensures that gradient-based optimization algorithms (like gradient descent, see Chapter 9.4) converge reliably.

Figure 15 shows a comparison of the \(\ell_2\)-loss, the \(\ell_1\)-loss, and the Huber loss for a linear regression example with a few outliers (left panel), together with the associated loss functions (right panel). It is evident that the solution of the \(\ell_2\) objective is affected most by the few outliers. The solutions associated with \(\ell_1\)-loss and with Huber loss, respectively (which are almost identical in this case) are only marginally affected by the outliers.

5.6 Wiener filtering↩︎

In classical linear regression, we often treat the data as a collection of fixed points to be “fit.” However, in the presence of noise, a more powerful perspective is to treat the measurement as a corrupted signal. This leads us to the Wiener Filter, the mathematically optimal way to extract a signal from noise via a linear estimator by balancing experimental evidence against prior knowledge [@wiener1949extrapolation; @kailath2000linear].

In this section we will show that ridge regression is simultaneously a maximum a posteriori (MAP) estimator under a Gaussian prior, the solution to a properly normalized loss function whose two terms have direct statistical interpretations, and a componentwise Wiener filter applied in the SVD basis of the design matrix. The unifying object is the signal-to-noise ratio (SNR), which controls how strongly the estimator trusts the data relative to the prior. We will see that ridge regression acts as a component-wise Wiener filter in the SVD basis.

The Wiener filter addresses a fundamental question: Given a noisy measurement \(y = s + n\), how much should we trust \(y\) to estimate the true signal \(s\)? Let us first gain some intuition for the scalar case. If we assume the signal \(s\) and noise \(n\) are independent, zero-mean scalar random variables with finite variances \(\sigma_s^2\) and \(\sigma_n^2\), the optimal linear estimate \(\hat{s} = Gy\) is found by minimizing the MSE \[\label{eq:wiener1} \min_{G} \, \mathbb{E}[(s - Gy)^2].\tag{69}\]

To find the minimum mean squared error (MMSE) estimator, we expand the cost function 69 , use the fact the signal and the noise are uncorrelated, and set the derivative of 69 with respect to \(G\) to zero. We arrive at the solution, also called Wiener Gain, \[G = \frac{\sigma_s^2}{\sigma_s^2 + \sigma_n^2} = \frac{1}{1 + \frac{1}{\text{SNR}}},\] where \(\text{SNR} = \sigma_s^2 / \sigma_n^2\) represents the signal-to-noise ratio²⁰. For high SNR, i.e., if the signal is much stronger than the noise, \(G \approx 1\), we trust the measurement. For low SNR, i.e., if the noise dominates, \(G \approx 0\), we ignore the measurement to avoid a high-variance error, relying instead on our prior knowledge (that the signal is centered at zero).

We can now generalize this logic to the multivariate case. Consider the linear model \[\label{eq:linear95model} y = X\theta + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, \sigma^2 I_n),\tag{70}\] where \(y \in \mathbb{R}^n\) is the measurement vector, \(X \in \mathbb{R}^{n \times p}\) is the design matrix, \(\theta \in \mathbb{R}^p\) is the unknown parameter vector, and \(\varepsilon\) is Gaussian noise with known variance \(\sigma^2 > 0\). The likelihood of \(\theta\) given the observations is therefore \[\label{eq:wiener2} p(y \mid \theta) = (2\pi\sigma^2)^{-n/2} \exp \Big(-\frac{1}{2\sigma^2}\|y - X\theta\|^2\Big).\tag{71}\] We place a zero-mean isotropic Gaussian prior on \(\theta\): \[\label{eq:prior} \theta \sim \mathcal{N}(0, \tau^2 I_p),\tag{72}\] where \(\tau^2 > 0\) is the prior variance.

By Bayes’ theorem the posterior is \(p(\theta | y) \propto p(y | \theta)\,p(\theta)\). Taking the negative log-posterior and dropping terms that do not depend on \(\theta\), the MAP estimate is the minimizer of \[\label{eq:normalised95loss} \mathcal{L}(\theta) = \frac{1}{2\sigma^2}\|y - X\theta\|^2 + \frac{1}{2\tau^2}\|\theta\|^2.\tag{73}\] This is the normalized loss function. The first term is the data-fidelity term: it measures how well \(\theta\) explains the observations, weighted by the inverse noise variance \(1/\sigma^2\). The second term is the regularization term: it penalizes large coefficients, weighted by the inverse prior variance \(1/\tau^2\). Both terms now have identical units (nats of log-probability), making the trade-off between them dimensionally consistent.

Consider 73 and set \(\nabla_\theta \mathcal{L}(\theta) = 0\); this gives \[\label{eq:gradient} -\frac{1}{\sigma^2}X^\top(y - X\theta) + \frac{1}{\tau^2}\theta = 0.\tag{75}\] Rearranging terms and multiplying by \(\sigma^2\) yields \[\label{eq:ridge95normal1} \left(X^\top X + \lambda I_p\right)\hat{\theta} = X^\top y, \qquad \lambda = \frac{\sigma^2}{\tau^2}.\tag{76}\] Since \(X^\top X \succeq 0\) and \(\lambda > 0\), the matrix \(X^\top X + \lambda I_p \succ 0\) is strictly positive definite for every \(\lambda > 0\), regardless of the rank of \(X\). The unique MAP (and ridge) estimate is therefore \[\label{eq:ridge95estimate} \hat{\theta}_{\mathrm{MAP}} = \hat{\theta}_{\mathrm{ridge}} = \left(X^\top X + \lambda I_p\right)^{-1}X^\top y.\tag{77}\] We already confirmed earlier that this gives a strict global minimum.

Because both the likelihood 71 and the prior 72 are Gaussian, the posterior is also Gaussian: \[\label{eq:posterior} \theta \mid y \;\sim\; \mathcal{N}\!\left( \hat{\theta}_{\mathrm{MAP}},\; \sigma^2\left(X^\top X + \lambda I_p\right)^{-1} \right).\tag{78}\]

It is instructive to revisit the componentwise structure of the ridge estimator by using the SVD of the design matrix, \(X=U \Sigma V^\top\), which we derived in 66 . The ridge estimate takes the form \[\label{eq:ridge95svd1} \hat{\theta}_{\mathrm{ridge}} = V(\Sigma^2 + \lambda I_p)^{-1}\Sigma U^\top y = \sum_{i=1}^p \frac{\sigma_i}{\sigma_i^2 + \lambda} (u_i^\top y)\, v_i.\tag{79}\] Now we apply the rotations \[\label{eq:rotated95coords} \tilde{\theta} = V^\top\theta \in \mathbb{R}^p, \qquad \tilde{y} = U^\top y \in \mathbb{R}^p.\tag{80}\] In these coordinates the model 70 becomes \(\tilde{y}_i = \sigma_i \tilde{\theta}_i + \tilde{\varepsilon}_i\) with \(\tilde{\varepsilon}_i \sim \mathcal{N}(0, \sigma^2)\) i.i.d., and the ridge estimate 79 decouples into \(p\) independent scalar problems: \[\label{eq:ridge95component} \hat{\tilde{\theta}}_i = \frac{\sigma_i}{\sigma_i^2 + \lambda}\,\tilde{y}_i = \frac{1}{\sigma_i}\cdot\frac{\sigma_i^2}{\sigma_i^2 + \lambda}\,\tilde{y}_i, \qquad i = 1, \ldots, p.\tag{81}\]

We now introduce the per-component signal-to-noise ratio. Under the prior 72 , the \(i\)-th rotated coefficient satisfies \(\tilde{\theta}_i \sim \mathcal{N}(0, \tau^2)\), so the signal component \(\sigma_i\tilde{\theta}_i\) has variance \(\sigma_i^2\tau^2\). The corresponding observation \(\tilde{y}_i = \sigma_i\tilde{\theta}_i + \tilde{\varepsilon}_i\) has noise variance \(\sigma^2\). The SNR in the \(i\)-th direction is therefore \[\label{eq:snr95def} \mathrm{SNR}_i = \frac{\sigma_i^2\,\tau^2}{\sigma^2}.\tag{82}\] Substituting \(\lambda = \sigma^2/\tau^2\), note that \[\label{eq:lambda95snr} \frac{\lambda}{\sigma_i^2} = \frac{\sigma^2/\tau^2}{\sigma_i^2} = \frac{1}{\mathrm{SNR}_i},\tag{83}\] so the shrinkage factor in 81 can be written as \[\label{eq:shrinkage95snr} \frac{\sigma_i^2}{\sigma_i^2 + \lambda} = \frac{1}{1 + \lambda/\sigma_i^2} = \frac{1}{1 + 1/\mathrm{SNR}_i} = \frac{\mathrm{SNR}_i}{1 + \mathrm{SNR}_i}.\tag{84}\] The ridge estimate in the principal-component basis therefore takes the form \[\label{eq:wiener95form} \hat{\tilde{\theta}}_i = \underbrace{\frac{1}{\sigma_i}}_{\text{LS inversion}} \cdot \underbrace{\frac{1}{1 + 1/\mathrm{SNR}_i}}_{\text{Wiener filter gain}} \cdot\; \tilde{y}_i, \qquad i = 1, \ldots, p.\tag{85}\]

The connection to Wiener filtering deserves elaboration. Consider the scalar model \(\tilde{y}_i = \sigma_i\tilde{\theta}_i + \tilde{\varepsilon}_i\) with \(\tilde{\theta}_i \sim \mathcal{N}(0,\tau^2)\) and \(\tilde{\varepsilon}_i \sim \mathcal{N}(0,\sigma^2)\) independent. This is exactly the signal-plus-noise model underlying the Wiener filter that we introduced at the beginning of this section: we observe the signal \(\sigma_i\tilde{\theta}_i\) corrupted by noise \(\tilde{\varepsilon}_i\), and we wish to estimate \(\tilde{\theta}_i\).

The linear MMSE estimator of \(\tilde{\theta}_i\) given \(\tilde{y}_i\) is: \[\label{eq:mmse95component} \mathbb{E}[\tilde{\theta}_i \mid \tilde{y}_i] = \frac{\mathrm{Cov}(\tilde{\theta}_i,\tilde{y}_i)}{\mathrm{Var}(\tilde{y}_i)} \,\tilde{y}_i = \frac{\sigma_i\tau^2}{\sigma_i^2\tau^2 + \sigma^2} \,\tilde{y}_i = \frac{1}{\sigma_i}\cdot \frac{\sigma_i^2\tau^2}{\sigma_i^2\tau^2 + \sigma^2} \,\tilde{y}_i.\tag{86}\] Dividing numerator and denominator of the second term on the right-hand side of 86 by \(\sigma^2\) gives \[\label{eq:mmse95snr} \mathbb{E}[\tilde{\theta}_i \mid \tilde{y}_i] = \frac{1}{\sigma_i}\cdot\frac{\mathrm{SNR}_i}{1 + \mathrm{SNR}_i} \,\tilde{y}_i = \frac{1}{\sigma_i}\cdot\frac{1}{1 + 1/\mathrm{SNR}_i} \,\tilde{y}_i,\tag{87}\] which is identical to 85 . Since the model is jointly Gaussian, the linear MMSE estimator coincides with the MAP estimator, confirming that ridge regression achieves the minimum MSE in each principal direction simultaneously.

The Wiener filter gain \(W_i = \mathrm{SNR}_i/(1+\mathrm{SNR}_i)\) interpolates between two extreme regimes of high SNR and low SNR by continuously interpolating between the data and the prior. The fraction of the least squares estimate that is retained is exactly \(\mathrm{SNR}_i/(1 + \mathrm{SNR}_i) \in (0,1)\).

Reflecting on the derivations above, we have shown that the following three problems have identical solutions:

1. Regularized least squares (ridge regression). Minimize the standard ridge objective: \[\hat{\theta}_{\mathrm{ridge}} = \operatorname{argmin}_\theta\; \|y - X\theta\|^2 + \lambda\|\theta\|^2.\]

2. MAP estimation under a Gaussian prior. Minimize the normalized probabilistic loss: \[\hat{\theta}_{\mathrm{MAP}} = \operatorname{argmin}_\theta\; \frac{1}{2\sigma^2}\|y - X\theta\|^2 + \frac{1}{2\tau^2}\|\theta\|^2, \qquad \lambda = \frac{\sigma^2}{\tau^2}.\]

3. Componentwise Wiener filtering in the SVD basis. Apply the Wiener filter gain to each singular value component: \[\hat{\tilde{\theta}}_i = \frac{1}{\sigma_i}\cdot\frac{1}{1 + 1/\mathrm{SNR}_i}\cdot\tilde{y}_i, \qquad \mathrm{SNR}_i = \frac{\sigma_i^2\tau^2}{\sigma^2} = \frac{\sigma_i^2}{\lambda}.\]

Exercises↩︎

rng(100);
n = 500; p = 1000; X = randn(n, p);
theta_true = [ones(10, 1); zeros(p-10, 1)]; 
y = X * theta_true + 0.5 * randn(n, 1); 

6 Graphs, Networks, and Clustering↩︎

A crucial part of data science consists of the study of networks. Network science, or graph theory, unifies the study of diverse types of networks, such as social networks, protein-protein interaction networks, gene-regulation networks, and the internet. In this chapter we introduce graph theory and treat the problem of clustering, to identify similar data points, or vertices, in (network) data.

6.1 PageRank↩︎

Before we introduce the formalism of graph theory, we describe the celebrated PageRank algorithm. This algorithm is a principal component²¹ behind web search algorithms, in particular in Google. The goal of PageRank is to quantitatively rate the importance of each page on the web, allowing the search algorithm to rank the pages and thereby present to the user the more important pages first. Search engines such as Google have to carry out three basic steps:²²

• Crawl the web and locate all, or as many as possible, accessible webpages.

• Index the data of the webpages from step 1, so that they can be searched efficiently for relevant key words or phrases.

• Rate the importance of each page in the database, so that when a user does a search and the subset of pages in the database with the desired information has been found, the more important pages can be presented first.

Here, we will focus on the third step. We follow mainly the derivation in [@Bryan2006]. We aim to develop a score of importance for each webpage. A score will be a non-negative number. A key idea in assigning a score to any given webpage is that the page’s score is derived from the links made to that page from other webpages — “A person is important not if it knows a lot of people, but if a lot of people know that person”.

Suppose the web of interest contains \(n\) pages, each page indexed by an integer \(k\), \(1 \le k \le n\). A typical example is illustrated in Figure 16, in which an arrow from page \(k\) to page \(j\) indicates a link from page \(k\) to page \(j\). Such a web is an example of a directed graph. The links to a given page are called the backlinks for that page. We will use \(x_k\) to denote the importance score of page \(k\) in the web. \(x_k\) is nonnegative and \(x_j > x_k\) indicates that page \(j\) is more important than page \(k\).

A very simple approach is to take \(x_k\) as the number of backlinks for page \(k\). In the example in Figure 16, we have \(x_1 = 2, x_2 = 1, x_3 = 3\), and \(x_4 = 2\), so that page 3 is the most important, pages 1 and 4 tie for second, and page 2 is least important. A link to page \(k\) becomes a vote for page \(k\)’s importance. This approach ignores an important feature one would expect a ranking algorithm to have, namely, that a link to page \(k\) from an important page should boost page \(k\)’s importance score more than a link from an unimportant page. In the web of Figure 16, pages 1 and 4 both have two backlinks: each links to the other, but the second backlink from page 1 is from the seemingly important page 3, while the second backlink for page 4 is from the relatively unimportant page 2. As such, perhaps the algorithm should rate the importance of page 1 higher than that of page 4.

As a first attempt at incorporating this idea, let us compute the score of page \(j\) as the sum of the scores of all pages linking to page \(j\). For example, consider the web in our toy example. The score of page 1 would be determined by the relation \(x_1 = x_3 + x_4\). However, since \(x_3\) and \(x_4\) will depend on \(x_1\), this seems like a circular definition, since it is self-referential (it is exactly this self-referential property that will establish a connection to eigenvector problems!).

We also seek a scheme in which a webpage does not gain extra influence simply by linking to lots of other pages. We can do this by reducing the impact of each link, as more and more outgoing links are added to a webpage. If page \(j\) contains \(n_j\) links, one of which links to page \(k\), then we will boost page \(k\)’s score by \(x_j/n_j\), rather than by \(x_j\). In this scheme, each webpage gets a total of one vote, weighted by that web page’s score, that is evenly divided up among all of its outgoing links. To quantify this for a web of \(n\) pages, let \(L_k \subset \{1,2,\dots, n\}\) denote the set of pages with a link to page \(k\), that is, \(L_k\) is the set of page \(k\)’s backlinks. For each \(k\) we require \[x_k = \sum_{j \in L_k} \frac{x_j}{n_j},\] where \(n_j\) is the number of outgoing links from page \(j\).

If we apply these scheme to the toy example in Figure 16, then for page 1 we have \(x_1 = x_3/1 + x_4/2\), since pages 3 and 4 are backlinks for page 1 and page 3 contains only one link, while page 4 contains two links (splitting its vote in half). Similarly, \(x_2 = x_1/3, x_3 = x_1/3 + x_2/2 + x_4/2\), and \(x_4 = x_1/3 + x_2/2\). These conditions can be expressed as linear system of equations \(Ax =x\), where \(x=[x_1,x_2,x_3,x_4]^T\) and \[A = \begin{bmatrix} 0 & 0 & 1 & \frac{1}{2} \\ \frac{1}{3} & 0 & 0 & 0 \\ \frac{1}{3} & \frac{1}{2} & 0 & \frac{1}{2} \\ \frac{1}{3} & \frac{1}{2} & 0 & 0 \end{bmatrix}\] Thus, we end up with an eigenvalue/eigenvector problem: Find the eigenvector \(x\) of the matrix \(A\), associated with the eigenvalue 1. We note that \(A\) is a column-stochastic matrix, since it is a square matrix for which all of its entries are nonnegative and the entries in each column sum to 1. If we build a random walk on the internet where each link is clicked with equal probability then \(A_{ij}\) is the probability that a random walked in page \(j\) goes to page \(i\). Later, in Chapter 7, we will use random walks on (undirected) graphs in order to embed their nodes in euclidean space, we note that, in that context, we will mostly work with \(M=A^T\) the matrix for which \(M_{ij}\) corresponds to the probability of taking step from \(i\) to \(j\). Stochastic matrices arise in the study of Markov chains and in a variety of modeling problems in economics and operations research. See e.g. [@HJ90] for more details on stochastic matrices. The fact that 1 is an eigenvalue of \(A\) is not just coincidence in this example, but holds true in general for stochastic matrices.

In our example, we obtain as eigenvector \(x\) of \(A\) associated with eigenvalue 1 the vector \(x = [x_1,x_2,x_3,x_4]^T\) with entries \(x_1 = \frac{12}{31}, x_2 = \frac{4}{31}, x_3 = \frac{9}{31}\), and \(x_4 = \frac{6}{31}\). Hence, perhaps somewhat surprisingly, page 3 is no longer the most important one, but page 1. This can be explained by the fact, that the in principle quite important page 3 (which has three webpages linking to it) has only one outgoing link, which gets all its “voting power”, and that link points to page 1.

In reality, \(A\) can easily be of size a billion times a billion. Fortunately, we do not need compute all eigenvectors of \(A\), only the eigenvector associated with the eigenvalue 1, which, as we know, is also the largest eigenvalue of \(A\). This in turn means we can resort to standard power iteration to compute \(x\) fairly efficiently (and we can also make use of the fact that \(A\) will be a sparse matrix, i.e., many of its entries will be zero). The actual PageRank algorithms adds some minor modifications,²³ but the essential idea is as described above.

The idea to use eigenvectors for ranking dates back to the late 1800s to the work of Edmund Landau in the context of ranking in chess tournaments, it was actually Landau’s first mathematical paper when he was 18 years old [@Landau1895Turnierresultate] (with a follow-up later [@Landau1914Preisverteilung]). You can read about this nice story in [@sinn2022landau].

6.2 Graphs and clustering↩︎

We now introduce the formalism for undirected²⁴ graphs, one of the main objects of study in what follows. A graph \(G = (V,E)\) contains a set of nodes \(V=\left\{v_1,\dots,v_n\right\}\) and edges \(E\subseteq {V \choose 2 }\). An edge \((i,j)\in E\) if \(v_i\) and \(v_j\) are connected. Figure 17 depicts one of the graph theorists’ favorite examples, the Petersen graph²⁵.

Let us recall some concepts about graphs that we will need.

• A graph is connected if, for all pairs of vertices, there is a path between these vertices in the graph. The number of connected components is simply the size of the smallest partition of the nodes into connected subgraphs. The Petersen graph is connected (and thus it has only \(1\) connected component).

• A clique of a graph \(G\) is a subset \(S\) of its nodes such that the subgraph corresponding to it is complete. In other words \(S\) is a clique if all pairs of vertices in \(S\) share an edge. The clique number \(c(G)\) of \(G\) is the size of the largest clique of \(G\). The Petersen graph has a clique number of \(2\).

• An independent set of a graph \(G\) is a subset \(S\) of its nodes such that no two nodes in \(S\) share an edge. Equivalently it is a clique of the complement graph \(G^c \mathrel{\mathop{:}}=\left(V,E^c\right)\). The independence number of \(G\) is simply the clique number of \(S^c\). The Petersen graph has an independence number of \(4\).

A particularly useful way to represent a graph is through its adjacency matrix. Given a graph \(G=(V,E)\) on \(n\) nodes (\(|V|=n\)), we define its adjacency matrix \(A \in \mathbb{R}^{n\times n}\) as the symmetric matrix with entries \[\label{adjacencymatrix} A_{ij} = \left\{ \begin{array}{ll} 1 & \text{ if } (i,j)\in E, \\ 0 & \text{otherwise}. \end{array} \right.\tag{88}\]

Sometimes, we will consider weighted graphs \(G=(V,E,W)\), where edges may have weights \(w_{ij}\) that are non-negative \(w_{ij}\geq 0\) and symmetric \(w_{ij} = w_{ji}\).

Much of the sequel will deal with graphs. Chapter 7 will treat (network) data visualization, dimension reduction, and embeddings of graphs in Euclidean space. Chapter 13 will introduce and study important random graph models. The rest of this Chapter will be devoted to spectral graph theory, clustering, and graph metrics.

Clustering is one of the central tasks in machine learning. Given a set of data points, or nodes of a graph, the purpose of clustering is to partition the data into a set of clusters where data points assigned to the same cluster correspond to similar data points (depending on the context, it could be for example having small distance to each other if the points are in Euclidean space, or having high connectivity if on a graph). We will start with an example of clustering points in Euclidean space, and later move back to graphs.

6.3 \(k\)-means clustering↩︎

One the most popular methods used for clustering is \(k\)-means clustering. Given \(x_1,\dots,x_n\in\mathbb{R}^p\) the \(k\)-means clustering partitions the data points into clusters \(S_1\cup \dots \cup S_k\) with centers \(\mu_1,\dots,\mu_k\in\mathbb{R}^p\) as the solution to: \[\label{eq:3:kmeans95obj} \min_{\text{partition}\;\substack{S_1,\dots,S_k \\ \mu_1,\dots,\mu_k}}\sum_{l=1}^k \sum_{i\in S_l} \left\| x_i - \mu_l\right\|^2.\tag{89}\] Note that, given the partition, the optimal centers are given by \[\mu_l = \frac{1}{\left| S_l \right|} \sum_{i\in S_l}x_i.\]

Lloyd’s algorithm [@Lloyd:kmeans] (also sometimes known as the \(k\)-means algorithm), is an iterative algorithm that alternates between

• Given centers \(\mu_1,\dots,\mu_k\), assign each point \(x_i\) to the cluster \(S_l\) with nearest center \[l=\operatorname{argmin}_{l=1,\dots, k}\left\| x_i - \mu_l \right\|.\]

• Update the centers \(\mu_l = \frac{1}{\left| S_l \right|} \sum_{i\in S_l}x_i\).

Unfortunately, Lloyd’s algorithm is not guaranteed to converge to the solution of 89 . Indeed, Lloyd’s algorithm oftentimes gets stuck in local optima of 89 . In the following we will discuss convex relaxations for clustering, which can be used as an alternative algorithmic approach to Lloyd’s algorithm, but since optimizing 89 is \(NP\)-hard there is no polynomial time algorithm that works in the worst-case (assuming the widely believed conjecture \(P\neq NP\), see also Chapter [maxcutapprox])

While popular, \(k\)-means clustering has some potential issues:

• One needs to set the number of clusters a priori (a typical way to overcome this issue is by trying the algorithm for different number of clusters).

• The way 89 is defined it needs the points to be defined in an Euclidean space, oftentimes we are interested in clustering data for which we only have some measure of affinity between different data points, but not necessarily an embedding in \(\mathbb{R}^p\) (this issue can be overcome by reformulating 89 in terms of distances only).

• The formulation is computationally hard, so algorithms may produce suboptimal instances.

• The solutions of \(k\)-means are always convex clusters. This means that \(k\)-means may have difficulty in finding cluster such as in Figure 19.

6.4 Spectral clustering↩︎

A natural way to try to overcome the issues of \(k\)-means depicted in Figure 19 is by transforming the data into a graph and cluster the graph: Given the data points we can construct a weighted graph \(G = (V,E,W)\) using a similarity kernel \(K_{\varepsilon}\), such as \(K_{\varepsilon}(u) = \exp\left(-\frac{1}{2\varepsilon}u^2\right)\), by associating each point to a vertex and, for which pair of nodes, set the edge weight as \[w_{ij} = K_{\varepsilon}\left( \|x_i-x_j\| \right).\] Another popular procedure to transform data into a graph is by constructing a graph where data points are connected if they correspond to the nearest neighbors. We note that this procedure only needs a measure of distance, or similarity, of data points and not necessarily that they lie in an Euclidean space. Given this motivation, and the prevalence of network data, we will now address the problem of clustering the nodes of a graph.

Normalized Cut

Given a graph \(G=(V,E,W)\), the goal is to partition the graph into clusters in a way that keeps as many of the edges, or connections, within the clusters and has as few edges as possible across clusters. We will focus on the case of two clusters, and briefly address extensions at the end of this chapter. A natural way to measure a vertex partition \(\left(S,S^c\right)\) is \[\operatorname{cut}(S) = \sum_{i\in S}\sum_{j\in S^c} w_{ij}.\]

If we represent the partition by a vector \(y\in\{\pm 1\}^n\) where \(y_i=1\) if \(i\in S\), and \(y_i=-1\) otherwise, then the cut is a quadratic form on the Graph Laplacian.

Note that the entries of \(L_G\) are given by \[\left(L_G\right)_{ij} = \left\{ \begin{array}{ll} -w_{ij} & \text{ if } i\neq j \\ \deg(i) & \text{ if } i=j. \end{array} \right.\]

If \(S\subset V\) and \(y\in\{\pm 1\}^n\) such that \(y_i=1\) if \(i\in S\), and \(y_i=-1\) otherwise, then it is easy to see that \[\operatorname{cut}(S) = \frac{1}{4} \sum_{i<j}w_{ij}(y_i-y_j)^2.\] The following proposition establishes \[\label{eq:cut95Laplacian} \operatorname{cut}(S) = \frac{1}{4} y^TL_Gy,\tag{92}\] for \(y\in\{\pm 1\}^n\) such that \(y_i=1\) if and only if \(i\in S\).

Note that Proposition 6.2 implies that the graph Laplacian \(L_G\) is a positive semidefinite (PSD) matrix, since \(x^T L x \geq 0\) for all \(x\) (assuming non-negative weights \(w_{ij} \geq 0\)). This property can also be deduced directly from Definition 6.1 of \(L_G\) as the (non-negative) weighted sum of rank-1 PSD matrices.

While \(\operatorname{cut}(S)\) is a good way of measuring the fit of a partition, it has a major drawback: the minimum cut is achieved for \(S=\emptyset\) (since \(\operatorname{cut}(\emptyset) = 0\)) which is a rather meaningless choice of partition. Constraining the partition to be non-trivial would not overcome this drawback, because very unbalanced partitions would still be favored (e.g., partitions with \(|S|=1\) containing just a single node). Below we discuss how to promote more balanced partitions.

Asking for the partition to be exactly balanced is too restrictive in many cases. There are several ways to evaluate a partition that are variations of \(\operatorname{cut}(S)\) that take into account the intuition that one wants both \(S\) and \(S^c\) not to be too small (although not necessarily equal to \(|V|/2\)). A prime example is Cheeger’s cut.

A similar object is the Normalized Cut, \(\operatorname{Ncut}\), which is given by \[\operatorname{Ncut}(S) = \frac{\operatorname{cut}(S)}{\operatorname{vol}(S)} + \frac{\operatorname{cut}(S^c)}{\operatorname{vol}(S^c)}.\]

Note that \(\operatorname{Ncut}(S)\) and \(h(S)\) are tightly related, in fact it is easy to see that: \[h(S) \leq \operatorname{Ncut}(S) \leq 2h(S).\]

6.4.0.1 Normalized Cut as a spectral relaxation

Below we will show that \(\operatorname{Ncut}\) can be written in terms of a minimization of a quadratic form involving the graph Laplacian \(L_G\), analogously to the balanced partition as described in Remark 6.1.

Recall that the cut value of a balanced partition can be written as \[\frac{1}{4} \min_{\substack{y\in\{-1,1\}^n \\ \mathbf{1}^Ty = 0}} y^TL_Gy.\] An intuitive way to relax the balanced condition is to allow the labels \(y\) to take values in two different real values \(a\) and \(b\) (e.g. \(y_i = a\) if \(i\in S\) and \(y_j = b\) if \(i\notin S\)) but not necessarily \(\pm1\). We can then use the notion of volume of a set to ensure a less restrictive notion of balanced by asking that \[\label{eq:1TDy0} a\operatorname{vol}\left(S\right) + b\operatorname{vol}\left(S^c\right) = 0,\tag{94}\] where \[\label{eq:vol} \operatorname{vol}(S) = \sum_{i\in S} \deg(i).\tag{95}\] Thus 94 corresponds to \(\mathbf{1}^TDy = 0\).

We also need to fix a scale for \(a\) and \(b\): \[a^2\operatorname{vol}\left(S\right) + b^2\operatorname{vol}\left(S^c\right) = 1,\] which corresponds to \(y^TDy = 1\).

This suggests considering \[\min_{\substack{y\in\{a,b\}^n \\ \mathbf{1}^TDy = 0,\, y^TDy = 1}} y^TL_Gy.\]

As we will see below, this corresponds precisely to \(\operatorname{Ncut}\).

This means that finding the minimum \(\operatorname{Ncut}\) corresponds to solving

\[\label{eq:3:Ncut:norelax} \begin{array}{rl} \min & y^TL_Gy \\ \text{s. t.} & y\in \{a,b\}^n \text{ for some } a \text{ and } b \\ & y^TDy = 1 \\ & y^TD\mathbf{1} = 0. \end{array}\tag{96}\]

Since solving 96 is, in general, NP-hard, we consider a similar problem where the constraint that \(y\) can only take two values is removed:

\[\label{eq:3:Ncut:relaxed} \begin{array}{rl} \min & y^TL_Gy \\ \text{s. t.} & y\in \mathbb{R}^n \\ & y^TDy = 1 \\ & y^TD\mathbf{1} = 0. \end{array}\tag{97}\]

Given a solution of 97 we can round it to a partition by setting a threshold \(\tau\) and taking \(S = \left\{ i\in V:\, y_i \leq \tau \right\}\). We will see below that \(\eqref{eq:3:Ncut:relaxed}\) is an eigenvector problem (for this reason we call 97 a spectral relaxation).

In order to better see that 97 is an eigenvector problem (and thus computationally tractable), set \(z = D^{\frac{1}{2}}y\) and \[\label{eq:LLL95G} \mathcal{L}_G = D^{-\frac{1}{2}} L_G D^{-\frac{1}{2}},\tag{98}\] then 97 is equivalent

\[\label{eq:3:Ncut:relaxed95LLL} \begin{array}{rl} \min & z^T\mathcal{L}_Gz \\ \text{s. t.} & z\in \mathbb{R}^n \\ & \|z\|^2 = 1 \\ & \left(D^{\frac{1}{2}}\mathbf{1}\right)^Tz = 0. \end{array}\tag{99}\]

Note that \(\mathcal{L}_G = I - D^{-\frac{1}{2}}WD^{-\frac{1}{2}}\). We order its eigenvalues in increasing order \(0=\lambda_1\left( \mathcal{L}_G\right) \leq \lambda_2\left( \mathcal{L}_G\right) \leq \cdots \leq \lambda_n\left( \mathcal{L}_G\right)\). The eigenvector associated with the smallest eigenvalue is given by \(D^{\frac{1}{2}}\mathbf{1}\). By the variational interpretation of the eigenvalues, it follows that the solution of 99 is \(\lambda_2\left( \mathcal{L}_G\right)\) and the minimizer is given by the second smallest eigenvector of \(\mathcal{L}_G = I - D^{-\frac{1}{2}}WD^{-\frac{1}{2}}\), which we call \(v_2\). Note that this corresponds also to the second largest eigenvector of \(D^{-\frac{1}{2}}WD^{-\frac{1}{2}}\). This means that the optimal \(y\) in 97 is given by \(\varphi_2 = D^{-\frac{1}{2}}v_2\). This motivates Algorithm 20, which is often referred to as Spectral Clustering:

6.5 Cheeger’s Inequality↩︎

Because the relaxation 97 is obtained from 96 by removing a constraint we immediately have that \[\lambda_2\left(\mathcal{L}_G\right) \leq \min_{S\subset V} \operatorname{Ncut}(S).\] This means that \[\frac{1}{2}\lambda_2\left(\mathcal{L}_G\right) \leq h_G.\]

In what follows we will show a performance guarantee for Algorithm 20.

This implies in particular that \[h(S) \leq \sqrt{4h_G},\] meaning that Algorithm 20 is suboptimal at most by a square-root factor. Note that this also directly implies the famous Cheeger’s Inequality, stated next.

Cheeger’s inequality was first established for manifolds by Jeff Cheeger in 1970 [@JCheeger_1970], and the graph version is due to Noga Alon and Vitaly Milman [@NAlon_1986; @NAlon_VMilman_1986] in the mid 80s. Cheeger’s inequality for a closed manifold \(M\) gives a bound on the area of a hypersurface \(E\), that partitions \(M\) into two disjoint parts.

6.5.0.1 Bounding the diameter of a graph using \(\lambda_2(\mathcal{L}_G)\)

Before proving Cheeger’s inequality for graphs, we first discuss another eigenvalue inequality, that gives a flavor of the connection between the second eigenvalue of the normalized graph Laplacian and the geometry of the graph (see also [@FanChung_SpectralGraphTheory]).

Associate a cost \(c(u,v)\) for each pair of nodes inversely proportional to the weight between them: \[c(u,v) = \left\{\begin{array}{cc} \frac{1}{w_{uv}} & \text{for } (u,v)\in E, \\ & \\ +\infty & \text{for } (u,v)\notin E. \end{array} \right.\] A path of length \(k\) connecting \(u\) and \(v\) is of the form \[P = (u=v_1,v_2,\ldots,v_k=v).\] The cost associated with this path is \[c(P) = \sum_{i=1}^{k-1} c(v_i,v_{i+1}).\] The geodesic distance or shortest path between \(u\) and \(v\) is the minimum of the cost over all possible paths connecting them: \[d_g(u,v) = \min \{c(P) \,| \, P(1)=u,\, P(k)=v,\; |P|=k\}\] The diameter of the graph \(G\) is the largest geodesic distance \[\text{diam}(G) = \max_{u,v} d_g(u,v).\]

We demonstrate that for \(\lambda_2(\mathcal{L}_G)\), the second smallest eigenvalue of the normalized graph Laplacian \(\mathcal{L}_G = D^{-1/2}L_G D^{-1/2}\), \[\label{eq:diam-lambda} \text{diam}(G) \geq \frac{ 1}{\text{vol}(G)\lambda_2(\mathcal{L}_G)}.\tag{100}\] That is, a small second eigenvalue implies that the diameter of the graph must be large. We already know that \(\lambda_2(\mathcal{L}_G) = 0\) implies that the graph is disconnected and the diameter is infinite. The bound (100 ) is more quantitative from that standpoint.

6.5.0.2 The proof of Cheeger’s Inequality

The upper bound in Cheeger’s inequality (corresponding to Lemma 6.1) is more interesting but more difficult to prove, and it is often referred to as the “the difficult part” of Cheeger’s inequality. We will prove this Lemma in what follows. There are several proofs of this inequality (see [@Chung_CheegersIneq] for four different proofs!). The proof that follows is an adaptation of the proof in this blog post [@trevisan:blog:cheegerproof] for the case of weighted graphs.

Multiple Clusters

Much of the above can be easily adapted to multiple clusters. Algorithm 21 is a natural extension of spectral clustering to multiple clusters.²⁷

There is also an analogue of Cheeger’s inequality. A natural way of evaluating \(k\)-way clustering is via the \(k\)-way expansion constant (see [@JRLee_SOGharan_LTrevisan_2011]): \[\rho_G(k) = \min_{S_1,\dots,S_k} \max_{l=1,\dots,k} \left\{ \frac{\operatorname{cut}(S_l)}{\operatorname{vol}(S_l)} \right\},\] where the maximum is over all choice of \(k\) disjoint subsets of \(V\) (but not necessarily forming a partition).

Another natural definition is \[\varphi_G(k) = \min_{S: \operatorname{vol}{S}\leq \frac{1}{k} \operatorname{vol}(G)} \frac{\operatorname{cut}(S)}{\operatorname{vol}(S)}.\] It is easy to see that \[\varphi_G(k) \leq \rho_G(k).\] The following are analogues of Cheeger’s inequality for multiple clusters.

Exercises↩︎

- the smallest eigenvalue $\lambda_1(L)$ is always zero;

- the second smallest eigenvalue $\lambda_2(L)$ is zero if and only
  if $G$ is disconnected;

- the nullity of $L$ (i.e. dimension of the kernel $\ker(L)$) equals
  the number of components of $G$.

Exercise 6.11 (Tightness of the lower bound in Cheeger’s inequality). Let \(d \geq 2\) be an integer, \(G = (V,E)\) be the \(d\)-dimensional hypercube, and \(\mathcal{L}_G\) its normalized Laplacian. We index the \(n=2^d\) vertices of the hypercube by \(d\)-dimensional \(\left\{ 0, 1 \right\}\)-vectors, i.e. \(V = \left\{ 0,1 \right\}^d\), and given any \(x, y \in V\), we have \(\left\{ x,y \right\} \in E\) if and only if \(x\) and \(y\) differ in exactly one coordinate. The example for \(d = 3\) is given in Figure 23.

Given a (possibly empty) subset \(T \subseteq [d]\) let \(v_T \in \mathbb{R}^n\) be a vector, whose coordinates are indexed by \(n=2^d\) vertices of the hypercube and defined by \[v_T(x) = (-1)^{\sum_{i \in T} x_i},\] where \(x_i\) is the \(i\)-th coordinate of the vertex \(x\in\{0,1\}^d\). Also, let \(S_T \subseteq V\) be the subset of vertices given by \[S_T = \left\{ x \in V \colon v_T(x) = 1 \right\}.\] [When \(T = \emptyset\), we interpret the empty sum as zero, i.e. \(\sum_{i \in \emptyset}x_i = 0\).]

(a) Compute \(h(S_{\left\{ 1 \right\}})\), the Cheeger’s cut of the subset \(S_{\left\{ 1 \right\}}\).

(b) Show that every vertex in \(S_T\) has exactly \(d-|T|\) neighbours in \(S_T\).

(c) Show that for any \(T \subseteq [d]\), \(v_T\) is an eigenvector of \(\mathcal{L}_G\) with eigenvalue \(\frac{2|T|}{d}\).

(d) Show that if \(T' \subseteq [d]\) is distinct from \(T\), then \(v_T\) and \(v_{T'}\) are orthogonal.

(e) Conclude that for any \(0 \leq k \leq d\), the eigenspace corresponding to the eigenvalue \(\frac{2k}{d}\) has dimension \({d \choose k}\).

(f) Compute \(h_G\), the Cheeger’s constant of \(G\).

7 Nonlinear Dimension Reduction and Diffusion Maps↩︎

In Chapter 4 we discussed dimension reduction via Principal Component Analysis. Many datasets however have low dimensional structure that is not linear. In this chapter we will discuss nonlinear dimension reduction techniques. Just as with spectral clustering in Chapter 6 we will focus on graph data while noting that most types of data can be transformed into a weighted graph by means of a similarity kernel (Section 7.1.1). The goal of this chapter is to embed the nodes of a graph in Euclidean space in a way that best preserves the intrinsic geometry of the graph (or the data that gave rise to the graph).

7.1 Diffusion maps↩︎

Diffusion maps [@RRCoifman_SLafon_2006] will allow us to represent (weighted) graphs \(G=(V,E,W)\) in \(\mathbb{R}^d\), i.e. associating, to each node, a point in \(\mathbb{R}^d\). Before presenting diffusion maps, we will introduce a few important notions. The reader may notice the similarities with the objects described in the context of PageRank in Chapter 6, the main difference is that here the connections between graphs have no direction, meaning that the weight matrix \(W\) is symmetric²⁸; this will be crucial in the derivations below (as we will see in Section 7.7 this symmetry can be viewed as corresponding to a certain Markov Chain being reversible).

Given \(G=(V,E,W)\) we consider a random walk (with independent steps) on the vertices of \(V\) with transition probabilities: \[{\mathbb{P}}\left\{ X(t+1) = j | X(t) = i \right\} = \frac{w_{ij}}{\deg(i)},\] where \(\deg(i) = \sum_j w_{ij}\). Let \(M\) be the matrix of these probabilities, \[\label{eq:def:M} M_{ij} = \frac{w_{ij}}{\deg(i)}.\tag{105}\] It is easy to see that \(M_{ij}\geq 0\) and \(M\mathbf{1} = \mathbf{1}\) (indeed, \(M\) is a transition probability matrix). Recalling that \(D\) is the diagonal matrix with diagonal entries \(D_{ii} = \deg(i)\) we have \[M = D^{-1}W.\]

If we start a random walker at node \(i\) (\(X(0)=1\)) then the probability that, at step \(t\), is at node \(j\) is given by \[{\mathbb{P}}\left\{ X(t) = j | X(0) = i \right\} = \left(M^t\right)_{ij}.\] In other words, the probability cloud of the random walker at point \(t\), given that it started at node \(i\) is given by the row vector \[{\mathbb{P}}\left\{ X(t) | X(0) = i \right\} = e_i^TM^t = M^t[i,:].\]

\(M\) is not symmetric, but a matrix similar to M, \(S = D^{\frac{1}{2}} M D^{-\frac{1}{2}}\) is, indeed \(S = D^{-\frac{1}{2}}WD^{-\frac{1}{2}}\). We consider the spectral decomposition of \(S\) \[\label{eq:SisVlambdaVT} S = V\Lambda V^T,\tag{106}\] where \(V= \left[ v_1,\dots,v_n \right]\) satisfies \(V^TV = I_{n\times n}\) and \(\Lambda\) is diagonal with diagonal elements \(\Lambda_{kk} = \lambda_k\) (and we organize them as \(\lambda_1\geq \lambda_2 \geq \cdots \geq \lambda_n\)). Note that \(Sv_k = \lambda_k v_k\). Also, \[M = D^{-\frac{1}{2}} S D^{\frac{1}{2}} = D^{-\frac{1}{2}} V\Lambda V^T D^{\frac{1}{2}} = \left( D^{-\frac{1}{2}} V\right) \Lambda \left( D^{\frac{1}{2}} V \right)^T.\] We define \(\Phi = D^{-\frac{1}{2}} V\) with columns \(\Phi = \left[ \varphi_1,\dots,\varphi_n\right]\) and \(\Psi = D^{\frac{1}{2}} V\) with columns \(\Psi = \left[ \psi_1,\dots,\psi_n\right]\). Then \[\label{eq:MisPhiLambdaPsiT} M = \Phi \Lambda \Psi^T,\tag{107}\] and \(\Phi,\, \Psi\) form a biorthogonal system in the sense that \(\Phi^T\Psi = I_{n\times n}\) or, equivalently, \(\varphi_j^T\psi_k = \delta_{jk}\). Note that \(\varphi_k\) and \(\psi_k\) are, respectively right and left eigenvectors of \(M\), indeed, for all \(1\leq k\leq n\): \[M\varphi_k = \lambda_k \varphi_k \quad \text{ and } \psi_k^T M = \lambda_k \psi_k^T.\] Also, we can rewrite this decomposition as \[M = \sum_{k=1}^n \lambda_k \varphi_k \psi_k^T.\] and it is easy to see that \[\label{eq:Masrank1withphiandpsi} M^t = \sum_{k=1}^n \lambda_k^t \varphi_k \psi_k^T.\tag{108}\]

Let us revisit the embedding suggested on Remark 7.1. It would correspond to \[i \to M^t[i,:] = \sum_{k=1}^n \lambda_k^t \varphi_k(i) \psi_k^T,\] it is written in terms of the basis \(\psi_k\). The diffusion map will essentially consist of the representing a node \(i\) by the coefficients of the above map \[\label{eq:2:DM:withcoordinate1} i \to \left[ \begin{array}{c} \lambda_1^t \varphi_1(i) \\ \lambda_2^t \varphi_2(i) \\ \vdots \\ \lambda_n^t \varphi_n(i) \end{array} \right],\tag{109}\] Note that \(M\mathbf{1} = \mathbf{1}\), meaning that one of the right eigenvectors \(\varphi_k\) is simply a multiple of \(\mathbf{1}\) and so it does not distinguish the different nodes of the graph. We will show that this indeed corresponds to the the first eigenvalue.

By the proposition above we can take \(\varphi_1 = \mathbf{1}\), meaning that the first coordinate of 109 does not help differentiate points on the graph. This suggests removing that coordinate:

This map is still a map to \(n-1\) dimensions. But note now that each coordinate has a factor of \(\lambda_k^t\) which, if \(\lambda_k\) is small will be rather small for moderate values of \(t\). This motivates truncating the diffusion map by taking only the first \(d\) coefficients.

We note that the truncated diffusion map without scaling the coordinates by \(\lambda^t\) (or equivalently, taking \(t=0\)) is known as Laplacian Eigenmap, introduced by Belkin and Niyogi [@belkin2001laplacian; @belkin2003laplacian].

In the following theorem we show that the Euclidean distance in the diffusion map coordinates (called diffusion distance) meaningfully measures distance between the probability clouds after \(t\) iterations.

As simple example, let us consider the cycle graph (or ring graph) \(C_n\), which is a graph on \(n\) nodes \(\{1,\dots,n\}\) such that node \(k\) is connected to \(k-1\) and \(k+1\) and \(1\) is connected to \(n\). The diffusion map of \(C_{10}\) truncated to two dimensions, shown in Figure 24, gives a very natural way of displaying the graph (indeed, if one is asked to draw the ring graph, this is probably the drawing that most people would do). It is not difficult to analytically compute the diffusion map of \(C_n\) and confirm that it displays the points in a circle.

7.1.1 Diffusion maps of point clouds↩︎

Very often we are interested in embedding in \(\mathbb{R}^d\) a point cloud of points \(x_1,\dots,x_n \in \mathbb{R}^p\) and not necessarily a graph. One option is to use Principal Component Analysis (PCA), but PCA is only designed to find linear structure of the data and the low dimensionality of the dataset may be non-linear. For example, let us say that our dataset is images of the face of a person taken from different angles and lighting conditions, for example, the dimensionality of this dataset is limited by the amount of muscles in the head and neck and by the degrees of freedom of the lighting conditions but it is not clear that this low dimensional structure is linearly apparent on the pixel values of the images.

Let us consider a point cloud that is sampled from a two dimensional swiss roll embedded in three dimension (see Figure 25). In order to learn the two dimensional structure of this object we need to differentiate points that are near each other because they are close by in the manifold and not simply because the manifold is curved and the points appear nearby even when they really are distant in the manifold (see Figure 25 for an example). We will achieve this by creating a graph from the data points.

Our goal is for the graph to capture the structure of the manifold. To each data point we will associate a node. For this we should only connect points that are close in the manifold and not points that maybe appear close in Euclidean space simply because of the curvature of the manifold. This is achieved by picking a small scale and linking nodes if they correspond to points whose distance is smaller than that scale. This is usually done smoothly via a kernel²⁹ \(K_\varepsilon\), and to each edge \((i,j)\) associating a weight \[w_{ij} = K_\varepsilon\left( \|x_i - x_j\|_2 \right),\] a common example of a kernel is \(K_\varepsilon(u) = \exp\left( -\frac{1}{2\varepsilon} u^2 \right)\), that gives essentially zero weight to edges corresponding to pairs of nodes for which \(\|x_i-x_j\|_2 \gg \sqrt{\varepsilon}\). We can then take the the diffusion maps of the resulting graph.

7.1.2 An illustrative simple example and noisy data↩︎

A simple and illustrative example is to take images of a blob on a background in different positions (imagine a white square on a black background and each data point corresponds to the same white square in different positions). This dataset is clearly intrinsically two dimensional, as each image can be described by the (two-dimensional) position of the square. However, we do not expect this two-dimensional structure to be directly apparent from the vectors of pixel values of each image; in particular we do not expect vectors of images with this type of structure to lie in a two dimensional affine subspace!³⁰

Let us start by experimenting with the above example for one dimension. In that case the blob is a vertical stripe and simply moves left and right. We think of our space as the one in many arcade games, if the square or stripe moves to the right all the way to the end of the screen, it shows up on the left side (and same for up-down in the two-dimensional case). Not only should this point cloud have a one dimensional structure but it should also exhibit a circular structure. Remarkably, this structure is completely apparent when taking the two-dimensional diffusion map of this dataset, see Figure 26.

For the two dimensional example, we expect the structure of the underlying manifold to be a two-dimensional torus. Indeed, Figure 27 shows that the three-dimensional diffusion map captures the toroidal structure of the data.

In many applications, the data are noisy and the Euclidean distances \(\|x_i - x_j\|\) would be affected by the measurement noise. It is therefore common to apply a denoising algorithm to the data prior to the application of diffusion maps. A popular choice is to first denoise and reduce dimensionality using PCA, followed by diffusion maps in the reduced space of principal components. The dimension of PCA is chosen large enough to retain the non-linear structure of the data and the signal variability, but small enough to effectively diminish the noise. See [@singer2013two; @gilles2025cryo] for the application of PCA followed by diffusion maps in biomedical imaging. Fitting a manifold in the presence of noise is a topic of ongoing research [@aizenbud2025estimation; @fefferman2025fitting].

7.1.3 Diffusion maps, images, and data manifolds↩︎

To illustrate the ability of the diffusion map to recover physically meaningful latent structure from image data, we experimentally study its performance on two standard image datasets, MNIST and the Yale Face Database.

In both experiments, the affinity between two images \(x_i, x_j\) is defined by the Gaussian kernel \[W_{ij} = \exp\!\left(-\frac{\|x_i - x_j\|^2}{2\varepsilon}\right),\] where \(\varepsilon > 0\) is the bandwidth parameter. Here we choose \(\varepsilon\) proportional to the median of the pairwise squared distances: \(\varepsilon = \frac{1}{2}\mathrm{median}\bigl\{\|x_i - x_j\|^2 : i \neq j\bigr\}\). The median heuristic (as well as the mean heuristic) is a standard data-adaptive choice that scales \(\varepsilon\) to the typical squared distance between points in the dataset, avoiding the need to tune the bandwidth by hand.

7.1.3.1 Example: Diffusion maps on MNIST handwritten digits

To exhibit the geometry recovered by the diffusion map, we apply it to a subset of the MNIST database. Each image is represented as a vector in \(\mathbb{R}^{784}\) by flattening its pixel intensities, so the data matrix \(X \in \mathbb{R}^{784 \times 600}\) has images as columns and pixel intensities as rows.

The resulting embeddings are shown in Figure 28 and respectively. Figure 28 depicts the diffusion map embedding obtained from the first two diffusion coordinates \(\varphi_1, \varphi_2\) for all 600 images, with a subset of the images displayed as a small thumbnail at its embedding location and colored borders indicating digit class.

The four digit classes separate quite well in the \((\varphi_1, \varphi_2)\) plane without any use of class labels during training: the diffusion map discovers the class structure from geometry alone. This gives support to the assumption that images of the same digit class lie on a lower-dimensional manifold within \(\mathbb{R}^{784}\), corresponding to the natural variability in handwriting style, stroke width, and slant within that class. Within each digit class the images are organized continuously: adjacent points in the embedding correspond to images that look similar, and one can trace smooth trajectories through each cluster along which the appearance changes gradually.

Figure 29 shows the the diffusion map embedding corresponding to the first three diffusion coordinates \(\varphi_1, \varphi_2, \varphi_3\), which reveals a good (but certainly not perfect) unsupervised clustering of the data points according to their digit labels (here, 2 = yellow, 3 = purple, 4 = blue, 6 = red). It is important to emphasize that two and three diffusion coordinates were chosen for visualization purposes, but depending on the task at hand more diffusion dimensions may be needed to accomplish the desired goal. See also the discussion in Chapter 7.4.1 about using diffusion maps as a preprocessing step.

7.1.3.2 Example: Diffusion maps on the Extended Yale B Face Database

In this experiment we apply the diffusion map to the Extended Yale B database [@georghiades2001few], a controlled face image dataset consisting of \(38\) subjects photographed under \(64\) distinct illumination conditions. Each image has been cropped and normalized to \(192 \times 168\) pixels. For this experiment we select a single subject and work with all \(64\) of their frontal-pose images, representing each as a vector in \(\mathbb{R}^{896}\) obtained by resizing to \(32 \times 28\) pixels and flattening.

The illumination conditions are parameterized by two angles: the azimuth \(\phi \in [-130^\circ, +130^\circ]\) and the elevation \(\theta \in [-40^\circ, +90^\circ]\) of the light source relative to the camera axis. As the light source moves around the face, the image intensity pattern changes dramatically: under extreme lateral illumination large portions of the face fall into shadow, while under frontal illumination the face appears nearly uniform. Despite this large photometric variation, all \(64\) images depict the same individual with the same geometry, so they lie on a smooth two-dimensional manifold in \(\mathbb{R}^{896}\) parameterized by \((\phi, \theta)\).

The embedding provided by the diffusion map reveals the intrinsic two-dimensional structure of the illumination manifold, see Figure 30. The first diffusion coordinate \(\varphi_1\) correlates strongly with the azimuth of the light source: moving along \(\varphi_1\) corresponds to sweeping the light from left to right across the face, with the shadow pattern shifting accordingly. The second diffusion coordinate \(\varphi_2\) captures the elevation of the light source: moving along \(\varphi_2\) corresponds to raising or lowering the light from below to above the face.

7.2 Connections between diffusion maps and spectral clustering↩︎

Diffusion maps are tightly connected to spectral clustering (described in Chapter 6). In fact, spectral clustering can be understood as simply performing \(k\)-means on the embedding given by diffusion maps truncated to \(k-1\) dimensions.

A natural way to try to overcome the issues of \(k\)-means depicted in Figure 19 is by using diffusion maps: Given the data points we construct a weighted graph \(G = (V,E,W)\) using a kernel \(K_{\varepsilon}\), such as \(K_{\varepsilon}(u) = \exp\left(-\frac{1}{2\varepsilon}u^2\right)\), by associating each point to a vertex and, for which pair of nodes, set the edge weight as \[w_{ij} = K_{\varepsilon}\left( \|x_i-x_j\| \right).\]

Recall the construction of a matrix \(M = D^{-1}W\) as the transition matrix of a random walk \[{\mathbb{P}}\left\{ X(t+1) = j | X(t) = i \right\} = \frac{w_{ij}}{\deg(i)} = M_{ij},\] where \(D\) is the diagonal with \(D_{ii} = \deg(i)\). The \(d\)-dimensional diffusion map is given by \[\mathcal{D}_t^{(d)}(i) = \left[ \begin{array}{c} \lambda_2^t \varphi_2(i) \\ \vdots \\ \lambda_{d+1}^t \varphi_{d+1}(i) \end{array} \right],\] where \(M = \Phi \Lambda \Psi^T\) where \(\Lambda\) is the diagonal matrix with the eigenvalues of \(M\) and \(\Phi\) and \(\Psi\) are, respectively, the right and left eigenvectors of \(M\) (note that they form a bi-orthogonal system, \(\Phi^T\Psi = I\)).

If we want to cluster the vertices of the graph in \(k\) clusters, then it is natural to truncate the diffusion map to have \(k-1\) dimensions (since in \(k-1\) dimensions we can have \(k\) linearly separable sets). If indeed the clusters were linearly separable after embedding then one could attempt to use \(k\)-means on the embedding to find the clusters, this is precisely the motivation for spectral clustering.

In order to show that this indeed coincides with Algorithm 21, it is enough to show that \(\varphi_m = D^{-\frac{1}{2}}v_m\) where \(v_m\) is the eigenvector associated with the \(m\)-th smallest eigenvalue of \(\mathcal{L}_G\). This follows from the fact that \(S=D^{-\frac{1}{2}}WD^{-\frac{1}{2}}\) as defined in 106 is related to \(\mathcal{L}_G\) by \[\mathcal{L}_G = I - S,\] and \(\Phi = D^{-1/2}V\).

Proposition 7.3 below establishes a connection between \(\operatorname{Ncut}\) (as described in Chapter 6) and the random walks introduced above. Let \(M\) as defined in 105 denote the matrix of transition probabilities. Recall that \(M\mathbf{1} = \mathbf{1}\), corresponding to \(M\varphi_1 = \varphi_1\), which means that \(\psi_1^TM =\psi_1^T\), where \[\psi_1 = D^{\frac{1}{2}}v_1 = D\varphi_1 = \left[\deg(i) \right]_{1\leq i\leq n}.\] This means that \(\left[\frac{ \deg(i)}{ \operatorname{vol}(G)} \right]_{1\leq i\leq n}\) is the stationary distribution of this random walk. Indeed it is easy to check that, if \(X(t)\) has a certain distribution \(p_t\) then \(X(t+1)\) has a distribution \(p_{t+1}\) given by \(p_{t+1}^T = p_{t}^TM\)

7.3 Commute time distance↩︎

We saw that the diffusion distance is a way to equip the graph with a metric that has a nice probabilistic interpretation. A natural question is if there are other meaningful and useful metrics for graphs. Perhaps the first metric that comes to mind is the shortest path metric, or geodesic distance. A cost \(c_{ij}\) is associated with each edge, e.g., \(c_{ij} = w_{ij}^{-1}\). The cost of a path is the sum of costs along the path edges. The geodesic distance between \(i\) and \(j\) is the smallest possible cost over all possible paths connecting \(i\) and \(j\). It can be verified that the geodesic distance is indeed a metric (i.e., it satisfies all properties required from a metric), but it is not necessarily an \(\ell^2\) metric. Not every metric can be embedded in a Euclidean (or in general, Hilbert) space. For example, the \(\ell^1\) distance (Manhattan distance) cannot be embedded in a Euclidean space (it does not satisfy the parallelogram identity). Also, the geodesic distance has other shortcomings that will be discussed later in the course, such as its sensitivity to noise. We remark that the geodesic distance is the basis for a remarkable non-linear dimensionality reduction method called Isomap [@ISOMAP_paper_2000], which will be revisited later on.

Are there any other interesting/natural graph metrics? Given the random walk process over the graph, it is natural to consider the average time it takes to get from node \(i\) to node \(j\). It can be easily seen to satisfy the triangle inequality. However, it is not symmetric: the average time to get from \(i\) to \(j\) may not equal the time from \(j\) to \(i\). So, we can symmetrize it. Clearly, the time to get from a node to itself is 0, so overall we get a metric. It is called the commute time distance. Surprisingly, a kind of miracle happens: The commute time distance turns out to be a Euclidean distance corresponding to an embedding that is also given in terms of eigenvectors and eigenvalues of the graph Laplacian. Our goal is to prove this by starting from the probabilistic interpretation (see also [@doyle1984random; @lovasz1993random]).

Consider an irreducible random walk \(\{X(t)\}_{t=0}^\infty\) with jump probabilities \(M_{ij} = {\mathbb{P}}\left\{ X(t+1) = j | X(t) = i \right\} = w_{ij}/\deg(i)\) where \(w_{ij}\) are the weights of a undirected graph; irreducibility is then equivalent to the graph being connected. We define the first passage time and mean first passage time as follows:

First Passage Time

\[\tau_{ij}=\left\{\begin{array}{cc} \inf\left\{ t>0 \,|\, X(t)={x}_{j}, \; X(0)={x}_{i}\right\} & fori\neq j, \\ &\\ 0 & fori=j. \end{array} \right.\]

Mean First Passage Time

\[t_{ij}=\mathbb{E}[\tau_{ij}].\] The expected commute time between \({x}_{i}\) and \({x}_{j}\) is defined as \[t_{ij}+t_{ji}.\] Next, we show that the expected commute time is a metric. Let \[\tau'_{ij}=\inf\left\{ t>0 \,|\, X(t)={x}_{j},\,X(0)={x}_{i}\right\}.\] Notice that \(\tau_{ii} \neq \tau'_{ii}\) because \(\tau'_{ii} > 0.\) Let \[t'_{ij}=\mathbb{E}[\tau'_{ij}].\] The random walker can get from \(x_i\) to \(x_j\) either directly in one time step with probability \(M_{ij}=w_{ij}/\deg(i)\), or it first goes from \(x_i\) to \(x_k\) (\(k\neq j\)) with probability \(M_{ik}\) and the expected time to get to \(x_j\) is \(1+t'_{kj}\). Let \(T\) and \(T'\) denote the matrices with entries respectively \(t_{ij}\) and \(t'_{ij}\). We have,³¹ \[\begin{align} T'_{ij}&=&M_{ij}\cdot1+\sum_{k\neq j}M_{ik}(1+T'_{kj}) \nonumber \\ &=& \sum_{k}M_{ik} + \sum_{k\neq j} M_{ik}T'_{kj} \nonumber \\ &=& 1 + \sum_{k} M_{ik}T'_{kj} - M_{ij}T'_{jj} \end{align}\] that is, \[T'=\mathbf{1}\mathbf{1}^T+M(T'-\mathrm{ddiag}(T')),\label{eq:cmt1}\tag{110}\] where \(\mathbf{1}\mathbf{1}^T\) denotes the all-ones matrix and \(\mathrm{ddiag}(T')\) is the diagonal matrix obtained by setting to zero all non-diagonal entries of \(T'\).

Let us try to determine \(\mathrm{ddiag}(T')\) first. To that end, let \({\psi}_{1}^T=[\deg(i)]_{1\leq i\leq n}\) be be the left eigenvector of \(M\) \[{\psi}_{1}^{T}M={\psi}_{1}^{T}.\] Multiply (110 ) by \({\psi}_{1}^{T}\) on both sides to obtain \[\begin{align} {\psi}_{1}^{T}T' & = & {\psi}_{1}^{T}\mathbf{1}\mathbf{1}^T+{\psi}_{1}^{T}M(T'-\mathrm{ddiag}(T'))\\ & = & {\psi}_{1}^{T}\mathbf{1}\mathbf{1}^T+{\psi}_{1}^{T}(T'-\mathrm{ddiag}(T')).\nonumber \end{align}\] Thus, we have \[{\psi}_{1}^{T}\mathrm{ddiag}(T')={\psi}_{1}^{T}\mathbf{1}\mathbf{1}^T={\psi}_{1}^{T}\mathbf{1}\mathbf{1}^T=\left(\sum_{k=1}^{n}\deg(k)\right)\mathbf{1}^T, \label{eq:cmt2}\tag{111}\] or equivalently, \[T'_{ii}=\frac{1}{\deg(i)}\sum_{k=1}^{n}\deg(k) = \frac{1}{\deg(i)} \operatorname{vol}(G),\] recall that te volume of \(G\) is the sum of the node degrees. \(T'_{ii}\) is the mean return (or recurrence) time from node \(i\) to itself. We see that it is proportional to \(\operatorname{vol}(G)\): The larger the graph, the more time the random walker wanders (on average) around before returning back. \(T'_{ii}\) is inversely proportional to \(\deg(i)\): The larger \(\deg(i)\), the more possibilities the random walker has to return to node \(i\). Remarkably, \(T'_{ii}\) does not see much about the structure of the graph, besides the total edge weights and the specific node degree.

Now we will show uniqueness of the solution \(T'\) to (110 ). Indeed, assume that \(T_2'\) is another solution such that \[T_2'=\mathbf{1}\mathbf{1}^T+M(T_2'-\mathrm{ddiag}(T')),\] then \[T'-T_2'=M(T'-T_2').\] We see that the columns on \(T'-T_2'\) are right eigenvectors of \(M\) with eigenvalue 1. But if \(M\) is irreducible, the only eigenvector with eigenvalue is 1 is the all-ones vector. Therefore, \[T'-T_2'={\boldsymbol{1}}{u}^{T}\] for some \(u\); the fact that \(T'\) and \(T'_2\) need to agree on the diagonal means that \(u\) must be zero and we must have \(T'_2=T'\).

Now that we established that \(T'\) is unique let us solve for it. Notice that that the mean first passage time matrix \[T = [t_{ij}]_{1\leq i,j \leq n}\] can be written as \[T = T' - \mathrm{ddiag}(T').\] Therefore, (110 ) is equivalent to \(T+\mathrm{ddiag}(T')=\mathbf{1}\mathbf{1}^T+MT\) which we can rewrite as \[(I-M)T=\mathbf{1}\mathbf{1}^T-\mathrm{ddiag}(T').\label{eq:cmt3}\tag{112}\] Unfortunately, \((I-M)\) is not invertible, so we cannot obtain \(T\) by simply inverting \((I-M)\).

Recall that the graph Laplacian is defined as \[L=D-W=D(I-M).\label{eq:graph95laplacian}\tag{113}\] Multiplying (112 ) by \(D\) gives \[LT=D\mathbf{1}\mathbf{1}^T-D\mathrm{ddiag}(T').\label{eq:cmt4}\tag{114}\] We also recall that \(L\mathbf{1}=0\) and \(L\succeq 0\).

The eigen-decomposition of \(L\) can be written as \[L=\sum_{l=1}^{n}\mu_{l}{\rho}_{l}{\rho}_{l}^{T} = \sum_{l=2}^{n}\mu_{l}{\rho}_{l}{\rho}_{l}^{T}\] with \(0=\mu_{1}\leq\mu_{2}\leq\cdots\leq\mu_{n}\), and \(\rho_1 = \frac{1}{\sqrt{n}}{\boldsymbol{1}}\). Note that since \(G\) is connected then \(\mu_{2}>0\). The pseudo-inverse of the graph Laplacian is given by \[L^\dagger=\sum_{l=2}^{n}\frac{1}{\mu_{l}}{\rho}_{l}{\rho}_{l}^{T}.\] Then, \(L^\dagger\) has the properties listed in 28 . Moreover, \[L^\dagger L=\sum_{l=2}^{n}{\rho}_{l}{\rho}_{l}^{T}=I-{\rho}_{1}{\rho}_{1}^{T} = I - \frac{1}{n}{\boldsymbol{1}}{\boldsymbol{1}}^T.\label{eq:cmt5}\tag{115}\]

Multiplying (114 ) by \(L^\dagger\) gives \[\left(L^\dagger L\right)T=L^\dagger D\mathbf{1}\mathbf{1}^T -L^\dagger DM_{d}.\label{eq:cmt6}\tag{116}\] Combining (115 ) and (116 ) for each column of \(T\), we have \[T=L^\dagger D\mathbf{1}\mathbf{1}^T-L^\dagger D\mathrm{ddiag}(T')+{\boldsymbol{1}}{ u}^{T}\] for some vector \(u\), that is, \[t_{ij}=\sum_{k=1}^{n}L_{ik}^\dagger \deg(k)-\operatorname{vol}(G)L_{ij}^\dagger +u_{j}\label{eq:cmt7}\tag{117}\] (here we used that \(D\,\mathrm{ddiag}(T')= \operatorname{vol}(G) I\)). From \(t_{ii}=0\) we find \[u_{i}=-\sum_{k=1}^{n}L_{ik}^\dagger \deg(k)+\text{vol}(G)L_{ii}^\dagger .\label{eq:cmt8}\tag{118}\] From (117 ) and (118 ), for \(i\neq j\), the average commute time between \({x}_{i}\) and \({x}_{j}\) is \[\begin{align} t_{ij}+t_{ji} & = & \sum_{k=1}^{n}L_{ik}^\dagger \deg(k)-\text{vol}(G)L_{ij}^\dagger -\sum_{k=1}^{n}L_{jk}^\dagger \deg(k)+\text{vol}(G)L_{jj}^\dagger \nonumber \\ & & +\sum_{k=1}^{n}L_{jk}^\dagger \deg(k)-\text{vol}(G)L_{ji}^\dagger -\sum_{k=1}^{n}L_{ik}^\dagger \deg(k)+\text{vol}(G)L_{ii}^\dagger \nonumber \\ & = & \text{vol}(G)\left(L_{ii}^\dagger +L_{jj}^\dagger -2L_{ij}^\dagger \right),\label{eq:cmt9} \end{align}\tag{119}\] where \[\begin{align} L_{ij}^\dagger & = & \sum_{l=2}^{n}\frac{1}{\mu_{l}}{\rho}_{l}(i){\rho}_{l}(j).\label{eq:cmt10} \end{align}\tag{120}\]

Define the commute time mapping \(\Psi : V \mapsto \mathbb{R}^{n-1}\) as \[\Psi(x_{i}) = \left(\frac{1}{\sqrt{\mu_{2}}}{\rho}_{2}(i),\,\frac{1}{\sqrt{\mu_{3}}}{\rho}_{3}(i),\,\cdots,\,\frac{1}{\sqrt{\mu_{n}}}{\rho}_{n}(i)\right)^{T}.\label{eq:commute95time95embedding}\tag{121}\] We see that the commute time between \(x_i\) and \(x_j\) is equivalent to the Euclidean distance between \(\Psi(x_i)\) and \(\Psi(x_j)\): \[t_{ij}+t_{ji}=\text{vol}(G)\left\Vert \Psi({x}_{i})-\Psi({x}_{j})\right\Vert ^{2}.\] Also, \[t_{ij}+t_{ji}=\text{vol}(G) (e_i - e_j)^T L^\dagger (e_i - e_j),\] where \(e_i\) is the \(i\)’th standard unit vector.

7.3.1 Comparison between the diffusion distance and the commute time distance↩︎

Let us assume that the graph is regular, that is, \(\deg(i)=d\) for all \(i=1,2,\ldots,n\), or equivalently \(D=dI\). Then, \[L = D - W = D(I-M) = d(I-M),\] so \(M\) and \(L\) share the same eigenvectors \[\phi_l = \rho_l,\] and the eigenvalues are related by \[\mu_l = d(1-\lambda_l).\]

The commute time distance can become meaningless for some families of graphs. We see that the filter corresponding to the commute time distance associates a similar non-zero weight to all small eigenvalues. If all eigenvalues are assigned the same weight (i.e, constant filter), then the resulting distance between nodes \(i\) and \(j\) is simply \(\frac{1}{d_i} + \frac{1}{d_j}\). This is quite boring, as it tells us nothing about the connectivity between the two nodes. It turns out that in certain cases the commute time distance has exactly this boring behavior [@von2014hitting].

We comment that many other filter/transfer functions are possible in practice, and it all depends on the problem at hand.

7.3.2 Effective resistance↩︎

The commute time distance is also proportional to the resistance distance, where the graph is viewed as an electrical network of resistors, with the resistance of each edge given by \[r_{ij} = w_{ij}^{-1}\] Using Ohm’s Law and Kirckhoff’s Law it is possible to compute the effective resistance between any two nodes in the graph. The effective resistance turns out to be exactly the same as the commute time distance just without the factor of \(\text{vol}(G)\). This provides an alternative physical interpretation of the commute time distance.

To see this, suppose that \(\mathbf{v} \in \mathbb{R}^n\) is the vector of electric potentials at the vertices. Recall Ohm’s Law that relates the potential with the resistance \(R\) and the electric current \(I\) \[V = RI\] Let \(\mathbf{i} \in \mathbb{R}^{2|E|}\) be the vector of signed edge currents, with \[\mathbf{i}_{i \to j} = \frac{\mathbf{v}_i - \mathbf{v}_j}{r_{ij}} = w_{ij}(\mathbf{v}_i - \mathbf{v}_j)\] We see that \[\mathbf{i} = \Delta U\mathbf{v},\] where \(U\) is an \(|E|\times n\) signed edge adjacency matrix of \(G\), meaning that the row of \(U\) corresponding to \((i,j)\in E\) is given by \(e_i-e_j\). The signing is arbitrary but needs to be set. \(\Delta=\mathop{\mathrm{diag}}(\mathrm{vec}[w_{ij}]_{(i,j)\in E})\) is the \(|E|\times |E|\) matrix whose diagonal elements are the edge weights.

Also, Kirckhoff’s Law says that a resistor network cannot hold current: all flow entering vertex \(i\) from edges in the graph must exit \(i\) to an external source. Let \(\mathbf{i}_{\text{ext}} \in \mathbb{R}^n\) denote the vector of external currents, where \(\mathbf{i}_{\text{ext}}(i)\) is the amount of current entering the graph through node \(i\). We then have \[\mathbf{i}_{\text{ext}}(i) = \sum_{j : (i,j)\in E} \mathbf{i}_{i\to j}.\] Equivalently, in matrix form \[\mathbf{i}_{\text{ext}} = U^T \mathbf{i} = U^T \Delta U\mathbf{v} = L\mathbf{v}\] Thus, by applying the Laplacian on the vector of potentials we get the external source current. In particular, we distinguish between internal nodes \(i\) for which \(\mathbf{i}_{\text{ext}}(i)=0\), and boundary nodes with \(\mathbf{i}_{\text{ext}}(i) \neq 0\). Notice that for internal nodes \[L \mathbf{v}(i)=0\] is equivalent to \[\mathbf{v}_i = \frac{1}{\deg(i)} \sum_{j:(i,j)\in E} w_{ij} \mathbf{v}_j\] that is, the potential at the node is the weighted average of the potential at its neighbors.

Now, in order to calculate the effective resistance between nodes \(i\) and \(j\) we need to solve \[L\mathbf{v} = \mathbf{i}_{\text{ext}}\] with nodes \(i\) and \(j\) being boundary nodes and all other nodes being internal nodes. That is, we let \[\mathbf{i}_{i,j,\text{ext}}(a) = I (e_i - e_j) = \left\{\begin{array}{rc} I & a=i \\ -I & a=j \\ 0 & \text{otherwise} \end{array} \right.\] Notice that \(\mathbf{1}^T \mathbf{i}_{i,j,\text{ext}}=0\), therefore if the graph is connected (the eigenvalue 0 is simple), the solution for the potential vector \(\mathbf{v}\) is given by \[\mathbf{v} = L^\dagger \mathbf{i}_{i,j,\text{ext}} = I L^\dagger(e_i - e_j)\] The potential drop between nodes \(i\) and \(j\) is \[\mathbf{v}_i - \mathbf{v}_j = I (e_i - e_j)^T L^\dagger (e_i - e_j)\] Which shows that the effective resistance is \[R^{\text{eff}}(i,j) = \frac{\mathbf{v}_i - \mathbf{v}_j}{I} = (e_i - e_j)^T L^\dagger(e_i - e_j)\] Therefore, the effective resistance and the commute time distance are equivalent (up to the factor of \(\text{vol}(G)\)).

7.3.3 Time synchronization on graphs↩︎

Here we demonstrate another application of the pseudo-inverse of the Laplacian and the commute time distance / effective resistance. Consider the time-synchronization problem, in which \(n\) sensing nodes are part of a multi-hop communication network. Each node has a local clock, but each pair of clocks differs by a constant offset. However, nodes that communicate directly can estimate the difference between their local clocks by exchanging messages that are time stamped with local clock times. For example, suppose that nodes \(i\) and \(j\) can communicate directly with each other and have clock offsets \(t_i\) and \(t_j\) with respect to a reference clock. By passing messages back and forth, the nodes can measure the relative clock offset \(t_i-t_j\) with the real-valued measurement \[\label{offsets} s_{ij} = t_i - t_j + \xi_{ij},\qquad (i,j)\in E\tag{122}\] where \(\xi_{ij}\) is the measurement noise. The task is to estimate the clock offsets \(t_1,\ldots, t_n\) (up to a global offset) from such pair-wise measurements. It is also an instance of a ranking problem, in which \(t_i\) may represent the “strength" of the \(i\)’th player, and \(s_{ij}\) is the outcome or score of a match between players \(i\) and \(j\). Clearly, we need that the graph of measurements is connected. Also, suppose that the random variables \((\xi_{ij})_{(i,j)\in E}\) are independent zero-mean Gaussians, i.e., \[\xi_{ij} \sim \mathcal{N}(0,\sigma_{ij}^2),\qquad (i,j)\in E\] Clearly, using the offset measurements we can only hope to estimate \(t_1,\ldots,t_n\) up to a global shift. We may therefore assume that they are centered, i.e., \(\sum_{i=1}^n t_i = 0\), or \[{\mathbf{1}}^T t = 0\] We seek to find a solution to the over-determined system of linear equations (122 ) which in vector notation can be written as \[s = Ut + \xi\] and \(\xi \sim \mathcal{N}(0,\mathop{\mathrm{diag}}(\sigma)^2)\), where \(\mathop{\mathrm{diag}}(\sigma)\) is the diagonal \(|E|\times |E|\) matrix with \(\sigma_{ij}\) in the diagonal. We pre-whiten the noise by dividing each equation with the standard deviation \(\sigma_{ij}\) to get \[\mathop{\mathrm{diag}}(\sigma)^{-1} s = \mathop{\mathrm{diag}}(\sigma)^{-1} Ut + z\] where \[z \sim \mathcal{N}(0,I)\] The least squares estimator \(\hat{t}\) satisfies \[U^T \mathop{\mathrm{diag}}(\sigma)^{-2} U\hat{t} = U^T \mathop{\mathrm{diag}}(\sigma)^{-2} s\] We therefore define weights \[w_{ij} = \sigma_{ij}^{-2}\] for which the least squares estimator satisfies \[L \hat{t} = U^T \Delta s\], where \(\Delta = \mathop{\mathrm{diag}}(\mathrm{vec}[w_{ij}]_{(i,j)\in E})\) is the \(|E|\times |E|\) matrix whose diagonal elements are the edge-weights. The solution (assuming \(\mathbf{1}^T t = 0\)) is then given by \[\hat{t} = L^\dagger U^T \Delta s\] Notice that the least squares estimator is unbiased: \[\mathbb{E}[\hat{t}] = L^\dagger U^T \Delta \mathbb{E}[s] = L^\dagger U^T \Delta Ut = L^\dagger L t = t\] The error is given by \[\hat{t} - t = L^\dagger U^T \Delta (s - \mathbb{E}[s]) = L^\dagger U^T \Delta \xi\] Hence, the covariance of the estimation error is given by \[\mathbb{E}[(\hat{t} - t)(\hat{t}-t)^T] = L^\dagger U^T \Delta \Delta^{-1} \Delta U L^\dagger = L^\dagger L L^\dagger = L^\dagger\] In particular, the mean squared error (MSE) of the \(i\)’th clock is the \(i\)’th diagonal term of \(L^\dagger\) \[\mathbb{E}[(\hat{t}_i - t_i)^2] = L^\dagger_{ii},\] and the MSE of all \(n\) estimators is \[\mathbb{E}[\|\hat{t} - t\|^2] = \sum_{i=1}^n L^\dagger_{ii} = \operatorname{Tr}(L^\dagger).\] Finally, the offset \(t_i-t_j\) between nodes \(i\) and \(j\) is estimated as \(\hat{t}_i - \hat{t}_j\). The MSE of the offset is \[\begin{align} \mathbb{E}[((\hat{t}_i - \hat{t}_j) - (t_i - t_j))^2] &=& \mathbb{E}[((\hat{t}_i - t_i) - (\hat{t}_j - t_j))^2] \nonumber \\ &=& \mathbb{E}[(\hat{t}_i - t_i)^2] + \mathbb{E}[(\hat{t}_j - t_j)^2] - 2\mathbb{E}[(\hat{t}_i - t_i)(\hat{t}_j - t_j)] \nonumber \\ &=& L^\dagger_{ii} + L^\dagger_{jj} - 2L^\dagger_{ij}, \end{align}\] which is the effective resistance distance. More accurate offset estimates are reflected by shorter commute time distances.

7.4 Other non-linear dimensional reduction techniques↩︎

There are several other similar non-linear dimensional reduction methods. A particularly popular one is Isomap [@ISOMAP_paper_2000]. The idea is to find an embedding in \(\mathbb{R}^d\) for which Euclidean distances in the embedding correspond as much as possible to geodesic distances in the graph. This can be achieved by, between pairs of nodes \(v_i,\,v_j\) finding their geodesic distance and then using, for example, Multidimensional Scaling to find points \(y_i\in\mathbb{R}^d\) that minimize (for example) \[\min_{y_1,\dots,y_n \in \mathbb{R}^d}\sum_{i,j} \left( \|y_i-y_j\|^2 - \delta_{ij}^2 \right)^2,\] which can be done with spectral methods (it is a good exercise to compute the optimal solution to the above optimization problem).

Locally linear embedding (LLE) [@roweis2000nonlinear] and Isomap are often considered landmarks in the area of nonlinear dimension reduction. Beyond diffusion maps [@RRCoifman_SLafon_2006] and Laplacian eigenmaps [@belkin2001laplacian; @belkin2003laplacian], there are several other nonlinear dimension reduction methods worth mentioning here, including Hessian LLE [@donoho2003hessian], local tangent space alignment [@zhang2004principal], and UMAP [@mcinnes2018umap]. There are other variants of diffusion maps such as anisotropic diffusion maps based on the Mahalanobis distance [@singer2008non] and alternating diffusion [@lederman2018learning] for learning the geometry of common latent variables. Here, we choose to present in more detail t-Distributed Stochastic Neighbor Embedding (t-SNE) [@maaten2008visualizing], a popular method for visualization of high-dimensional data, and later in Chapter 12.4 vector diffusion maps [@ASinger_HTWu_2011_VDM] that uses the graph connection Laplacian and generalizes diffusion maps from scalar functions to vector fields.

7.4.1 t-Distributed Stochastic Neighbor Embedding↩︎

t-SNE is a nonlinear dimension reduction method designed primarily for the visualization of high-dimensional data [@maaten2008visualizing]. Unlike spectral methods such as diffusion maps, which are based on eigenvalue problems and have a clear operator-theoretic interpretation, t-SNE is formulated as a stochastic optimization problem whose objective is to preserve local neighborhood relations in a probabilistic sense.

The central idea of t-SNE is that points that are close in high-dimensional space should remain close in the embedding, while distant points are modeled with much weaker constraints. To that end we will define two distances, one suitable for the original high-dimensional data space, and one for the low-dimensional embedding space.

Given a set of points \(\{x_1, \dots, x_n\} \subset \mathbb{R}^d\), we define the similarity of point \(x_j\) to \(x_i\) as the conditional probability \(p_{i|j}\) given by \[p_{i|j} = \frac{\exp(-\|x_i - x_j\|^2 / 2\sigma_i^2)}{\sum_{k \neq i} \exp(-\|x_i - x_k\|^2 / 2\sigma_i^2)} \quad \text{and} \quad p_{ij} = \frac{p_{j|i} + p_{i|j}}{2n}.\] In contrast to graph Laplacians and diffusion maps, the bandwidth \(\sigma_i\) is not global. Instead it is determined such that it is adapting the scale to the local density of the data.

Let \(\{y_1, \dots, y_n\} \subset \mathbb{R}^m\) (where \(m \ll d\), typically \(m=2\) or \(3\)) denote the embedded points. In the low-dimensional space, we use a heavy-tailed Student t-distribution with one degree of freedom (which is a Cauchy distribution) to model similarities \[q_{ij} = \frac{(1 + \|y_i - y_j\|^2)^{-1}}{\sum_{k \neq l} (1 + \|y_k - y_l\|^2)^{-1}}.\] The choice of the t-distribution in the embedding space alleviates the so-called crowding problem, which arises when high-dimensional neighborhoods are projected into low-dimensional spaces. The heavy tails of the t-distribution allow distant points in \(\mathbb{R}^d\) to be modeled by even larger distances in \(\mathbb{R}^m\), preventing the “crowding” at the center of the map often seen in standard SNE [@hinton2002stochastic].

The embedding is determined by minimizing the Kullback-Leibler (KL) divergence between the joint probability distribution \(P\) (in \(\mathbb{R}^d\)) and \(Q\) (in \(\mathbb{R}^m\)) \[\label{KL-tsne} C = \text{KL}(P \| Q) = \sum_{i \neq j} p_{ij} \log \frac{p_{ij}}{q_{ij}}\tag{123}\] with respect to the variables \(y_1,\dots,y_n\). This cost function is non-convex and is typically optimized via gradient descent³². The gradient has a particularly intuitive “spring-force” interpretation \[\frac{\partial C}{\partial y_i} = 4 \sum_j (p_{ij} - q_{ij})\frac{y_i - y_j}{1 + \|y_i - y_j\|^2}\] with two opposing forces:

• Attraction: The term \(p_{ij}\) pulls \(y_i\) toward \(y_j\) if they are close in the original space.

• Repulsion: The term \(q_{ij}\) pushes \(y_i\) away from \(y_j\) if they are close in the embedding space but not in the original space.

We note that the embedding provided by t-SNE is invariant under global rotations and translations but not under scaling, and global distances and large-scale geometry are generally distorted. Thus, t-SNE may therefore better be viewed as a local neighborhood-preserving visualization method rather than a general-purpose dimension reduction technique.

We demonstrate the effectiveness of t-SNE by applying it to the MNIST (Modified National Institute of Standards and Technology) database [@lecun1998mnist] of handwritten digits. Each datapoint (image) \(x_i\) is an \(28 \times 28\) black-and-white image containing a handwritten digit. See the left panel of Figure 34 for some example images. The corresponding 2-D embedded points \(y_i\) are computed via t-SNE. Since the t-SNE objective is non-convex, the final embedding is sensitive to initialization and hyperparameters. Here, we use the implementation proposed in [@linderman2019clustering]. The resulting embedding is shown in the right panel of Figure 34. It is evident that t-SNE is quite effective in separating the different digits.

There are some interesting connections between t-SNE and spectral clustering, as uncovered in [@linderman2019clustering] and [@cai2022theoretical], which provides some of the first rigorous mathematical justifications for why t-SNE behaves so effectively as a clustering tool. First note that the attractive term resembles a graph Laplacian \((L_P y)_i = \sum_j p_{ij}(y_i - y_j)\), where \(P\) acts as a weighted adjacency matrix. Also, the gradient descent flow during the initial phase behaves similarly to power iteration on a graph Laplacian [@cai2022theoretical]. Specifically, the early iterations essentially perform a form of spectral clustering. The clusters are formed based on the dominant eigenvectors of the transition matrix, much like the initial stages of a diffusion process. It is furthermore shown in [@linderman2019clustering] that if the data contains well-separated clusters (quantified by the spectral gap of the Laplacian), t-SNE is provably guaranteed to keep these clusters separated in the low-dimensional embedding. These connections should not be interpreted as placing t-SNE within the framework of spectral dimension reduction methods; rather, they provide insight into the local dynamics of its optimization landscape.

Diffusion maps are sometimes used as a preprocessing step for t-SNE [@wolf2018scanpy; @moon2019visualizing; @haghverdi2015diffusion]. The standard pipeline, common in computational biology [@wolf2018scanpy], is to first reduce the data to a moderate number of diffusion coordinates which removes noise and compresses the intrinsic geometric structure into a tractable representation, and then to apply t-SNE (or the closely related UMAP algorithm [@mcinnes2018umap]) to this \(k\)-dimensional intermediate representation to produce the final 2D visualization. This two-stage approach aims to combine the strengths of both methods: diffusion maps provide a geometrically meaningful and noise-robust intermediate representation, while t-SNE optimizes that representation specifically for visual interpretability in two dimensions.

7.5 Kernel learning↩︎

Linear methods are computationally efficient and well-understood, but they fail when the relationship between variables is inherently non-linear. The kernel approach allows us to lift our data into a high-dimensional (potentially infinite-dimensional) feature space where linear boundaries can capture complex structures, all without ever explicitly computing the coordinates in that space, a procedure referred to as the kernel trick.

One of the primary advantages of the kernel approach is its ability to bypass the often difficult task of explicit feature engineering. In many domains—network analysis being a prime example—practitioners may struggle to define a concrete feature vector for an object, yet they can intuitively provide a measure of affinity or similarity between pairs of objects.

From a mathematical perspective, this “user-defined” similarity measure can be accompanied with some rigorous guarantees, as we will demonstrate in this section: If the practitioner’s affinity function qualifies as a positive semi-definite kernel (a property verifiable via Bochner’s Theorem, Theorem 7.6), then Mercer’s Theorem (Theorem 7.5) guarantees that this measure corresponds to a valid inner product in some high-dimensional feature space. Consequently, the kernel trick allows us to perform sophisticated non-linear learning via a similarity metric, rather than requiring the explicit construction of a feature map.

7.5.1 Positive semidefinite kernels and reproducing kernels↩︎

In practice, choosing a kernel is a combination of principled heuristics and empirical validation. The linear kernel \(K(x,x') = x^\top x'\) assumes the data is already linearly separable (or nearly so) in its original coordinates. The polynomial kernel \(K(x,x') = (\langle x, x' \rangle + r)^d\), is useful in image processing where pixel interactions matter. The Gaussian kernel \(K(x,x') = \exp (-\|x-x'\|^2/2\varepsilon^2)\), which is sometimes also called RBF kernel (where RBF stands for Radial Basis Function) serves as a powerful default for most general-purpose applications, provided its hyperparameter is tuned carefully. The Laplacian kernel \(K(x,x') = \exp\bigl(-\gamma\|x-x'\|\bigr)\) for \(\gamma>0\) is well suited for modeling physical processes or signals that have sharp changes or disjointed behavior (e.g., “edges” in natural images) that a Gaussian kernel might over-smooth.

One property that turns out to be very useful for a kernel is positive-definiteness.

Condition 124 is purely algebraic and requires no topological structure on \(\mathcal{X}\). The set \(\mathcal{X}\) can be discrete, a manifold, a function space, or any abstract set.

Reproducing Kernel Hilbert Spaces (RKHS) provide a convenient and rigorous functional analysis framework for kernel learning.

The following proposition (whose proof is left to the reader) summarizes basic properties of the reproducing kernel.

A kernel \(K\) can often be understood via an explicit feature map \(\phi: \mathcal{X} \to \mathcal{F}\), where \(\mathcal{F}\) is some (possibly high-dimensional) feature space, such that \[\label{eq:kerneltrick} K(x, x') = \langle \phi(x),\, \phi(x') \rangle_{\mathcal{F}}.\tag{126}\] This means the kernel computes the inner product between the feature representations of \(x\) and \(x'\), without needing to construct \(\phi\) explicitly. This is the kernel trick: algorithms that depend only on inner products between inputs can be kernelized by replacing \(\langle x, x' \rangle\) with \(K(x, x')\), implicitly working in \(\mathcal{F}\).

As a simple example, consider \(x \in \mathbb{R}^2\) and the polynomial kernel given by \(K(x, x') = (x^\top x')^2\). Expanding this out reveals it equals \(\langle \phi(x), \phi(x') \rangle\) for \(\phi(x) = (x^2,\, (x')^2,\, \sqrt{2} x x')^\top\), a three-dimensional feature map computed via a scalar operation.

7.5.2 Mercer’s Theorem↩︎

While the kernel trick allows us to compute similarities without knowing the coordinates of the data in the feature space, Mercer’s Theorem provides the explicit mathematical construction of those coordinates. It provides a constructive answer to the question “which functions \(K\) are valid kernels, i.e., correspond to an inner product in some feature space?”

The closely related but more general Moore-Aronszajn Theorem [@aronszajn1950theory] states that every positive definite kernel defines a unique RKHS (existence and structure). In comparison, Mercer’s Theorem establishes under additional analytic assumptions, that the kernel admits an eigenfunction expansion which comes with an explicit feature map given by \[\phi(x) = \bigl(\sqrt{\lambda_1}\psi_1(x),\,\sqrt{\lambda_2}\psi_2(x),\,\ldots\bigr).\] While this concrete construction of the feature map is useful for theoretical purposes, in practice we avoid computing it by using the kernel trick.

7.5.3 Bochner’s Theorem↩︎

While Mercer’s Theorem provides a link between a PSD kernel and a feature map in \(L_2\), Bochner’s Theorem provides a powerful characterization of PSD kernels for the case when the kernel is translation-invariant, i.e., a kernel that can be expressed in the form \(K(x, x') = \kappa(x - x')\) for some function \(\kappa\).

Consider a finite Borel measure \(\mu\) of \(\mathbb{R}^p\). That is, \(\mu(\mathbb{R}^p) < \infty\) and \(\mu(A) = \int_{A \subset \mathbb{R}^p} \,d\mu\). For example, if \(\mu\) is absolutely continuous with respect to the Lebesgue measure \(dx\) then \(d\mu = f(x)\,dx\) where \(f\) is the probability density function. The Fourier transform of \(\mu\) is \[\hat{\mu}(\xi) = \int_{\mathbb{R}^p} e^{-2\pi \imath \langle \xi, x \rangle}\,d\mu(x).\] Observe that \[|\hat{\mu}(\xi)| \leq \mu(\mathbb{R}^p) < \infty.\] That is, \(\hat{\mu}\) is bounded. Also, \(\hat{\mu}\) is a continuous function. Moreover, for every \(n=1,2,\ldots\) we can associate with \(\hat{\mu}\) a PSD quadratic form: \[\begin{align} \sum_{i,j=1}^n \hat{\mu}(\xi_i-\xi_j) \alpha_i \bar{\alpha}_j &=& \sum_{i,j=1}^n \alpha_i \bar{\alpha}_j \int_{\mathbb{R}^p} e^{-2\pi \imath \langle \xi_i - \xi_j, x\rangle}\,d\mu(x) \\ &=& \int_{\mathbb{R}^p} \left| \sum_{i=1}^n \alpha_i e^{-2\pi \imath \langle \xi_i, x \rangle }\right|^2 \,d\mu(x) \geq 0. \end{align}\] We conclude that \(\hat{\mu}\) is bounded, continuous, and the matrix \((\hat{\mu}(\xi_i-\xi_j))_{i,j=1}^n\) is PSD for every \(n=1,2,\ldots\) and for all \(\xi_1,\ldots,\xi_n \in \mathbb{R}^p\).

But the converse is also true: If \(f\) is bounded, continuous, and the matrix \((f(x_i-x_j))_{i,j=1}^n\) is PSD for all choices of \(n=1,2,\ldots\) and \(x_1,\ldots,x_n\in \mathbb{R}^p\), then \(f\) is the Fourier transform of a finite Borel measure.

From the convolution theorem and Parseval’s theorem we get \[\langle \hat{K} \hat{\alpha}, \hat{\alpha} \rangle \geq 0.\] It follows that \[\int_{\mathbb{R}^p} |\hat{\alpha}|^2(\xi) \hat{K}(\xi)\,d\xi \geq 0.\] Since we can choose \(\alpha(x)\) (or \(\hat{\alpha}(\xi)\)), we must have \(\hat{K}(\xi) \geq 0\).

We thus have proved Bochner’s theorem [@bochner1932vorlesungen]:

When \(K\) is normalized so that \(K(0) = 1\), the measure \(\mu(\xi)\) becomes a probability distribution over frequencies \(\xi \in \mathbb{R}^d\). For the Gaussian kernel \[K(x, x') = \exp\!\left(-\frac{\|x - x'\|^2}{2\varepsilon^2}\right),\] the spectral measure \(\mu(\xi)\) is a Gaussian distribution with variance \(1/\varepsilon^2\) (because the inverse Fourier transform of the Gaussian is also a Gaussian, see e.g. [@Gro01]), which immediately confirms that the Gaussian kernel is PSD.

This has immediate consequences for diffusion maps. Let \(x_1,x_2,\ldots, x_n \in \mathbb{R}^p\) and consider the \(n\times n\) weight matrix \(W\) with entries \(W_{ij} = e^{-\frac{\|x_i - x_j \|^2}{2\varepsilon^2}}\). Then \(W\) is PSD for all values of \(n\), \(p\), and \(\varepsilon\). More generally, Bochner’s theorem implies that for positive kernels the eigenvalues of \(W\) and the adjacency matrix \(A\) are non-negative and we can use the diffusion map with any positive real \(t>0\).

7.5.4 The Representer Theorem↩︎

The following theorem, known as Representer Theorem, is fundamental in kernel methods because it guarantees that the optimal solution to a regularized cost functional lives in a finite-dimensional subspace, regardless of how large the underlying Hilbert space is.

The representer theorem acts as a kind of dimensionality reduction principle for optimization: although learning occurs in a potentially infinite-dimensional RKHS, the optimal solution lies in the finite-dimensional span of kernel evaluations at the training data. Consequently, the optimization problem depends on the number of training observations rather than the dimensionality of the feature space.

7.5.5 Applications of kernel learning↩︎

Kernel ridge regression: The kernelized version of ridge regression replaces \(\min \|X\theta - y\|^2 + \lambda\|\theta\|^2\) with a minimization over functions \(f\) in a RKHS \(\mathbb{H}\), i.e., we consider \[\min_{f \in \mathbb{H}} \,\frac{1}{m}\sum_{i=1}^m (f(x_i) - y_i)^2 + \lambda \|f\|^2_{\mathbb{H}}.\] By Theorem 7.7, the solution takes the form \(f^*(x) = \sum_{i=1}^m \alpha_i K(x_i, x)\), and substituting this into the objective reduces the problem to \[\min_{\alpha \in \mathbb{R}^m} \,\frac{1}{m}\|K\;\alpha - y\|^2 + \lambda\, \alpha^\top K \alpha,\] where \(G_K\) is the \(m \times m\) Gram matrix with entries \((G_K)_{ij} = K(x_i, x_j)\). Setting the gradient to zero yields the closed-form solution: \[\alpha = (K + \lambda m I)^{-1}y.\] Predictions at a new point \(x\) are then given by \(f^*(x) = K(x)^\top \alpha\), where \(K(x) = (K(x_1, x), \dots, K(x_m, x))^\top\). The key insight is that the algorithm never requires explicit access to the feature map \(\phi(x)\); all computations are performed through the kernel matrix \(G_K\).

Kernel regression: While kernel ridge regression is an optimization-based approach, kernel regression is a localized density estimation approach. In particular, it is a non-parametric regression method that estimates the relationship between an input variable \(x\) and an output \(y\) by taking a weighted average of observed responses, where the weights depend on the similarity between inputs measured by a kernel function.

Suppose we observe data \((x_1,y_1),\dots,(x_n,y_n)\), to predict the value of \(y\) at a new point \(x\), kernel regression assigns larger weights to observations whose inputs are close to \(x\). The resulting prediction is a weighted average of the form \[\hat{f}(x) = \frac{\sum_{i=1}^{n} K(x,x_i) y_i}{\sum_{i=1}^{n} K(x,x_i)}.\] This formula is known as the Nadaraya–Watson kernel regression estimator [@nadaraya1964estimating; @watson1964smooth].

See [@gilles2025cryo] for an application of kernel regression in cryogenic electron microscopy (cryo-EM).

Kernel PCA: Standard PCA finds directions of maximum variance by computing the eigenvectors of the sample covariance matrix. Kernel PCA [@scholkopf1998nonlinear] extends this to nonlinear settings by implicitly performing PCA in the feature space \(\mathcal{F}\). By exploiting the kernel trick, kernel PCA does not need to set up a feature map explicitly. The (centered) covariance operator in \(\mathcal{F}\) is: \[C = \frac{1}{m}\sum_{i=1}^m \langle \phi(x_i), \phi(x_i)\rangle,\] and one seeks eigenvectors \(v \in \mathcal{F}\) satisfying \(Cv = \mu v\). Since any such eigenvector must lie in the span of \(\{\phi(x_i)\}\), we can write \(v = \sum_{i=1}^m \alpha_i \phi(x_i)\), and the eigenvalue problem reduces to: \[\tilde{G}\alpha = m\mu\, \alpha,\] where \(\tilde{G}\) is the centered Gram matrix \[\tilde{G} = G_K - \mathbf{1}_m G_K - G_K\mathbf{1}_m + \mathbf{1}_m G_K \mathbf{1}_m,\] with \(\mathbf{1}_m = \frac{1}{m}\mathbf{1}\mathbf{1}^\top\) the \(m\times m\) matrix with all entries \(1/m\). The projection of a new point \(x\) onto the \(r\)-th kernel principal component is: \[z_r(x) = \sum_{i=1}^m \alpha_i^{(r)}\, K(x_i, x),\] where \(\alpha^{(r)}\) is the \(r\)-th eigenvector of \(\tilde{G}\). Kernel PCA is particularly effective for dimensionality reduction and feature extraction when the data lies on a nonlinear manifold, as the feature space \(\mathcal{F}\) may “unfold” structure that is invisible to linear PCA in the input space.

Both kernel PCA and diffusion maps are spectral manifold learning techniques that embed high-dimensional data into a lower-dimensional space by solving an eigenvalue problem derived from a kernel matrix \(K(x,x')\). But kernel PCA operates on a centered version of \(K\). It seeks the directions of maximum variance in a high-dimensional feature space, hence it is concerned with certain geometric projections in Hilbert space. It is also highly sensitive to data density. In contrast, diffusion maps operate on a normalized version of \(K\), so it can be made more invariant to density via renormalization. It views \(K\) as a transition matrix for a Markov chain on the data and emphasizes intrinsic geometry and connectivity, rather than variance in feature space.

Random Fourier Features: Assume we are given a dataset \(\{x_1, \dots , x_n\} \in \mathbb{R}^d\). Since the measure \(\mu(\xi)\) arising in Bochner’s Theorem is a probability distribution (for a properly normalized kernel \(K\)), we can write \[K(x, x') = \mathbb{E}_{\xi \sim \mu}\!\left[e^{2\pi \imath \xi^\top(x - x')}\right] = \mathbb{E}_{\omega \sim \mu}\!\left[\cos\!\left(\omega^\top(x - x')\right)\right],\] where the second equality holds when \(\mu\) is symmetric. This means we can approximate \(K(x, x')\) by drawing random frequencies \(\xi_1, \dots, \xi_r \sim \mu\) with \(r \ll n\) and constructing the \(r\)-dimensional random Fourier feature map \[\phi(x) = \sqrt{\frac{2}{d}}\Bigl[\cos(\xi_1^\top x + b_1),\; \dots,\; \cos(\xi_r^\top x + b_r)\Bigr]^\top,\] where \(b_i \sim \text{Uniform}[0, 2\pi]\), so that \(K(x, x') \approx \langle \phi(x), \phi(x') \rangle\). This approximation (due to Rahimi & Recht [@rahimi2007random]) allows kernel methods to scale to large datasets by replacing expensive \(\mathcal{O}(n^2)\) kernel matrix computations with cheaper inner products in a \(r\)-dimensional feature space.

The Random Fourier Features method can be used for instance to create a low-rank approximation of the kernel matrix for kernel PCA. This technique makes kernel PCA scalable to large datasets. Instead of working with the kernel matrix \(K \in \mathbb{R}^{n \times n}\), we work with the explicit feature matrix \(\Phi \in \mathbb{R}^{n \times D}\): \[\label{eq:feature95matrix} \Phi = \begin{bmatrix} \phi(x_1)^\top\\ \phi(x_2)^\top\\ \vdots\\ \phi(x_n)^\top \end{bmatrix} \in \mathbb{R}^{n \times D}\tag{128}\] Then one can show that with a certain probability \(K \approx \Phi\Phi^\top\), see [@rahimi2007random]. This is a form of matrix sketching, see also Chapter 8.3.

Kernel SVM: Another natural use case of the kernel approach arises in connection with Support Vector Machines, a classification method which will be discussed in detail (together with its kernel variant) in Chapter 10.4.

7.6 Semi-supervised learning↩︎

Classification is a central task in machine learning which we will explore in detail in Chapters 10 and 11. In this section, we discuss an application of graph Laplacians and diffusion maps to classification in the semi-supervised setting.

In a supervised learning setting we are given many labeled examples and want to use them to infer the label of a new, unlabeled example. For simplicity, let us focus on the case of two labels, \(\{-1,+1\}\).

Let us consider the task of labeling the point “?” in Figure 35 given the labeled points. The natural label to give to the unlabeled point would be \(1\).

However, if we are given not just one unlabeled point, but many, as in Figure 36; then it starts being apparent that \(-1\) is a more reasonable guess.

Intuitively, the unlabeled data points allowed us to better learn the intrinsic geometry of the dataset. That is the idea behind Semi-Supervised Learning (SSL). While supervised learning relies exclusively on labeled examples, semi-supervised learning combines labeled and unlabeled data to leverage the underlying structure of the input space when labeled data are limited.

So, let us explore how we can exploit the structure of the input data to improve classification when labels are scarce. Just as above, we will use the data points to construct (via a kernel \(K_\varepsilon\)) a graph \(G=(V,E,W)\) where nodes correspond to points. More precisely, let \(l\) denote the number of labeled points with labels \(f_1,\dots,f_l\), and \(u\) the number of unlabeled points (with \(n=l+u\)), the first \(l\) nodes \(v_1,\dots,v_l\) correspond to labeled points and the rest \(v_{l+1},\dots,v_n\) are unlabaled. We want to find a function \(f:V\to \{-1,1\}\) that agrees on labeled points: \(f(i)=f_i\) for \(i=1,\dots,l\) and that is “as smooth as possible” on the graph. A way to pose this is the following [@zhu2003semi] \[\min_{f:V\to \{-1,1\}:\, f(i)=f_i\, i=1,\dots,l} \sum_{i<j}w_{ij}\left( f(i) - f(j) \right)^2.\] Instead of restricting ourselves to giving \(\{-1,1\}\) we allow ourselves to give real valued labels, with the intuition that we can “round” later by, e.g., assigning the sign of \(f(i)\) to node \(i\).

We thus are interested in solving \[\min_{f:V\to \mathbb{R}:\, f(i)=f_i\, i=1,\dots,l} \sum_{i<j}w_{ij}\left( f(i) - f(j) \right)^2.\]

If we denote by \(f\) the vector (in \(\mathbb{R}^n\) with the function values) then, recalling Proposition 6.2, we are can rewrite the problem as \[\begin{align} \sum_{i<j}w_{ij}\left( f(i) - f(j) \right)^2 = f^TL_Gf. \end{align}\]

Let us consider our problem \[\min_{f:V\to \mathbb{R}:\, f(i)=f_i\, i=1,\dots,l} f^TL_G f.\] We can write \[D = \left[ \begin{array}{cc} D_L & 0 \\ 0 & D_U \end{array} \right], \quad W = \left[ \begin{array}{cc} W_{LL} & W_{LU} \\ W_{UL} & W_{UU} \end{array} \right], \quad L_G = \left[ \begin{array}{cc} D_L - W_{LL} & -W_{LU} \\ -W_{UL} & D_U - W_{UU} \end{array} \right],\] and \(f = \left[ \begin{array}{c} f_L \\ f_U \end{array} \right]\), where \(L\) and \(U\) correspond to the set of indices of labeled and unlabeled nodes, respectively (not to be confused with the graph Laplacian \(L_G\)).

Then we want to find (recall that \(W_{UL}=W_{LU}^T\)) \[\min_{f_U\in \mathbb{R}^{u}} f_L^T\left[ D_L - W_{LL}\right] f_L - 2f_U^TW_{UL}f_L + f_U^T\left[ D_U - W_{UU}\right] f_U.\] By first-order optimality conditions³³, it is easy to see that the optimal satisfies \[\left( D_U - W_{UU}\right) f_U = W_{UL}f_L.\] If \(D_U-W_{UU}\) is invertible³⁴, then \[f_U^\ast = \left( D_U - W_{UU}\right)^{-1}W_{UL}f_L.\]

7.6.0.1 An interesting experiment and the Sobolev Embedding Theorem

Let us try a simple experiment. Let us say we have a grid on \([-1,1]^d\) dimensions (with say \(m^d\) points for some large \(m\)) and we label the center as \(+1\) and every node that is at distance larger or equal to \(1\) to the center, as \(-1\). We are interested in understanding how the above algorithm will label the remaining points, hoping that it will assign small numbers to points far away from the center (and close to the boundary of the labeled points) and large numbers to points close to the center.

See the results for \(d=1\) in Figure 37, \(d=2\) in Figure 38, and \(d=3\) in Figure 39. While for \(d\leq 2\) it appears to be smoothly interpolating between the labels, for \(d=3\) it seems that the method simply learns essentially \(-1\) on all points, thus not being very meaningful. Let us turn to \(\mathbb{R}^d\) for intuition:

Let us say that we want to find a function in \(\mathbb{R}^d\) that takes the value \(1\) at zero and \(-1\) at the unit sphere, that minimizes \(\int_{B_{0}(1)}\|\nabla f(x)\|^2 dx\). Let us consider the following function on \(B_0(1)\) (the ball centered at \(0\) with unit radius) \[f_\varepsilon(x) = \left\{ \begin{array}{cc} 1-2\frac{|x|}{\varepsilon} & \text{if} |x|\leq \varepsilon\\ -1 & \text{otherwise.} \end{array} \right.\] A quick calculation suggest that \[\int_{B_{0}(1)}\|\nabla f_\varepsilon(x)\|^2 dx = \int_{B_{0}(\varepsilon)}\frac{1}{\varepsilon^2}dx = \operatorname{vol}(B_0(\varepsilon)) \frac{1}{\varepsilon^2}dx \approx \varepsilon^{d-2},\] meaning that, if \(d>2\), the performance of this function is improving as \(\varepsilon\to 0\), explaining the results in Figure 39.

One way of thinking about what is going on is through the Sobolev Embedding Theorem. \(H^m\left( \mathbb{R}^d\right)\) is the space of function whose derivatives up to order \(m\) are square-integrable in \(\mathbb{R}^d\), Sobolev Embedding Theorem says that if \(m > \frac{d}{2}\) then, if \(f\in H^m\left( \mathbb{R}^d\right)\) then \(f\) must be continuous and bounded (i.e., \(f \in C_0(\mathbb{R}^d)\)), which would rule out the behavior observed in Figure 39. It also suggests that if we are able to control also second derivatives of \(f\) then this phenomenon should disappear (since \(2>\frac{3}{2}\)). While we will not describe it here in detail, there is, in fact, a way of doing this by minimizing not \(f^TLf\) but \(f^TL^2f\) instead, Figure 40 shows the outcome of the same experiment with the \(f^TLf\) replaced by \(f^TL^2f\) and confirms our intuition that the discontinuity issue should disappear (see, e.g., [@Nadler_SSL_sobolev_09] for more on this phenomenon).

7.7 A brief primer on reversible Markov chains↩︎

With the material already described in this chapter one is able to very easily describe the basics of the theory of reversible Markov chains, so we will take the opportunity to do this (while pointing the interested reader to [@levinperesmarkovchains2017] for a thorough treatment). We will focus here on Markov chains on finite state spaces³⁵.

A Markov chain on \(N\) states is a random walk on the state space, identified with \([N]\), whose transition probabilities depend only on the current state and not on the past (the so-called Markovian property). Let us describe the transition probabilities by \[\label{eq:MarkovChainTransitionMatrix} M_{ij} = \mathbb{P}\left( X(t+1)=j | X(t) = i\right).\tag{129}\]

Let \(\pi\in[0,1]^{N}\) (with \(\pi^\top\mathbf{1}=1\)) be a stationary measure of the Markov chain, meaning that \(\pi^\top M = \pi^\top\).³⁶ A Markov chain is said to be irreducible if there is no proper subset of states where the probability mass can be fully trapped, in other words: A Markov chain is irreducible if for all \(i,j\in[N]\), \(\exists_{t>0}\) such that \(\mathbb{P}\left( X(t)=j | X(0) = i\right)>0\). One can show that irreducible Markov chains have unique stationary measures (and the stationary measure needs to be strictly positive).³⁷

In the rest of this section we will focus on Markov chains with a unique stationary \(\pi\) that furthermore satisfies \(\pi>0\) (we will identify the Markov chain with the transition matrix \(M\) given in 129 ). Since \(\pi\) is a stationary measure of \(M\) we must have that, for each state \(i\), the probability inflow at the: \(\sum_{j\in[N]}\pi_jM_{ji}\) is the same as the probability outflow: \(\sum_{j\in[N]}\pi_iM_{ij}\). Note that this is indeed equivalent to stationarity as \(\sum_{j\in[N]}\pi_jM_{ji}=\left(\pi^\top M\right)_i\) and \(\sum_{j\in[N]}\pi_iM_{ij} = \pi_i\sum_{j\in[N]}M_{ij}=\pi_i\). A natural stronger condition is to request that the inflow and outflow are the same per edge, not just total per state: for \(\pi>0\) the stationary measure and for all \(i,j\in[N]\), \[\label{eq:reversibleMC} \pi_jM_{ji} = \pi_i M_{ij}.\tag{130}\] Condition 130 is known as the detailed balance condition. Markov chains satisfying detailed balance 130 are called reversible Markov chains, they form a very important class of Markov chains. A reversible Markov chain can be thought of as a Markov chain that, at stationarity (or equilibrium), has the same dynamics when viewed running forward or backwards in time.

Reversibility (condition 130 ) can be viewed as \(M\) being self-adjoint in the inner product \(\langle \cdot, \cdot \rangle_\pi\) (given by \(\langle x, y \rangle_\pi = \sum_{i=1}^N \pi_i x_i y_i\)). Namely \[\langle x, My \rangle_\pi = \sum_{i=1}^N\sum_{j=1}^N \pi_i x_i M_{ij}y_j = \sum_{i=1}^N\sum_{j=1}^N \pi_j x_i M_{ji}y_j = \langle M^\top x, y \rangle_\pi.\]

Given \(M\) a reversible Markov chain with a stationary measure \(\pi>0\), let us define \[\label{eq:WisdiagpiM} W = \mathop{\mathrm{diag}}(\pi)M,\tag{131}\] where \(\mathop{\mathrm{diag}}(\pi)\) is the diagonal matrix with diagonal entries given by \(\pi\). Reversibility is equivalent to \(W\) being a symmetric matrix, since \(W\geq 0\) it can be viewed as the weight matrix of a graph on \(N\) nodes. This means that any reversible Markov chain with a strictly positive stationary measure can be viewed as the diffusion map from some weighted graph \(G=([N],E,W)\) and thus the theory developed in this chapter (including connections with spectral clustering and Cheeger’s Inequality) can be directly used to understand reversible Markov chains.

A very important concept in the study of Markov Chains is that of mixing time: how many iterations of the Markov Chain are needed, when starting from an initial distribution \(\mu_0\), to reach a distribution that is at most \(\varepsilon\) away from \(\pi\) in Total Variation (TV) distance. The TV distance between two probability distributions \(p\) and \(q\) (defined on the same probability space) is given by \(\sup_E\|p(E)-q(E)\|\) where the supremum is taken over all measurable sets in the probability space. Since in our case the distributions are finite (they are distributions over [N]) we have \[\|\mu-\nu\|_{TV} =\frac{1}{2}\sum_{i\in [N]}\left| \mu_i - \nu_i \right| = \frac{1}{2}\|\mu-\nu\|_{1},\] where \(\mu_i\) and \(\nu_i\) are the probabilty of \(i\) under, respectively, \(\mu\) and \(\nu\), and \(\|\cdot\|_1\) denotes the \(\ell_1\)-norm in \(\mathbb{R}^N\).

A core idea in the study of Markov Chains is to bound mixing times via functional inequalities that involve the spectrum of the transition matrix. A central object in this context is the Dirichlet Form which measures how a function varies under the transitions of the Markov Chain: given a function \(f:[N]\to\mathbb{R}\) (here viewed as a vector \(f\in\mathbb{R}^N\)) the Dirichlet Form of \(M\) is given by \[\mathcal{E}_M(f) \mathrel{\mathop{:}}=\frac{1}{2} \sum_{i,j=1}^N \pi_i f_i (I-M)_{ij}f_j.\] The Markov Chain is said to satisfy a Poincaré inequality with constant \(C_M\) if, for all \(f\in\mathbb{R}^{N}\), \[\label{eq:poincare-for-MC} \mathop{\mathrm{Var}}_{\pi}(f) \leq C_M \mathcal{E}_M(f),\tag{132}\] where \(\mathop{\mathrm{Var}}_{\pi}(f)\) is the variance of \(f\) under the measure \(\pi\). In other words: \(\mathop{\mathrm{Var}}_{\pi}(f) = \sum_{i=1}^n \pi_i\left( f_i - \pi^\top f \right)^2\). Note that \(\pi^\top f=\mathbb{E}_\pi f\) corresponds to the mean of \(f\) under the distribution \(\pi\).

The best constant in 132 is known as the spectral gap of the Markov Chain: \[\label{eq:Gap-for-MC} \mathrm{gap}_M \mathrel{\mathop{:}}=\inf_{\substack{f\in\mathbb{R}^N \\ \mathop{\mathrm{Var}}_{\pi}(f)\neq 0 }}\frac{\mathcal{E}_M(f)}{\mathop{\mathrm{Var}}_{\pi}(f)}.\tag{133}\]

Note that we have \[\label{eq:DirichletAndLaplacian} \mathcal{E}_M(f) = \frac{1}{2} f^\top \mathop{\mathrm{diag}}(\pi) (I-M) f = \frac{1}{2} f^\top (D_W-W) f = \frac{1}{2} f^\top L_W f,\tag{134}\] where \(W\) is the weight matrix of a graph, as in 131 , and \(D_W=\mathop{\mathrm{diag}}(\pi)\) and \(L_W=D_W-W\) are the associated degree and graph Laplacian matrices, respectively. Thus, \(\mathop{\mathrm{Var}}_{\pi}(f)=0\Rightarrow \mathcal{E}_M(f)=0\) and \(\mathrm{gap}_M\) in 133 indeed corresponds to the best Poincaré constant in 132 . Similarly, noting that \(\pi^\top f = \mathbf{1}^\top D_W f\), \[\begin{align} \mathop{\mathrm{Var}}_\pi(f) & = & \sum_{i=1}^n \pi_i\left( f_i - \pi^\top f \right)^2 = \mathbb{E}_\pi f^2 - \left( \mathbb{E}_\pi f \right)^2 \nonumber\\ & = & f^\top D_W f - \left(\mathbf{1}^\top D_W f\right)^2. \label{eq:VarAndDW} \end{align}\tag{135}\]

Recall that, by construction \(\operatorname{Tr}D_W = \sum_{i=1}^N \pi_i = 1\) and that we assume in this section that \(\pi_i>0\) for all \(i\in[N]\).

Recall from Chapter 6 that the normalized graph Laplacian is given by \(\mathcal{L}_W = I - D_W^{-\frac{1}{2}} W D_W^{-\frac{1}{2}}\) and its second smallest eigenvalue is given (as we saw in 97 and 99 ) by \[\lambda_2(\mathcal{L}_W) = \inf_{\substack{y\in\mathbb{R}^N\setminus\{0\} \\ \mathbf{1}^\top D_W y = 0}} \frac{y^\top L_W y}{y^\top D_W y}.\]

Note that adding a multiple of \(\mathbf{1}\) to \(f\in\mathbb{R}^N\) does not affect either \(\mathcal{E}_M(f)\) or \(\mathop{\mathrm{Var}}_\pi(f)\) (see 134 ). This means we can write \[\mathrm{gap}_M = \inf_{\substack{f\in\mathbb{R}^N\setminus \{0\}\\ \mathbf{1}^\top D_W f=0 }}\frac{\mathcal{E}_M(f)}{\mathop{\mathrm{Var}}_{\pi}(f)} = \inf_{\substack{f\in\mathbb{R}^N\setminus \{0\}\\ \mathbf{1}^\top D_W f=0 }}\frac{f^\top L_W f}{f^\top D_W f} = \lambda_2(\mathcal{L}_W) = 1-\lambda_2(M),\] where \(\lambda_2(M)\) is the second largest eigenvalue of \(M\). The second inequality used 134 and 135 . The fourth inequality follows from the fact that \(M=\mathop{\mathrm{diag}}(\pi)^{-1}W\) (recall 131 ) is similar to \(\mathop{\mathrm{diag}}(\pi)^{-\frac{1}{2}}W\mathop{\mathrm{diag}}(\pi)^{-\frac{1}{2}}\) and \(\mathcal{L}_W = I - \mathop{\mathrm{diag}}(\pi)^{-\frac{1}{2}}W\mathop{\mathrm{diag}}(\pi)^{-\frac{1}{2}}\).

Similarly to the situation in Remark 7.2 it is simpler to focus on lazy walks. We say a Markov Chain \(M\) is lazy if \(M_{ii}\geq \frac{1}{2}\) for all \(i\in[N]\) (this ensures that all eigenvalues of \(M\) are non-negative). Given a transition matrix \(M\) we can consider a lazy version of it by taking \(M_{\text{lazy}}=\frac{1}{2}M+\frac{1}{2}I\). This does not affect reversibility and reduces the spectral gap by half.

We now show a central result in this area: a control on the spectral gap (or, equivalently on \(\lambda_2(\mathcal{L}_W)\)) provides a bound on the mixing time of a Markov Chain.

Theorem 7.8 shows that mixing times are governed by a Fiedler eigenvalue \(\lambda(\mathcal{L}_W)\) of the underlying weighted graph. Cheeger’s inequality. Theorem 6.5, further shows that these are directly related to bottlenecks (small cuts) in the graph. These results together can be interpreted as saying that if there are no bottlenecks in the graph, the associated Markov Chain necessarily mixes fast.

It is generally possible to improve the dependency on \(\pi_{\min}\) on the control of mixing times. A (now) standard way is through so-called Log-Sobolev inequalities. The idea is to develop analogues to Poincaré inequalities that control the entropy of a test function, rather than the variance. This tends to allow one to replace the \(\log(\pi_{\min}^{-1})\) term by a \(\log\log(\pi_{\min}^{-1})\) term. See, e.g., [@Diaconisetal-LogSobolev96].

As we will see below, many Markov Chains of interest have exponential state spaces (such as \(\{\pm1\}^n\), meaning that \(N=2^n\)) and the goal is to have polynomial in \(n\) mixing times.

7.7.1 A word about Markov Chain Monte Carlo methods↩︎

An important use of Markov Chains is to sample from high dimensional distributions that are easy to access up to normalization. This strategy usually goes by the name of Markov Chain Monte Carlo [@hastings1970monte; @robert1999monte]. These are central objects in both Statistical Physics [@krauth2006statistical] and high dimensional Bayesian Statistics [@gelman1995bayesian].

Let us focus on a (rich) class of examples for the sake of exposition: The goal is to sample from a probability distribution over the hypercube \(\{\pm1\}^n\) (meaning that \(N=2^n\)), which \(\pm1\) entry is usually called a spin. Let \(J\in\mathbb{R}^{n\times n}\) be a symmetric matrix. Given \(\beta\geq 0\) we define the following probability distribution over the hypercube, known as a Gibbs measure³⁸: \[\label{eq:gibbsmeasurexTJx} \pi_x = \frac{1}{Z_n}\exp{\left(\beta x^TJx\right)},\tag{140}\] for \(x\in\{\pm1\}^n\), where \(Z_n = \sum_{x\in\{\pm1\}^n}\exp{\left(\beta x^TJx\right)}\) is known as a the normalization constant, or partition function.

Distributions of the form of 140 are central in both Statistical Physics and Bayesian Statistics. In Statistical Physics, by taking different choices of \(J\), it includes very important models such as Curie-Weiss, Ising, Sherrington-Kirkpatrick, and Edward-Anderson models. It is also central in Statistics. In fact the posterior distribution of the node memberships given the graph edges in the Stochastic Block Model (in two communities) of Chapter 13 is naturally written in the form 140 , see [@Abbe_SBM_survey].

The goal will be to construct a Markov Chain whose unique stationary measure is the target distribution. A key property of these measures is that while they are hard to compute (because computing \(Z_n\) requires summing over exponentially many summands) they are easy to access up to normalization, meaning that it is easy to compute relative probabilities, such as \(\pi_x/\pi_y\) for two states \(x\) and \(y\).

Motivated by the example above let us show the construction of a reversible Markov Chain to sample from a (positive) probability distribution \(\pi\) in \([N]\). We start by building a graph on \([N]\), and for the sake of exposition, let us say this is a simple unweighted \(d\)-regular graph with adjacency \(W\). Each step of the Markov Chain proceeds as follows: From state \(x\in[N]\), the chain picks a neighbor uniformly at random, let us call that neighbor \(y\in[N]\); then the Markov Chain moves to \(y\in[N]\) with probability \(\frac{\pi_y}{\pi_x+\pi_y}\) and stays in \(x\) with probability \(\frac{\pi_x}{\pi_x+\pi_y}\) (note that while computing \(\pi_x\) is difficult for, e.g. 140 , computing these update probabilities is easy). The reader can check that this Markov Chain is reversible and that \(\pi\) is indeed a stationary measure. This essentially corresponds to the celebrated Metropolis algorithm (see, e.g. [@Newman-Barkema-MCMC99]).³⁹

Returning to the hypercube and the target distribution 140 , a classical Markov Chain is obtained by taking the contribution above for the graph whose nodes are elements of \(\{\pm1\}^n\) and two nodes are connected if they differ in only one entry (often called a spin). In this case \(N=2^n\) and \(d=n\). This gives rise to the celebrated Glauber Dynamics (see [@levinperesmarkovchains2017]).⁴⁰ In many situations this process allows one to sample distributions of the form of 140 in time that is polynomial in \(n\) (and so logarithmic in \(N\)). There has been exciting recent progress on understanding mixing times of these (and other related) chains [@AnariLiuOveisGharanFOCS2020; @ChenEldan2025LocalizationSchemes].

Exercises↩︎

Prove that $F$ diagonalizes $C$, meaning that $FCF^*$ is a diagonal
matrix.

8 Linear Dimension Reduction via Random Projections↩︎

In Chapters 4 and 7 we saw both linear and non-linear methods for dimension reduction. In this chapter we will see one of the most fascinating consequences of the phenomenon of concentration of measure in high dimensions, one of the blessings of high dimensions described in Chapter 3. When given a data set in high dimensions, we will see that the projection to a lower dimensional space, taken at random, can preserve (under conditions made precise in the next section) certain geometric features of the data set. The remarkable aspect here is that this “lower” dimension can be strikingly lower. This gives rise to significant computational savings in many data processing tasks by including a random projection as a pre-processing step. Similar ideas are also the driving force behind many algorithms in randomized linear algebra, as we will demonstrate in this chapter via the randomized SVD. There is however another less obvious implication of this phenomenon with important practical implications: since the projection is agnostic of the data, it can be leveraged even when the data set is not explicit, such as the set of all natural images or the set of all “possible” brain scans; this is at the heart of Compressive Sensing, which we will explore in Chapter 16

8.1 The Johnson-Lindenstrauss Lemma↩︎

Suppose one has \(n\) points, \(X = \{x_1,\dots,x_n\}\), in \(\mathbb{R}^p\) (with \(p\) large). If \(p>n\), the points actually lie in a subspace of dimension \(n\), so the projection \(f:\mathbb{R}^p\to \mathbb{R}^n\) of the points to that subspace acts without distorting the geometry of \(X\). In particular, for every \(x_i\) and \(x_j\), \(\|f(x_i)-f(x_j)\|^2 = \|x_i-x_j\|^2\), meaning that \(f\) is an isometry in \(X\). Suppose instead we allow a bit of distortion, and look for a map \(f:\mathbb{R}^p\to \mathbb{R}^d\) that is an \(\varepsilon-\)isometry⁴¹, meaning that \[\label{def95epsilonisometry} (1-\varepsilon) \|x_i-x_j\|^2 \leq \|f(x_i)-f(x_j)\|^2 \leq (1+\varepsilon) \|x_i-x_j\|^2.\tag{143}\] Can we do better than \(d=n\)?

In 1984, Johnson and Lindenstrauss [@Johnson_Lindenstrauss_84] showed a remarkable lemma that answers this question affirmatively.

Note that by considering \(d\) slightly larger one can get a good projection on the first random attempt with high confidence. In fact, it is trivial to adapt the proof above to obtain the following lemma:

Proposition 8.2 is quite remarkable. Consider the situation where we are given a high-dimensional data set in a streaming fashion – meaning that we get each data point at a time, consecutively. To run a dimension-reduction technique like PCA or Diffusion maps we would need to wait until we received the last data point and then compute the dimension reduction map (both PCA and Diffusion Maps are, in some sense, data adaptive). Using Lemma 8.2 one can just choose a projection at random in the beginning of the process (all one needs to know is an estimate of the logarithm of the size of the data set) and just map each point using this projection matrix which can be done online – we do not need to see the next point to compute the projection of the current data point. Proposition 8.2 ensures that this (seemingly naïve) procedure will, with high probably, not distort the data by more than \(\varepsilon\).

One might wonder if such low-dimensional embeddings as provided by the Johson-Lindenstrauss Lemma also extend to other norms besides the Euclidean norm. For the \(\ell_1\)-norm there exist lower bounds which prevent such a dramatic dimension reduction (see [@lee2004embedding]), and for the \(\ell_{\infty}\)-norm one can easily construct examples that demonstrate the impossibility of dimension reduction. Hence, the Johnson-Lindenstrauss Lemma seems to be a subtle result about the Euclidean norm.

8.1.1 The Fast Johnson-Lindenstrauss transform↩︎

Let us continue thinking about our example of high-dimensional streaming data. After we draw the random projection matrix⁴³, say \(M\), for each data point \(x\), we still have to compute \(Mx\) which has a computational cost of \(\mathcal{O}(\varepsilon^{-2}\log(n)p)\) since \(M\) has \(\mathcal{O}(\varepsilon^{-2}\log(n)p)\) entries (since \(M\) is a random matrix, generically it will be a dense matrix). In some applications this might be too expensive, raising the natural question of whether one can do better. Moreover, storing a large-scale dense matrix \(M\) is not very desirable either. There is no hope of significantly reducing the number of rows in general, as it is known that the Johnson-Lindenstrauss Lemma is orderwise optimal [@Alon_JL_lowerbound_2003; @Larsen_Nelson_JL].

We might hope to replace the dense random matrix \(M\) by a sparse matrix \(S\) to speed up the matrix-vector multiplication and to reduce the storage requirements. This method was proposed and analyzed in [@Ailon_Chazelle_fast_JL]. Here we discuss a slightly simplified version, see also [@cohen2015optimal].

We let \(S\) be a very sparse \(d \times p\) matrix, where each row of \(S\) has just one single non-zero entry of value \(\sqrt{p/d}\) at a location drawn uniformly at random. Then, for any vector \(x \in \mathbb{R}^p\) \[\underset{i}{\mathbb{E}} [(Sx)_i^2]=\sum\limits_{j=1}^p \mathbb{P}(S_{ij}\neq 0) \cdot \frac{p}{d} \cdot x^2_j = \frac{1}{d} \|x\|^2,\] hence \(\mathbb{E}[\|Sx\|^2] = \mathbb{E}[\sum\limits^d_{i = 1} (Sx_i)^2]=\|x\|^2\). This result is satisfactory with respect to expectation (even for \(d=1\)), but not with respect to the variance of \(\|Sx\|^2\). For instance, if \(x\) has only one non-zero entry we need \(d \sim \mathcal{O}(p)\) to ensure that \(\|Sx\|^2 \neq 0\) with non-negligible probability. More generally, if one coordinate of \(x\) is much larger (in absolute value) than all its other coordinates, then we will need a rather large value for \(d\) to guarantee that \(\|Sx\| \approx \|x\|\).

A natural way to quantify the “peakiness” of a vector \(x\) is via the peak-to-average ratio⁴⁴ measured by the quantity \(\|x\|_\infty/\|x\|\). It is easy to see that we have (assuming \(x\) is not the zero-vector) \[\frac{1}{\sqrt{p}} \le \frac{\|x\|_\infty}{\|x\|} \le 1.\] The upper bounds is achieved by vectors with only one non-zero entry, while the lower bound is met by constant-modulus vectors. Thus, if \[\label{approxpar} \frac{\|x\|_\infty}{\|x\|} \approx \frac{1}{\sqrt{p}},\tag{145}\] we can hope that sparse subsampling of \(x\) will still preserve its Euclidean norm.

Thus, this suggests to include a preprocessing step by applying a rotation so that sparse vectors become non-sparse in the new basis, thereby reducing their \(\infty\)-norm (while their 2-norm remains invariant under rotation). Two natural choices for such a rotation are the Discrete Fourier transform (which maps unit-vectors into constant modulus vectors) and its \(\mathbb{Z}_2\)-cousin, the Walsh-Hadamard matrix⁴⁵. But since the Fast Johnson-Lindenstrauss Transform (FJLT) has to work for all vectors, we need to avoid that this rotation maps dense vectors into sparse vectors. We can address this issue by “randomizing” the rotation, thereby ensuring with overwhelming probability that dense vectors are not mapped into sparse vectors. This can be accomplished in a numerically efficient manner (thus maintaining our overall goal of numerical efficiency) by first randomizing the signs of \(x\) before applying the rotation. Putting these steps together we arrive at the following map.

We can carry out the matrix-vector multiplication by the DFT matrix via the Fast Fourier Transform (FFT) in \(\mathcal{O}(p \log p)\) operations; a similar algorithm exists for the Walsh-Hadamard matrix. The FJLT allows for a dimension reduction that is competitive with the Johnson-Lindenstrauss Lemma as manifested by the following theorem.

For convenience we will prove Theorem 8.3 for the case when \(F\) is the Walsh-Hadamard matrix so that the random variables are real-valued and the exposition is lighter. But it is straighforward to adapt the argument to the case when \(F\) is the DFT matrix. Keeping this minor simplification in mind, the proof of Theorem 8.3 follows from the two lemmas below. We first show that with high probability the random rotation \(FD\) produces vectors with a sufficiently low peak-to-average ratio.

Combining Lemmas 8.1 and 8.2, union bounding over the two events, and recalling that if \(\|x\|=1\) then \(\|y\|=1\) since both \(D\) and \(F\) are unitary, establishes Theorem 8.3.

Besides the potential speedup, another advantage of the FJLT is that it requires significantly less memory compared to storing an unstructured random projection matrix as is the case for the standard Johnson-Lindenstrauss approach.

There are many applications where one wants to perform dimension reduction while preserving the geometry of certain sets with infinite points, such as subspaces. Via an \(\varepsilon\)-net argument it is possible to show that random matrices that are good as Johnson-Lindenstrauss maps for \(N\) points, also approximately preserve the geometry of a (\(\log N\))-dimensional subspace. These go by the name of Oblivious Subspace Embeddings (see, e.g., Section 5.2.2. in [@KireevaTropp] where the fast JL construction we introduced above is shown to be an Oblivious Subspace Embedding). A similar idea will be explored below (Section 8.2) where a random projection will allow us to still approximate the subspace spanned by the singular vectors corresponding to the largest singular values, resulting in a randomized SVD algorithm that is workhorse of large scale linear algebra. In Chapter 14 we will develop analogues of the Johnson-Lindenstrauss Lemma for more general sets, through the celebrated Gordon’s espace through a mesh theorem.

8.2 Randomized SVD↩︎

Randomized linear algebra has gained significant relevance in recent years due to its ability to efficiently handle large-scale data problems that traditional methods struggle with. The Johnson-Lindenstrauss dimension reduction discussed in the previous section is one prominent example. The techniques encountered in the Johnson-Lindenstrauss transform naturally lend themselves as key ingredient of the randomized SVD, a method to compute an approximation of the SVD of a matrix. It is especially useful when the matrix is large and dense and can be well approximated by a low-rank matrix.

In this section we will give a detailed description of the randomized SVD. It is based on using random projections to reduce the dimensionality of the matrix before performing a more traditional SVD on the reduced matrix.

Recall from Chapter 4 that the SVD of a matrix \(A\in\mathbb{R}^{m\times n}\) is given by \[\label{eq:SVDrand} A = U \Sigma V^{\ast},\tag{147}\] where \(U \in \mathbb{R}^{m \times m}\) is a unitary matrix containing the left singular vectors, \(\Sigma \in \mathbb{R}^{m \times n}\) is a diagonal matrix with singular values (in decreasing order) on the diagonal, and \(V \in \mathbb{R}^{n \times n}\) is a unitary matrix containing the right singular vectors.

In practice \(A\) is often large, and we are interested in only a small number \(s\) of the largest singular values and corresponding singular vectors, giving a rank-\(k\) approximation of \(A\): \[A \approx U_k \Sigma_k V_k^{\ast}\] where \(U_k \in \mathbb{R}^{m \times k}\), \(\Sigma_k \in \mathbb{R}^{k \times k}\), and \(V_k \in \mathbb{R}^{n \times k}\).

We define the projection \(P_k\) onto the the span of the leading \(k\) left singular vectors of \(A\) via \[P_k = \sum_{i=1}^k u_i u_i^{\ast}.\] We know from Theorem 4.5 that the best rank-\(k\) approximation of \(A\) with respect to the Frobenius norm is given by \[\min_{\operatorname{rank}(B) \leq k} \|A - B \|_F = \|A - P_k A \|_F = \sqrt{\sum_{i > k} \sigma_i^2}.\] Thus, we can compute the best approximation by projecting the target matrix onto the \(k\)-dimensional subspace that captures “most of the action” of the matrix \(A\).

In many applications it suffices to compute the best approximation with moderate accuracy. Hence, following [@HMT11], we let \(s=k+p\) for small \(p\), and seek a rank-\(s\) approximation \(\hat{A}_s\) that is comparable to the best rank-\(k\) approximation, i.e., \[\|A - \hat{A}_s\|_F^2 \le (1+\varepsilon) \|A - P_k A \|_F^2 = (1+\varepsilon)\sum_{i > k} \sigma_i^2\] for some chosen tolerance \(\varepsilon> 0\).

The aim of the Randomized SVD (RSVD) is essentially to generate an efficient approximation for the best rank-\(k\) approximation of the matrix \(P_k A\). We proceed in two steps. In the first step we compute a low-dimensional subspace that captures the range of \(A\). In the second step we project \(A\) onto the subspace and compute the (standard SVD) of the resulting smaller matrix. Regarding the first step, recall that \(P_k\) is the orthogonal projection onto the \(k\) leading left singular vectors of \(A\). But since we do not know the subspace spanned by \(P_k\), a sensible choice, as suggested by the Johnson-Lindenstrauss Lemma, is to draw some probing random vectors from a rotationally invariant distribution to estimate the directions of the vectors \(u_i\) that make up \(P_k\). We collect these probing random vectors in a matrix \(\Psi\).

Hence, let \(\Psi \in \mathbb{R}^{n \times (k + p)}\) be a random matrix, where \(k\) is the desired number of singular values/vectors, and \(p\) is an oversampling parameter (typically \(p\) is small, like 5 or 10). A promising choice for \(\Psi\) is to select each entry of \(\Psi\) as a sample from a standard normal distribution \(\mathcal{N}(0, 1)\). We then form the matrix \(Y \in \mathbb{R}^{m \times (k + p)}\) as \[Y = A \Psi.\] This matrix \(Y\) captures the action of \(A\) on a random subspace of dimension \(k + p\). We will see (in Theorem 8.4) that, with high probability, the subspace spanned by the columns of \(Y\) contains a good approximation of the subspace spanned by the top \(k\) left singular vectors of \(A\).⁴⁶ Next, we find an orthonormal basis for the range of \(Y\). Using the standard QR decomposition, we compute \[Y = Q R,\] where \(Q \in \mathbb{R}^{m \times (k + p)}\) has orthonormal columns, and \(R \in \mathbb{R}^{(k + p) \times (k + p)}\) is an upper triangular matrix. Now, we project \(A\) onto the lower-dimensional subspace spanned by the columns of \(Q\) by computing \[B = Q^{\ast} A.\] Here, \(B \in \mathbb{R}^{(k + p) \times n}\) is a much smaller matrix that approximates \(A\) (in the sense that \(QB\) approximates \(A\)).

Next, we compute the exact SVD of the smaller matrix \(B\) \[B = U_B \Sigma_B V_B^{\ast}.\] where \(U_B \in \mathbb{R}^{(k + p) \times (k + p)}\), \(\Sigma_B \in \mathbb{R}^{(k + p) \times n}\), and \(V_B \in \mathbb{R}^{n \times (k + p)}\). Finally, the approximate SVD of \(A\) is given by \[A \approx \tilde{ U} \tilde{ \Sigma} \tilde{V}^{\ast}\] where \(\tilde{U} = Q U_B \in \mathbb{R}^{m \times (k + p)}\), \(\tilde{\Sigma} = \Sigma_B \in \mathbb{R}^{ (k + p) \times (k + p)}\), and \(\tilde{V} = V_B \in \mathbb{R}^{n \times (k + p)}\).

We summarize these steps in Algorithm 41.

The random matrix \(\Psi\) used in the random projection step ensures that with high probability, the subspace spanned by the columns of \(Y = A \Psi\) contains a good approximation to the dominant singular subspace of \(A\).

The error on the right hand side consists of the best rank-\(k\) approximation error and an additional factor due to using a randomized projection instead of the optimal projection \(P_k\).

For the proof of Theorem 8.4 we need some preparation. In what follows it will be convenient to rewrite the SVD of \(A\), cf. 147 , in block matrix form: \[A = \begin{bmatrix} U_k & U_\perp \end{bmatrix} \begin{bmatrix} \Sigma_k & 0 \\ 0 & \Sigma_\perp \end{bmatrix} \begin{bmatrix} V_k^{\ast} \\ V_\perp^{\ast} \end{bmatrix},\] where \(U_\perp \in \mathbb{R}^{m \times m-k}\) denotes the orthogonal complement of \(U_k \in \mathbb{R}^{m \times k}\) with respect to \(U \in \mathbb{R}^{m \times m}\). Furthermore, we set \[\Psi_k : = V_k^\ast \Psi \qquad \text{and} \qquad \Psi_\perp : = V_\perp^\ast \Psi .\] Here, \(\Psi_k\) captures the alignment of \(\Psi\) with the leading right singular vectors of \(A\), while \(\Psi_\perp\) captures the alignment with the trailing right singular vectors of \(A\). The matrix \(P_Y\) denotes the orthogonal projection onto the range of \(Y\). Since this projection is unique, we have that \(P_Y = Q Q^{\ast}\).

The following deterministic error bound (the bound does not require \(\Psi\) to be Gaussian) will be instrumental:

Lemma 8.3 holds for any rank-\(k\) matrix \(\Psi\). As pointed out in [@KireevaTropp], the first term on the right-hand side of 149 captures the trailing singular values of the matrix, which is independent of the specific choice of \(\Psi\), since this loss occurs for every rank-\(k\) approximation. The second term is the error caused by the dimension reducing matrix \(\Psi\). We want the component \(\Psi_k\) in the leading directions to be well conditioned (ideally, an identity matrix). We want the component \(\Psi_\perp\) in the trailing directions to be as small as possible.

The bound 149 suggests that in order to establish Theorem 8.4 we need to get a handle on \(\Psi_k \Psi_\perp^{\dagger}\) when \(\Psi_{ij} \overset{\text{i.i.d.}}{\sim} {\mathcal{N}}(0,1).\) To that end, the following lemma will be helpful.

With these lemmata in place, we are now ready to prove Theorem 8.4.

8.2.1 Computational complexity of the Randomized SVD↩︎

What is the computational complexity of the RSVD? The dominant costs for the RSVD are:

1. \({\mathcal{O}}(m n (k + p))\) for the matrix-matrix multiplication \(A \Psi\), assuming that \(\Psi\) is an unstructured, non-sparse matrix (like a Gaussian matrix).

2. \({\mathcal{O}}(m (k + p)^2)\) for the QR decomposition of \(Y\).

3. \({\mathcal{O}}(mn (k + p))\) for computing \(Q^{\ast} A\).

4. \({\mathcal{O}}((k + p)^2 n)\) for computing the SVD of the small matrix \(B\).

Thus, the overall complexity of the randomized SVD is \({\mathcal{O}}(m n (k + p))\), which is significantly smaller than the \({\mathcal{O}}(m n^2)\) cost of the classical SVD when \(k\) is much smaller than \(n\).

We can further reduce the cost of the matrix-matrix product \(A \Psi\) by introducing structure or sparsity into \(\Psi\). We already encountered one such choice in Definition 8.1 in the form of the Fast Johnson-Lindenstrauss transform. The improved computational efficiency usually comes at the cost of additional log-factors in terms of required samples to guarantee the same accuracy.

When the singular values of \(A\) decay slowly, Algorithm 41 is not very accurate in the regime of interest, namely when \(k \ll n\). The reason is that to achieve a small approximation error the matrix \(Q\) must align well with the first \(k + p\) left singular vectors of \(A\). As a consequence, the product \(A\Psi\) must reveal the singular vectors, which in turn would require the singular vectors associated with small singular values not to interfere with the calculation.

How can we remedy this issue? We can take inspiration from the power iteration. By raising our matrix to a power \(q\) via \((AA^{\ast} )^q A\) prior to applying it to \(\Psi\) we can enhance the decay of the spectrum of our matrix without effecting the direction of the eigenvectors. Hence, we replace in Algorithm 41 the step \(Y = A \Psi\) by \(Y = (AA^{\ast} )^q A \Psi\). A small power of \(q\) (such as \(q=1\) or 2) is often sufficient. In this modified RSVD it is also advisable to increase the oversampling factor \(p\). We refer to [@HMT11] for an approximation error bound for this modified RSVD as well as to [@tropp2023randomized] for a more detailed discussion.

The key ideas behind the RSVD also form the main ingredients for some other algorithms, such as randomized least squares and randomized Choleksy factorization, see [@rokhlin2008fast; @HMT11; @KireevaTropp].

8.3 Matrix sketching↩︎

The methods we encountered in this chapter are closely intertwined with the concept of matrix sketching [@liberty2013simple; @mahoney2011randomized; @rudelson2007sampling; @HMT11], a family of randomized techniques for constructing a low-dimensional representation of a large matrix while approximately preserving its essential linear-algebraic structure.

Given a data matrix \(A \in \mathbb{R}^{n \times d}\), the goal is to replace \(A\) by a much smaller matrix—called a sketch—that enables accurate approximation of quantities such as norms, inner products, or dominant singular subspaces.

Two complementary formulations are common. In the projection-based viewpoint, one forms a sketch \(SA\), where \(S \in \mathbb{R}^{m \times n}\) is a random embedding (e.g., Gaussian or subsampled randomized Hadamard), generalizing the Johnson–Lindenstrauss lemma from vectors to entire subspaces [@woodruff2014sketching; @tropp2017practical].

In the sampling-based viewpoint, one constructs a sketch by randomly selecting and rescaling columns (or rows) of \(A\), often using leverage-score sampling to retain the most informative directions [@frieze2004fast; @mahoney2009cur; @drineas2006subspace]. Both approaches can be understood as approximately preserving the geometry of the column space of \(A\).

Matrix sketching does not only underlie the randomized SVD, it provides a principled way to trade controlled approximation error for potentially substantial reductions in computational and memory costs. As such, it enables distributed linear algebra [@woodruff2014sketching] and it is the fundamental mathematical too that makes streaming possible for high-dimensional data [@harvey2008sketching; @mahoney2011randomized]. Here, streaming refers to a computational model in which data are processed sequentially, in one (or very few) passes, using limited memory, without storing the entire dataset [@muthukrishnan2005data; @woodruff2014sketching]. Streaming is essential when the dataset is too large to store in memory, data arrive continuously (logs, sensor data, social networks), or computation must be done online.

Let us explore some exemplary results in sketching beyond the Johnson-Lindenstrauss transform and the randomized SVD.

We will occasionally abuse notation and say that \(S\) is an \(\ell_2\)-subspace embedding for \(A\) itself, although it should be clear from the definition that this property is independent from the choice of a particular basis for the representation of the column space of \(A\).

It is easy to see from our previously derived concentration bounds that a Gaussian random matrix and a Subsampled Randomized Hadamard (or Fourier) Transform as introduced in Definition 8.1 satisfy the conditions of Theorem 8.5.

Another interesting and useful construction of a sketching matrix is the so-called CountSketch [@charikar2002finding; @thorup2012tabulation]. CountSketch is extremely sparse, making it an appealing choice for processing massive datasets or high-speed data streams. It is often used to estimate frequencies of elements in data streams.

A CountSketch matrix \(S \in \mathbb{R}^{m \times n}\) is defined by two independent hash functions:

1. \(h: \{1, \dots, n\} \to \{1, \dots, m\}\): This bucket hash function maps each column index to one of the \(m\) rows.

2. \(g: \{1, \dots, n\} \to \{+1, -1\}\): This sign hash function assigns a random sign to each column.

The entries of \(S\) are constructed as follows: \[S_{i,j} = \begin{cases} g(j) & \text{if } h(j) = i, \\ 0, & \text{otherwise}. \end{cases}\]

The CountSketch construction enjoys the following useful properties: (i) Sparsity: Each column of \(S\) has exactly one non-zero entry. (ii) Efficiency: To compute \(Sx\), we do not perform a standard matrix-vector product. Instead, for each entry \(x_j\), we simply update one row of the sketch: \(sketch[h(j)] \leftarrow sketch[h(j)] + g(j)x_j\). (iii) Time Complexity: Computing \(SA\) for \(A \in \mathbb{R}^{n \times d}\) takes only \(\mathcal{O}(\text{nnz}(A))\) (the number of non-zero entries in \(A\)) operations.

A typical realization of \(S\) might look like this: \[S = \begin{bmatrix} 0 & 0 & 0 & 0 & -1 & 0 \\ +1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & +1 & 0 & 0 \\ 0 & -1 & 0 & 0 & 0 & +1 \\ 0 & 0 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}.\] In this specific example, the two underlying hash functions performed the following actions: (i) The bucket hash (\(h\)): \(h(1)=2, h(2)=4, h(3)=5, h(4)=3, h(5)=1, h(6)=4\). Note that row 4 received two assignments (from columns 2 and 6, respectively), while row 6 is empty. This “collision” is expected in hashing. (ii) The sign hash (\(g\)): \(g(1)=+1, g(2)=-1, g(3)=-1, g(4)=+1, g(5)=-1, g(6)=+1\).

Before proving subspace embeddings, we establish that CountSketch preserves the inner product.

The proof is left to the reader in Exercise 8.7.

The attentive reader will notice a couple of places where this argument is rather suboptimal: the tail bound is obtained with Chebyshev’s inequality, and the quantity controlled is the Frobenius norm rather than the spectral norm (indeed, many \(d\times d\) random matrices tend to have Frobeniues norm that is about \(\sim\sqrt{d}\) times larger than the spectral norm). Using random matrix tools similar to the ones we will develop in Chapter 15, it is possible to significantly improve Theorem 8.6 (see e.g. [@nelson2013osnap]). In any case, while CountSketch is efficient, it requires a larger sketch size than Gaussian matrices to achieve the same \(\ell_2\)-subspace embedding guarantee. But CountSketch is the preferred choice when \(n\) and \(d\) are so large that even a single dense matrix multiplication is prohibitive.

The following theorem illustrates the use of sketching to speed up the solution of large-scale least squares problems, as they arise for example in high-dimensional linear regression.

8.3.1 Fast approximate matrix multiplication via sketching?↩︎

Computing the product of two matrices is one of the most common and basic procedures in numerical linear algebra. Let \(A,\,B\) be two arbitrary \(n\times n\) matrices. The cost for the classical way (the standard inner product method) to compute the product \(AB\) is \({\mathcal{O}}(n^3)\).

In 1969 Strassen took the linear algebra community by surprise when he showed that the standard matrix multiplication algorithm is not optimal from a computational viewpoint by presenting an algorithm that uses only \({\mathcal{O}}(n^{2.808})\) operations [@Strassen1969]. Strassen’s algorithm proceeds as follows. Let us decompose the product of \(A\) and \(B\) into blocks in the following way \[\begin{bmatrix} A_{11}&A_{12} \\ A_{21}&A_{22} \end{bmatrix} \begin{bmatrix} B_{11}&B_{12} \\ B_{21}&B_{22} \end{bmatrix} = \begin{bmatrix} C_{11}&C_{12} \\ C_{21}&C_{22} \end{bmatrix},\] where each block is of size \(\frac{n}{2}\times\frac{n}{2}\). Following the standard way of multiplication, this requires eight \(\left(\frac{n}{2}\times\frac{n}{2}\right)\) matrix multiplications and four additions of \(\frac{n}{2}\times\frac{n}{2}\) matrices. Strassen devised a formula that computes this block matrix multiplication using only seven \(\frac{n}{2}\times\frac{n}{2}\) matrix multiplication and eighteen \(\frac{n}{2}\times\frac{n}{2}\) additions. The cost for multiplying two \(\frac{n}{2}\times\frac{n}{2}\) matrices is \({\mathcal{O}}\left(\frac{n}{2}\right)^3\) and the cost for adding two \(\frac{n}{2}\times\frac{n}{2}\) matrices is \({\mathcal{O}}\left(\frac{n}{2}\right)^2\). Since matrix addition is cheaper than matrix multiplication, there is a net gain in computational cost. Proceeding recursively, with blocks of size \(\frac{n}{2}\), \(\frac{n}{4}\), …, (up to blocks of size 16) we end up with Strassen’s algorithm. The total cost now is \(\mathcal{O}(n^{\log_2 7})\sim \mathcal{O}(n^{2.8074})\).

Coppersmith and Winograd improved Strassen’s algorithm further, currently sitting at around \({\mathcal{O}}(n^{2.37})\) operations at the cost of an increased constant.

An important open question is whether we can we construct a fast, practical algorithm⁴⁷ that can multiply two arbitrary \(n\times n\) matrices in \({\mathcal{O}}(n^{2+\varepsilon})\) (or, say, \({\mathcal{O}}(n^{2}\log n)\)) operations.

Yet, in certain cases one may be content with computing \(AB\) approximately instead of exactly. By sacrificing precision, this should allow for a reduced complexity algorithm. Randomized algorithms are a natural candidate to realize this goal.

Let \(A \in \mathbb{R}^{m \times n}\) and \(B \in \mathbb{R}^{n \times p}\). An intuitive approach, popularized by Drineas, Kannan, and Mahoney [@drineas2006fast], treats the product \(AB\) as a sum of \(n\) rank-one outer products \[AB = \sum_{i=1}^n A^{(i)} B_{(i)}.\] Instead of summing all \(n\) terms, we pick \(c\) indices \(\{i_1, i_2, \dots, i_c\}\) from \(\{1, \dots, n\}\) with replacement, based on a probability distribution \(\{p_i\}\). We then form the approximation \[\widetilde{C} = \sum_{t=1}^c \frac{1}{c p_{i_t}} A^{(i_t)} B_{(i_t)}\] This is equivalent to \(AB \approx (AS)(S^T B)\), where \(S\) is a sparse sampling and scaling matrix.

The quality of the approximation depends heavily on the sampling probabilities \(\{p_i\}\). If we use probabilities that are proportional to the Euclidean norms of the columns and rows, \[p_i = \frac{\|A^{(i)}\| \|B_{(i)}\|}{\sum_{j=1}^n \|A^{(j)}\| \|B_{(j)}\|},\] the expected error is bounded by \[E[\|AB - \widetilde{C}\|_F] \leq \frac{1}{\sqrt{c}} \|A\|_F \|B\|_F.\] While \(\|A\|_F \|B\|_F\) can be large, this bound may be sufficient for applications where the matrices have a are already noisy or when the dimension \(n\) is in the order of millions and \(c\) is in the order of 10000, say. However, the cost for computing all the probabilities \(p_i\) must be factored into the total cost of matrix-matrix multiplication.

Random projections (like the subsampled Randomized Hadamard Transform or Gaussian sketches) lead to bounds of the form \[\|AB - (AS)(S^T B)\| \leq \|AB - (AB)_k\| + \varepsilon\|A\|_F \|B\|_F,\] where \((AB)_k\) is the best rank-\(k\) approximation of the product, see e.g. [@HMT11; @cohen2015optimal]. This indicates that the randomized product effectively captures the most significant singular values of the true product.

In conclusion, randomized sketching can approximate large matrix products efficiently when the matrices share a very large dimension. However, the benefits are limited when the output size is large, when strong spectral guarantees are required, or when the matrices have large stable rank. Consequently, sketching is most effective when the matrix product serves as an intermediate step in a larger algorithm, such as regression, low-rank approximation, or graph computations.

Exercises↩︎

You can use the mentioned fact **without proof**. Using (b), show
the following result: Let $\mathcal{X}$ be a set of $n$ points in
$\mathbb{R}^n$ and let $\varepsilon> 0$. Suppose that
$$m \geq \frac{C}{\varepsilon^2} \log n.$$

Consider a random subspace $E$ of dimension $m$ chosen uniformly
from all $m$-dimensional subspaces in $\mathbb{R}^n$, and let $P$ be
an orthogonal projection on this set. Then with probability at least
$1 - 2 \exp(-c \varepsilon^2 m)$, the scaled projection
$Q = \sqrt{\tfrac{n}{m}} P$ is an $\varepsilon$-approximate isometry
for $\mathcal{X}$, i.e., for all $x, y \in \mathcal{X}$,
$$(1 - \varepsilon) \| x - y \| \le \| Q x - Q y \| \le (1 + \varepsilon) \| x - y\|.$$
Here $C,c > 0$ are universal constants.

9 Optimization for Data Science↩︎

Optimization lies at the heart of modern data science. Whether fitting a linear regression model, clustering data, finding a sparse solution to an underdetermined system of equations, or training a deep neural network, the underlying task is almost always the same: find the parameters that minimize (or maximize) a well-defined objective function.

This chapter provides a brief, rigorous and practical introduction to optimization from a data science perspective. We begin by formalizing what an optimization problem is: a function to be minimized, variables to be chosen, and constraints that restrict which solutions are allowed. From there, we discuss convexity, introduce the necessary tools for analyzing optimality, including gradients, Hessians, the Karush–Kuhn–Tucker conditions, duality, and rate of convergence. Finally, we study gradient descent and popular variants, which play a central role in large-scale machine learning. For a comprehensive study of optimization we recommend the sources [@LVanderberghe_SBoyd_book; @calafiore2014optimization; @nesterov2013introductory; @nocedal2006numerical].

A mathematical optimization problem can be expressed as⁴⁸ \[\begin{align} \label{eq:optim} \begin{aligned} \underset{x \in \mathcal{D}}{\operatorname{minimize}} & \qquad f(x) \\ \text{subject to} & \qquad g_i(x) \le b_i, \quad \, i = 1, \dots, m, \\ & \qquad h_j(x) = c_j, \quad j=1,\dots,r. \end{aligned} \end{align}\tag{157}\] Here, the vector \(x \in {\mathcal{D}}\) is the optimization variable of the problem and \({\mathcal{D}}\) is the search space or domain of \(x\); often we set \({\mathcal{D}}= \mathbb{R}^n\). The function \(f: {\mathcal{D}}\to \mathbb{R}\) is the objective function; in data science and machine learning \(f\) is commonly referred to as the cost function or loss function when minimization is sought, and the utility function when maximization is sought. The functions \(g_i: {\mathcal{D}}\to \mathbb{R}, i=1,\dots,m\), \(h_j: {\mathcal{D}}\to \mathbb{R}, j=1,\dots,r\) are the equality and inequality constraint functions, and the constants \(b_i, c_j\) are the limits, or bounds, for the constraints, respectively.⁴⁹ If there are no constraints in 157 , the problem is called unconstrained, otherwise it falls in the realm of constrained optimization.

A vector \(x^*\) is called optimal, or a solution of the problem 157 , if it has the smallest objective value among all vectors that satisfy the constraints; that is, for any \(z\in{\mathcal{D}}\) with \(g_i(z) \le b_i, i=1,\dots,m\), \(h_j(z) = c_j, j=1,\dots,r\), we have \(f(z) \ge f(x^*)\). The set of all points satisfying the constraints is known as the feasible region, which we denote by \({\mathcal{F}}\).

Many foundational data science models naturally fit into this structure. For example, linear regression minimizes the mean squared error over all parameter vectors, usually without constraints. Logistic regression minimizes cross-entropy loss, again as an unconstrained problem unless we explicitly add regularization or parameter bounds. Data clustering and community detection can be expressed in the form 157 where the constraints are of discrete type (cf. Chapter 6 and Chapter 13, respectively). Understanding this basic formulation gives the language and framework needed for all the optimization techniques that follow.

The structure of an optimization problem determines both the difficulty of finding a solution and the type of algorithms we can use. We have already seen several optimization problems in the preceding chapter, and will see many more in chapters that follow. For example, the PCA objective 33 is an optimization problem on continuous variables, the \(k\)-means objective 89 was an optimization problem with both discrete variables (the assignments to clusters) and continuous variables (the cluster centers), and the Minimimum bisection problem 93 is an optimization problem on discrete variables. Below we will focus on continuous optimization. Indeed, while this book covers several combinatorial optimization problems the approach we take (as in Chapter 6) is often to relax to a suitable continuous optimization problems. We point the reader to [@wolsey1999integer; @schrijver2003combinatorial] for more on combinatorial optimization.

We have already seen (and will see more) optimization problems for which there are efficient algorithms (such as the PCA example) and optimization problems for which no such algorithm is believed to exist (such as minimum bisection). While the reader might be tempted to conjecture that discrete optimization problems are harder than continuous ones, the situation is more subtle: it is true that combinatorial optimization problems are often formulated over exponential sized sets (such as \(\{\pm1\}^n\)) but that on itself does not make them computationally hard (just take \(\min_{x\in\{\pm1\}^n} \mathbf{1}^Tx\) as an easy example). Also, note that while the Max-Cut problem 269 is an NP-hard combinatorial problem, it can be equivalently formulated as a continuous one by replacing \(x\in\{\pm1\}^n\) by \(-1\leq x_i\leq 1\) (see Exercise 9.1).

9.1 Convexity in Optimization↩︎

Convexity is one of the most mathematically important concepts in continuous optimization as it tends to directly impact the complexity and solvability of the problem.

Geometrically, this means that the line segment connecting any two points on the graph of the function lies entirely on or above the graph⁵⁰. We call \(f\) strictly convex if the inequality in 158 becomes a strict inequality.

If in addition \(f\) is continuously differentiable, convexity is equivalent to the first-order condition (we leave the proof of this equivalence as an exercise): \[\label{def:geomconvex} f(y) \ge f(x) + \langle \nabla f(x), y - x\rangle, \quad \forall x,y \in \mathbb{R}^d.\tag{159}\] This means that the function \(f\) is above its tangent at \(x\), as illustrated in Figure 42.

Strengthening convexity by requiring an additional positive quadratic term on the right-hand side of 159 leads to strong convexity.

This assumption will later allow us to establish exponential convergence rates.

9.1.1 Convex Optimization Problems↩︎

Recall that a set \(\mathcal{F}\) is a convex set if, for any two points \(x, y \in {\mathcal{F}}\) and any scalar \(\lambda \in [0, 1]\), the convex combination \(\lambda x + (1-\lambda) y\) also belongs to \({\mathcal{F}}\).