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) Lecture notes from these classes also contain a list of mathematical open problems [1], which has now helped to inspire a mathematical blog of open problems at https://randomstrasse101.math.ethz.ch/ . The open problems of 2024 are available at [2].. 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.

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.

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.

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.

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.

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 [3] 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 [2] (or [1] for a less updated list).

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
specialChar{39}{i}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 [4]. 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.

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 [5].

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.

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

the \( d\) -dimensional hypersphere (or \( d\) -sphere) of radius \( R\) is given by

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]\) :

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.2.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.2.1.1 Volume of the hyperball

Proposition 1

The volume of \( B^d(R)\) is given by

\[ \begin{equation} \rm{Vol}(B^d(R)) = \frac{\pi^{\frac{d}{2}} R^d}{\frac{d}{2}~\Gamma(\frac{d}{2})}. \end{equation} \]

(1)

Proof

The volume of \( B^{d}(R)\) is given by

\[ \begin{equation} \rm{Vol}(B^d(R)) = \int_0^R s_d r^{d-1} dr = \frac{s_d R^d}{d}, \end{equation} \]

(2)

where \( s_d\) denotes the (hyper-)surface area of a unit \( d\) -sphere. A unit \( d\) -sphere must satisfy

\[ s_d \int_0^{\infty} e^{-r^2} r^{d-1} dr = \underbrace{\int_{-\infty}^{\infty} \cdots \int_{-\infty}^{\infty}}_{\text{\( d\) times}} e^{-(x_1^2 + \cdots +x_d^2)} dx_1 \cdots d x_d = \Big( \int_{-\infty}^{\infty} e^{-x^2} dx \Big)^d. \]

Recall that the Gamma function is given by

\[ \Gamma(n) = \int_0^\infty r^{n-1} e^{-r} dr = 2 \int_0^\infty e^{-r^2} r^{2n-1} dr, \]

hence \( \frac{1}{2} s_d \Gamma(\frac{d}{2}) = \Big[ \Gamma(\frac{1}{2}) \Big]^d = \big( \pi^{\frac{1}{2}}\big)^d,\) and thus \( s_d = \frac{2 \pi^{\frac{d}{2}} }{\Gamma(\frac{d}{2})}.\) Plugging this expression into (2) gives

\[ \begin{equation} \rm{Vol}(B^d(R)) = \frac{\pi^{\frac{d}{2}} R^d}{\frac{d}{2}~\Gamma(\frac{d}{2})}. \end{equation} \]

(3)

Recall that for even \( d\) , we have \( \Gamma(d/2) = (d/2 - 1)!\) . To understand the asymptotic behavior \( \rm{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} \rm{Vol}(B^d) \approx \frac{1}{\sqrt{d\pi}} \Big( \frac{2\pi e}{d} \Big)^{\frac{d}{2}}, \end{align} \]

(4)

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 (4) 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}}\) .

Pan/Zoom

3.2.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 To define the “equator” of a the \( d\) -dimensional ball, we need to pick a “north pole” as reference. Without loss of generality we could pick the unit vector in the \( x_1\) -direction as defining “north”.” 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.

Fact 1

Almost all the volume of \( B^d(R)\) lies near its equator.

Proof

It suffices to prove the result for the unit \( d\) -ball. Without loss of generality we pick as “north” the direction \( x_1\) . The intersection of the sphere with the plane \( x_1=0\) forms our equator, which is formally given by the \( d-1\) -dimensional region \( \{x: \|x\| \le 1, x_1=0\}\) , where \( \|x\|\) is the Euclidean norm Throughout the book, \( \|\cdot \|\) denotes the Euclidean norm if we are dealing with vectors and it denotes the operator norm if we are dealing with matrices or operators. of \( x\) . This intersection is a sphere of dimension \( d-1\) with volume \( \rm{Vol}(B^{d-1})\) given by the \( (d-1)\) -analog of formula (3) with \( R=1\) .

We now compute the volume of \( B^d\) that lies between \( x_1=0\) and \( x_1=p_0\) . Let \( P_0=\{x: \|x\|\le 1, x_1\ge p_0\}\) be the “polar cap”, i.e., part of the sphere above the slab of width \( 2p_0\) around the equator. To compute the volume of the cap \( P\) we will integrate over all slices of the cap from \( p_0\) to 1. Each such slice will be a sphere of dimension \( d-1\) and radius \( \sqrt{1-p^2}\) , hence its volume is \( (1-p^2)^{\frac{d-1}{2}} \rm{Vol}(B^{d-1})\) . Therefore

\[ \rm{Vol}(P) = \int_{p_0}^1 (1-p^2)^{\frac{d-1}{2}} \rm{Vol}(B^{d-1}) \, dp = \rm{Vol}(B^{d-1}) \int_{p_0}^1 (1-p^2)^{\frac{d-1}{2}} \, dp. \]

Using \( e^x \ge 1+x\) for all \( x\) we can upper bound this integral by

\[ \begin{eqnarray*} \rm{Vol}(P) &\le& \rm{Vol}(B^{d-1}) \int_{p_0}^\infty e^{-\frac{d-1}{2}p^2} \, dp = \rm{Vol}(B^{d-1})\sqrt{\frac{2}{d-1}} \int_{p_0 \sqrt{\frac{d-1}{2}}}^{\infty} e^{-u^2}\,du \\ &=& \rm{Vol}(B^{d-1}) \sqrt{\frac{\pi}{2(d-1)}}\rm{erfc}\left(p_0 \sqrt{\frac{d-1}{2}} \right), \end{eqnarray*} \]

where \( \rm{erfc}(x) = \frac{2}{\sqrt{\pi}}\int_x^{\infty} e^{-u^2}\,du\) is the complementary error function. The upper bound \( \rm{erfc}(x) \leq \frac{e^{-x^2}}{\sqrt{\pi} x}\) gives

\[ \begin{equation} \rm{Vol}(P) \le \rm{Vol}(B^{d-1}) \sqrt{\frac{\pi}{2(d-1)}}\frac{e^{- \frac{d-1}{2}p_0^2}}{\sqrt{\pi} p_0 \sqrt{\frac{d-1}{2}}} = \frac{\rm{Vol}(B^{d-1})}{d-1}\frac{e^{- \frac{d-1}{2}p_0^2}}{ p_0} \end{equation} \]

(5)

Recall, from (3) that \( \rm{Vol}(B^{d}) = \frac{\pi^{\frac{d}{2}} }{\frac{d}{2} \Gamma(\frac{d}{2})}\) , so, for \( d\) large enough (since \( \frac{\Gamma(\frac{d}{2})}{\Gamma(\frac{d-1}{2})} \approx \sqrt{\frac{d}2}\) ),

\[ \begin{equation} \rm{Vol}(B^{d-1}) = \frac{\pi^{-1/2}}{\frac{d-1}{d}} \frac{\Gamma(\frac{d}{2})}{\Gamma(\frac{d-1}{2})}\rm{Vol}(B^{d}) \leq \frac{d-1}2\rm{Vol}(B^{d}). \end{equation} \]

(6)

Finally, combining (5) and (6) shows that the ratio between the volume of the polar caps and the entire hypersphere is bounded by

\[ \frac{2 \rm{Vol}(P)}{\rm{Vol}(B^{d})} \le (d-1)\frac{ \rm{Vol}(P)}{\rm{Vol}(B^{d-1})} \le \frac{1}{p_0}\exp\left(-\frac{d-1}{2}p_0^2\right). \]

The expression above shows that this ratio decreases exponentially as both \( d\) and \( p\) increase, proving our claim that the volume of the sphere concentrates strongly around its equator.

3.2.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 (3), the ratio of their volumes is

\[ \begin{align} \frac{\rm{Vol}\big(B^d(1-\varepsilon)\big)}{\rm{Vol}\big(B^d(1)\big)} = (1-\varepsilon)^d. \\\end{align} \]

(7)

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{\rm{Vol}(B^d(1-\varepsilon))}{\rm{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.2.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.

Proposition 2

The hypercube \( C^d\) has volume 1 and diameter \( \sqrt{d}\) .

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.

Pan/Zoom

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.

Pan/Zoom

Pan/Zoom

The considerations above suggest the following observation:

Fact 2

Most of the volume of the high-dimensional cube is located in its corners.

Proof

Let \( C^d\) denote the unit hypercube. The subset \( Q_R\) of \( C^d\) of points \( x\in\mathbb{R}^d\) whose distance to the center is smaller that \( R\) (\( Q_R = \{x\in C^d: \|x\|\leq R\}\) ) is given by \( Q = C^d\cap B^d(R)\) . Thus

\[ \frac{\operatorname{vol}(Q_R)}{\operatorname{vol}(C^d)} = \operatorname{vol}(Q) \le B^d(R) \approx \frac{1}{\sqrt{d\pi}} \Big( \frac{2\pi e R^2}{d} \Big)^{\frac{d}{2}}, \]

which goes to zero as long as \( R<\sqrt{\frac{d}{2\pi e}}\) . This shows that most of the volume of \( C^d\) lies in points at distance of order \( \sqrt{d}\) of the center.

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 [6, 7]).

The two most basic concepts in probability associated with a random variable \( X\) are expectation (or mean) and variance, denoted by

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

The two terms have a natural interpretation:

• \( \mathbb{E}[\rm{Var}(X \mid Y)]\) is the unexplained variance — the variability in \( X\) that persists even after knowing \( Y\) .
• \( \rm{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

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.,

with the usual extension to \( p=\infty\) by setting

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.,

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

respectively. Note that the standard deviation \( \sigma(X):=\sqrt{\rm{Var}(X)}\) of \( X\) can be written as

The covariance of the random variables \( X\) and \( Y\) is

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

The special case \( p=q=2\) is the Cauchy-Schwarz inequality

Jensen’s inequality states that for any random variable \( X\) and a convex function \( \phi: \mathbb{R} \to \mathbb{R}\) , we have

Since \( \phi(x) = x^{q/p}\) is a convex function for \( q\geq p \geq 0\) , it follows immediately from Jensen’s inequality that

Minkovskii’s inequality states that for any \( p \in [0,\infty]\) and any random variables \( X, Y\) , we have

which can be viewed as the triangle inequality.

The cumulative distribution function of \( X\) is defined by

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.3.1 Tail bounds

3.3.2 Sub-gaussian random variables

3.3.3 Sub-exponential random variables

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 [9, 10].

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. [11]). 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.4.1 Large deviation inequalities

3.4.2 The geometry of the hypercube revisited

3.4.3 Tail bounds for the chi-squared distribution

3.4.4 How to generate random points on a sphere

Hint Work with the shifted random variable \( Y = X - \mathbb{E} X + u\) , for arbitrary shift \( u\) .

Hint Note \( \mathbb{E}\left[ X \right] = \mathbb{E}\left[ X \mathbf{1}_{\left\{ X \leq \theta \, \mathbb{E}\left[ X \right] \right\}} \right] + \mathbb{E}\left[ X \mathbf{1}_{\left\{ X > \theta \, \mathbb{E}\left[ X \right] \right\}} \right]\) . How could one obtain the term \( \mathbb{E}\left[ X^2 \right]\) ?

Hint Introduce randomness on \( \varepsilon\) and take expectation.

Hint This inequality might be useful: \( 2 + t^2 \leq e^{t} + e^{-t}\) .

Hint After applying the Chernoff bound, choose \( \lambda = \frac{t}{\nu^2 + Rt}\) .

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 [15] and [16] 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

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

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

4.1.2 Existence of the 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

\[ \begin{equation} \|A - A_s \|_2 = \inf_{B \in \mathbb{R}^{m\times n}, \operatorname{rank}(B) \leq s} \|A - B \|_2 = \sigma_{s+1}. \end{equation} \]

(42)

A similar result holds for the Frobenius norm, with the approximation error given in terms of the remaining \( r-s\) smallest singular values as

\[ \begin{equation} \|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}. \end{equation} \]

(43)

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 [17].

We prove the low rank approximation property for the operator norm case Both, Theorem 9 and Theorem 11 are also known as Eckart-Young-Mirsky Theorem [16]..

Theorem 9

Let the SVD of a matrix \( A\in \mathbb{R}^{m\times n}\) be \( A = U \Sigma V^\top \) . If \( s < r = \operatorname{rank}(A)\) and we have the rank-\( s\) matrix \( A_s\) ,

\[ A_s = \sum_{k=1}^s \sigma_k u_kv_k^\top \]

then

\[ \min_{\operatorname{rank}(B) \leq s} \|A - B \|_2 = \|A - A_s \|_2 = \sigma_{s+1}. \]

Proof

Since \( A_s\) and \( A\) have the same left and right singular vectors and the same \( s\) leading singular values, it follows that \( U^\top A_s V = \operatorname{diag}(\sigma_1,…,\sigma_s,0,…,0)\) and \( U^\top (A - A_s) V = \operatorname{diag}(0,…,0,\sigma_{s+1},…,\sigma_p)\) , where \( p=\min\{m,n\}\) . This implies \( \|A-A_s\|_2 = \sigma_{s+1}\) , the largest remaining singular value. Now, suppose we have a matrix \( B \in \mathbb{R}^{m\times n}\) with \( \operatorname{rank}(B) = s\) . We can find an orthonormal basis for the null space of \( B\) , i.e. orthonormal vectors \( x_1, …, x_{n-s}\) such that \( \operatorname{null}(B) = \operatorname{span}\{x_1,…,x_{n-s}\}\) . By a dimension argument, we have that

\[ \operatorname{span}\{x_1,\ldots,x_{n-s}\}\cap \operatorname{span}\{v_1,\ldots,v_{s+1}\} \neq \{0 \}, \]

since these are two subspaces of \( \mathbb{R}^n\) of dimension \( n-s\) and \( s+1\) . Take a unit vector \( z\) in this intersection; for this vector, we know that \( Bz = 0\) and \( Az = \sum_{k=1}^{s+1} \sigma_k u_k (v_k^\top z)\) . Therefore,

\[ \|A-B\|_2^2 \geq \|(A-B)z\|^2 = \|Az\|^2 = \sum_{k=1}^{s+1} \sigma_k^2 (v_k^\top z)^2 \geq \sigma_{s+1}^2 \sum_{k=1}^{s+1} (v_k^\top z)^2 = \sigma_{s+1}^2 \]

which completes the proof.

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

Theorem 10

If \( A\) and \( B\) are complex-valued \( n\times n\) matrices with singular values \( \sigma_1(A),…, \sigma_n(A)\) and \( \sigma_1(B),…, \sigma_n(B)\) , then

\[ |\operatorname{Tr}(AB)| \leq \sum_{k=1}^n \sigma_k(A) \sigma_k(B). \]

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

Theorem 11

Let the SVD of a matrix \( A\in \mathbb{R}^{m\times n}\) be \( A = U \Sigma V^\top \) . If \( s < r = \operatorname{rank}(A)\) and we have the rank-\( s\) matrix \( A_s\) ,

\[ A_s = \sum_{k=1}^s \sigma_k u_k v_k^\top \]

then

\[ \min_{\operatorname{rank}(B) \leq s} \|A - B \|_F = \|A - A_s \|_F = \sqrt{\sum_{k=s+1}^r \sigma_k^2}. \]

Proof

Let \( B\) be any matrix with \( \operatorname{rank}(B) = s < r = \operatorname{rank}(A)\) . Notice that

\[ \|A - B\|_F^2 = \|A\|_F^2 + \|B\|_F^2 - 2\operatorname{Tr}(A^\top B) = \sum_{k=1}^r \sigma_k^2(A) + \sum_{k=1}^s \sigma_k^2(B) - 2\operatorname{Tr}(A^\top B). \]

By von Neumann’s trace inequality (whose statement for square matrices implies the general case of rectangular matrices)

\[ \operatorname{Tr}(A^\top B) \leq \sum_{k=1}^s \sigma_k(A) \sigma_k(B). \]

Therefore,

\[ \begin{eqnarray*} \|A - B\|_F^2 &\geq& \sum_{k=1}^r \sigma_k^2(A) + \sum_{k=1}^s \sigma_k^2(B) - 2\sum_{k=1}^s \sigma_k(A) \sigma_k(B) \\ &=& \sum_{k=s+1}^r \sigma_k^2(A) + \sum_{k=1}^s (\sigma_k(A) - \sigma_k(B))^2 \\ &\geq& \sum_{k=s+1}^r \sigma_k^2(A), \end{eqnarray*} \]

which establishes the lower bound. Finally, the lower bound is clearly attained for \( B = A_s\) .

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.

Pan/Zoom

Pan/Zoom

Pan/Zoom

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 9. 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.

Pan/Zoom

Remark 8

The computational complexity of computing the SVD of a matrix of size \( m\times n\) with \( m \geq n\) is \( \mathcal{O}(mn^2)\) . This cubic scaling could be prohibitive for massive data matrices. In Chapter 8.2 we discuss numerical algorithms that use randomization to efficiently compute the low rank approximation of such large matrices more efficiently.

Remark 9

The SVD plays a critical role in sensitivity analysis of linear systems of the form \( Ax=b\) . In particular, the error is bounded in terms of the condition number \( \kappa(A)\) which is the ratio of the largest singular value and the smallest singular value. A rigorous treatment of this important topic is beyond the scope of this book and can be found in other texts on numerical linear algebra [16].

4.1.4 Spectral Decomposition

4.1.5 Quadratic forms

4.1.6 SVD and the Moore-Penrose Pseudoinverse

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 [22] 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\‘eve 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

and its sample covariance as

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

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

4.2.3 Finding the Principal Components

4.2.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) [26]. 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 See [27] for more information about MNIST. The book also contains a lot of useful information about other benchmark datasets..

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 10.

Pan/Zoom

4.2.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 11(a). Figure 11(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.

Pan/Zoom

Pan/Zoom

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

Pan/Zoom

Pan/Zoom

Pan/Zoom

Pan/Zoom

Pan/Zoom

Pan/Zoom

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.

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

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\) . In this case, \( S_n\) is actually the maximum likelihood estimator for \( \Sigma\) .

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 21 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 (60), that we will discuss below.

Pan/Zoom

Pan/Zoom

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. [28, 29, 30]). This particular limiting distribution was first established in 1967 by Marčenko and Pastur [31] 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

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

\[ \begin{equation} x\sim \mathcal{N}(0,I + \beta uu^\top ). \end{equation} \]

(61)

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 (61). 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 24 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\) :

Pan/Zoom

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 25 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.

Pan/Zoom

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é [33]). There are many very nice papers about this and similar phenomena, including [34, 35, 33, 36, 37, 38, 39, 40]. Notice that the Marčenko-Pastur theorem does not imply that all eigenvalues are actually in the support of the Marčenko-Pastur distribution, it just rules out that a non-vanishing proportion are. However, it is possible to show that indeed, in the limit, all eigenvalues will be in the support (see, for example, [34]).

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 [34]) 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 = \frac1n\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 Intuitively, \( \lambda\) is larger than all the eigenvalues of \( Z_n\) , because it corresponds to a perturbation of \( Z_n\) by a positive definite matrix \( \beta uu^\top \) ; yet, a formal justification is beyond the present discussion.

\[ \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

\[ \begin{equation} u^\top (\lambda I - Z_n)^{-1} Z_n u = \frac{1}{\beta}. \end{equation} \]

(62)

Suppose \( w_1,…,w_p\) are orthonormal eigenvectors of \( Z_n\) (with corresponding eigenvalues \( \lambda_1,…, \lambda_p\) ), and expand \( u\) in that basis:

\[ u = \sum_{i=1}^p \alpha_i w_i. \]

Plugging this expansion in (62) gives

\[ \begin{equation} \sum_{i=1}^p \frac{\lambda_i}{\lambda-\lambda_i} \alpha_i^2 = \frac{1}{\beta} \end{equation} \]

(63)

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 (63) becomes

\[ \begin{equation} \lim_{p\to \infty} \frac{1}{p}\sum_{i=1}^p \frac{\lambda_i}{\lambda-\lambda_i} = \frac{1}{\beta} \end{equation} \]

(64)

Since the eigenvalues \( \lambda_1,…,\lambda_p\) follow the Marčenko-Pastur distribution, the limit on the left hand side can be replaced by the integral

\[ \begin{equation} \int_{\gamma_{-}}^{\gamma_{+}} \frac{t }{\lambda-t}\,d F_\gamma(t) = \frac{1}{\beta} \end{equation} \]

(65)

Using an integral table (or an integral software), we find that

\[ \begin{equation} \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]. \end{equation} \]

(66)

For \( \lambda=\gamma_{+}\) , that is, when the top eigenvalue touches the right edge of the Marčenko-Pastur distribution, (66) 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

\[ \begin{equation} \frac{1}{\beta_c} = \frac{1}{4 \gamma} \left((1+\sqrt{\gamma})^2-(1-\sqrt{\gamma})^2\right). \end{equation} \]

(67)

Therefore, the critical SNR is

\[ \begin{equation} \beta_{c}= \sqrt{\gamma}=\sqrt{\frac{p}{n}}. \end{equation} \]

(68)

When \( \beta > \sqrt{\frac{p}{n}}\) one can observe the pop out of the top eigenvalue from the bulk.

Eq. (68) 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 (66) (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

\[ \begin{equation} \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. \end{equation} \]

(69)

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

\[ \begin{equation} |\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. \end{equation} \]

(70)

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

4.3.3 Rank and covariance estimation

Hint Any \( V\) with \( \dim(V) = n-k+1\) intersects with the subspace spanned by the top \( k\) eigenvectors of \( A\) .

\begin{myexercise} This exercise is about PCA. We use the classical iris data set as an example (see https://en.wikipedia.org/wiki/Iris_flower_data_set for an explanation of that dataset). You can download a subset of the iris dataset from the UCI Machine Learning Repository).

The dataset contains fifty flowers and for each flower there are four measurements: sepal length, sepal width, petal length and petal width. We want to understand via PCA how these four measurements are related to each other. Are they essentially independent of each other or are some of them (weakly or strongly) correlated?

(a) Assume you loaded the iris data in Matlab into a data matrix \( X\) of size \( 4 \times 50\) . You can examine the structure of the data using plotmatrix:

[~,ax]=plotmatrix(X');
ax(1,1).YLabel.String='Sepal L';
ax(2,1).YLabel.String='Sepal W';
ax(3,1).YLabel.String='Petal L';
ax(4,1).YLabel.String='Petal W';
ax(4,1).XLabel.String='Sepal L';
ax(4,2).XLabel.String='Sepal W';
ax(4,3).XLabel.String='Petal L';
ax(4,4).XLabel.String='Petal W';

Which measurements are related and which ones are not-related? Compare your findings based on this visual inspection with your findings based on computing the variance and covariance of the data.

(b) Next carry out a principal component analysis by projecting the data onto the eigenvectors of the (centered) covariance matrix of the data. How much of the variation in the data is captured by the individual components? \end{myexercise}

Linear regression is the foundational supervised learning model in data science, statistics, and applied mathematics [53, 54, 16]. 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

This proposal immediately raises the question, why did we use the least-squares error in (73) 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,…,n} |y_i - x_i^\top\theta|\) .

The choice of the squared loss (73) is not arbitrary. It emerges from multiple complementary perspectives: geometric interpretation, computational tractability, and statistical optimality. We examine each justification in turn The dominance of the least squares approach is also deeply rooted in the history of astronomy and the development of the normal distribution. While Adrien-Marie Legendre was the first to publish the method in 1805 (as a purely algebraic way to fit curves to planetary data), Carl Friedrich Gauss published his version in 1809 and famously claimed he had been using it since 1795 [55]. Gauss provided the first probabilistic justification for the method, effectively inventing the field of statistical inference as we know it today [56]..

Pan/Zoom

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 {Total Least Squares}, a generalization of least squares, is a regression technique that assumes errors exist in {both} \( X\) and \( y\) , see e.g. [16]. The total least squares solution is equivalent to finding the first principal component of the combined matrix \( [X \, y]\) . 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\) ?”

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

Furthemore, \( s\) is the unique minimizer of

Let \( \rm{Col}(X)\) denote the column space of \( X\) and \( P_{\rm{Col}(X)}\) the orthogonal projection onto \( \rm{Col}(X)\) . Applying the subspace decomposition above to the case \( S = \rm{Col}(X)\) implies that for any \( y \in \mathbb{R}^n\) , there exists \( \hat \theta \in \mathbb{R}^p\) such that

Moreover, we have

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

if and only if

Equivalently,

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

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 \( \rm{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 may be computationally too expensive for large-scale problems. Conjugate gradient based iterative solvers are often the algorithms of choice in that case [16, 57]. See also the randomized SVD in Chapter 8.2 as a way to reduce the computational complexity of the SVD at the cost of reduced accuracy..

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 [16]. 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.

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: \( \rm{Var}(\varepsilon | X) = \sigma^2 I\) .
3. (A3) Full rank: \( \rm{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 Here, a causal interpretation means if we were to intervene and change \( x_j\) while holding other relevant factors fixed, \( y\) would change by \( \theta_j\) (on average)..

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 This is clearly a very idealized assumption. A more realistic assumption is heteroscedasticity. For example, homoscedasticity would imply that the variance of income is the same for people with 4 years of education as for people with 16 years. In contrast, heteroscedasticity would indicate that the income variance increases with education level, which is a much more realistic viewpoint since in general highly educated people have more diverse career outcomes.. 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 (73) 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)\) ,

The Maximum Likelihood Estimator (MLE) \( \hat{\theta}_{MLE}\) is defined as the argument that maximizes the likelihood function, i.e.,

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 [58, 59], i.e.

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

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,

Except for a few cases (like least squares, see also Theorem 13 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 (80) is biased; the unbiased estimator uses \( n-p\) instead of \( n\) in the denominator (see also Exercise 41).

Theorem 14 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?

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 \( \rm{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 17 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 17 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 (82) in Corollary 2 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\) ) For example if \( X\) is a Gaussian random matrix, we will see in Chapter 14 that the condition number of \( X\) grows quickly as \( p \to n\) . Indeed, let \( \gamma = p/n\) , then we have that \( \kappa(X) \approx \frac{\sqrt{n}+\sqrt{p}}{\sqrt{n}-\sqrt{p}}\) . See Theorem62 for a precise statement.. 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 2 also forms the basis for constructing prediction intervals, see e.g. [60, 53] 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

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 (84) 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 [61, 53].

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

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 (82) 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.

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 The term was introduced by Hoerl and Kennard in 1970, who used the word “ridge” specifically to refer to this flat elongated trough in the least squares objective surface that arises under multi-collinearity. They called their technique ridge regression because it was designed to deal with that ridge. 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 [62, 53]. In numerical analysis and optimization ridge regression is also known as Tikhonov regularization (named for Andrey Tikhonov [63]).

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 (89) with the least squares solution (when \( r = p\) ):

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 42).

Unlike in Theorem 15, 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 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 27 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 27, 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.

Pan/Zoom

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 [64] and can therefore be used to select hyperparameters such as the ridge penalty or the lasso regularization parameter (see (91)) 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 [65] 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 [66]).

Stein’s SURE has its roots in his earlier work [67, 68]. 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 [67] 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 [68]—that achieves a lower MSE for any true parameter value by pulling individual estimates toward a common average or the origin [69].

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) [70, 53], which replaces the \( \ell_2\) -penalty with an \( \ell_1\) -penalty, namely we consider the regularized regression problem

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

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 [71]. The Huber loss is defined as

where \( r_i = y_i - x_i^\top \theta\) is the \( i\) -th residual and \( \delta > 0\) is a tuning parameter (see also Figure 28(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 28 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.

Pan/Zoom

Pan/Zoom

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 [72, 73].

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

To find the minimum mean squared error (MMSE) estimator, we expand the cost function (92), use the fact the signal and the noise are uncorrelated, and set the derivative of (92) with respect to \( G\) to zero. We arrive at the solution, also called Wiener Gain,

where \( \text{SNR} = \sigma_s^2 / \sigma_n^2\) represents the signal-to-noise ratio Because many signals have a very wide dynamic range, signals are often expressed using the logarithmic decibel scale. Thus, in practice the SNR is often expressed in decibels as \( \text{SNR}_{\text{dB}} =10\log_{10}(\text{SNR})\) .. 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

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

We place a zero-mean isotropic Gaussian prior on \( \theta\) :

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

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 (96) and set \( \nabla_\theta \mathcal{L}(\theta) = 0\) ; this gives

Rearranging terms and multiplying by \( \sigma^2\) yields

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

We already confirmed earlier that this gives a strict global minimum.

Because both the likelihood (94) and the prior (95) are Gaussian, the posterior is also Gaussian:

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 (89). The ridge estimate takes the form

Now we apply the rotations

In these coordinates the model (93) 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 (102) decouples into \( p\) independent scalar problems:

We now introduce the per-component signal-to-noise ratio. Under the prior (95), 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

Substituting \( \lambda = \sigma^2/\tau^2\) , note that

so the shrinkage factor in (104) can be written as

The ridge estimate in the principal-component basis therefore takes the form

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:

Dividing numerator and denominator of the second term on the right-hand side of (109) by \( \sigma^2\) gives

which is identical to (108). 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, ~~ \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, ~~ \mathrm{SNR}_i = \frac{\sigma_i^2\tau^2}{\sigma^2} = \frac{\sigma_i^2}{\lambda}. \]

\begin{myexercise} This exercise is about ridge regression and cross-validation. Use the following Matlab-type code to generate a synthetic high-dimensional dataset where only a few features are truly informative:

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);

1. Recall that the ridge estimator is given by \( \hat{\theta}_{\lambda} = (X^T X + \lambda I)^{-1} X^T y\) . Write a function myridge(X, y, lambda) that returns the regression coefficients. Efficiency Hint: Since \( p > n\) , it is computationally cheaper to use the identity \( (X^T X + \lambda I)^{-1} X^T = X^T (XX^T + \lambda I)^{-1}\) . Prove this identity in your write-up.
2. Implement a 5-fold cross-validation routine to select the optimal \( \lambda\) from an equally spaced grid: lambdas = linspace(0.1, 5, 50). Partition the \( n=500\) observations into five folds of equal size. For each \( \lambda\) in the grid: Apply the ridge estimator to the four folds and calculate the Mean Squared Error (MSE) on the remaining fold. Average the MSE across all five folds to get the Cross-Validation Error, \( CV(\lambda)\) . Plot the \( CV(\lambda)\) curve against \( \log(\lambda)\) . Mark the \( \lambda^*\) that minimizes the error.
3. The “One-Standard-Error” Rule: Often, the \( \lambda\) that strictly minimizes \( CV(\lambda)\) is too small. Identify the largest \( \lambda\) whose error is within one standard error of the minimum. Why might a researcher prefer this “simpler” model?
4. Compute the SVD of \( X = U\Sigma V^T\) . Explain how \( \lambda\) modifies the singular values \( \sigma_i\) in the solution. Specifically, what happens to the components of the solution corresponding to very small \( \sigma_i\) as \( \lambda\) increases?

\end{myexercise}

Before we introduce the formalism of graph theory, we describe the celebrated PageRank algorithm. This algorithm is a principal component It is difficult to resist using this pun. 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: Another important component of modern search engines is personalization, which we do not discuss here.

• 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 [74]. 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 31, 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\) .

Pan/Zoom

A very simple approach is to take \( x_k\) as the number of backlinks for page \( k\) . In the example in Figure 31, 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 31, 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

where \( n_j\) is the number of outgoing links from page \( j\) .

If we apply these scheme to the toy example in Figure 31, 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

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. [15] 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, One such modification is the addition of very small weight edges between every pair of pages, corresponding to a “random teleport” of the random walker (or web crawler) that happens with very low probability. One advantage this brings is that the networks becomes irreducible. A reducible directed graph is one for which there are nodes \( u\) and \( v\) for which there is no directed path from \( u\) to \( v\) , if a networks is not reducible, it is called irreducible. 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 [75] (with a follow-up later [76]). You can read about this nice story in [77].

We now introduce the formalism for undirected The previous Chapter featured directed graphs, in which edges (links) have a meaningful direction. In what follows we will focus in undirected graphs in which an edge represents a connection, without meaningful direction. 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 32 depicts one of the graph theorists’ favorite examples, the Petersen graph The Petersen graph is often used as a counter-example in graph theory..

Pan/Zoom

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{ :}= \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

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.

Pan/Zoom

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:

Note that, given the partition, the optimal centers are given by

Lloyd’s algorithm [78] (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 = \frac1{\left| S_l \right|} \sum_{i\in S_l}x_i\) .

Unfortunately, Lloyd’s algorithm is not guaranteed to converge to the solution of (112). Indeed, Lloyd’s algorithm oftentimes gets stuck in local optima of (112). 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 (112) 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 13.1)

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 (112) 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 (112) 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 34.

Pan/Zoom

A natural way to try to overcome the issues of \( k\) -means depicted in Figure 34 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(-\frac1{2\varepsilon}u^2\right)\) , by associating each point to a vertex and, for which pair of nodes, set the edge weight as

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

6.4.2 Normalized Cut as a spectral relaxation

Because the relaxation (120) is obtained from (119) by removing a constraint we immediately have that

This means that

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

This implies in particular that

meaning that Algorithm 1 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 [79], and the graph version is due to Noga Alon and Vitaly Milman [80, 81] 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.1 Bounding the diameter of a graph using \( \lambda_2(\mathcal{L}_G)\)

6.5.2 The proof of Cheeger’s Inequality

Multiple Clusters

Pan/Zoom

Hint (c): Expand the square in (127) at the optimal choice of centers \( \mu_l = \frac{1}{|S_l|}\sum_{i\in S_l} x_i\) .

Hint There are finitely many partitions of \( n\) points.

Hint (b): Consider the quadratic form \( x^\top L_{C_n} x\) for the vector \( x \in \mathbb{R}^n\) given by \( x_i = |i-\frac{n}{2}|-\frac{n}{4}\) .

Exercise 60 (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 36.

Pan/Zoom

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\) .]

1. Compute \( h(S_{\left\{ 1 \right\}})\) , the Cheeger’s cut of the subset \( S_{\left\{ 1 \right\}}\) .
2. Show that every vertex in \( S_T\) has exactly \( d-|T|\) neighbours in \( S_T\) .
3. Show that for any \( T \subseteq [d]\) , \( v_T\) is an eigenvector of \( \mathcal{L}_G\) with eigenvalue \( \frac{2|T|}{d}\) .
4. Show that if \( T' \subseteq [d]\) is distinct from \( T\) , then \( v_T\) and \( v_{T'}\) are orthogonal.
5. Conclude that for any \( 0 \leq k \leq d\) , the eigenspace corresponding to the eigenvalue \( \frac{2k}{d}\) has dimension \( {d \choose k}\) .
6. Compute \( h_G\) , the Cheeger’s constant of \( G\) .

Hint (d): Try similar ideas as in the proof of the Gershgorin circle theorem, using \( \lambda_n(A) = -d\) .

Diffusion maps [86] 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 The goal in the context of PageRank was to do ranking, in which case a directed network is usually best. In this Chapter we will be interested in understanding similarities between data points, in which case it is more common to make use of undirected graphs.; 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:

where \( \deg(i) = \sum_j w_{ij}\) . Let \( M\) be the matrix of these probabilities,

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

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

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

\( M\) is not symmetric, but a matrix similar to M, \( S = D^{\frac12} M D^{-\frac12}\) is, indeed \( S = D^{-\frac12}WD^{-\frac12}\) . We consider the spectral decomposition of \( S\)

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 … \geq \lambda_n\) ). Note that \( Sv_k = \lambda_k v_k\) . Also,

We define \( \Phi = D^{-\frac12} V\) with columns \( \Phi = \left[ \varphi_1,\dots,\varphi_n\right]\) and \( \Psi = D^{\frac12} V\) with columns \( \Psi = \left[ \psi_1,\dots,\psi_n\right]\) . Then

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\) :

Also, we can rewrite this decomposition as

and it is easy to see that

Let us revisit the embedding suggested on Remark 20. It would correspond to

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

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 (133) 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 [87, 88].

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 37, 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.

Pan/Zoom

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 38). 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 38 for an example). We will achieve this by creating a graph from the data points.

Pan/Zoom

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 A detailed discussion of kernel based methods will be given in Chapter 7.5. For now it suffices to think of a kernel as a function. \( 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( -\frac1{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! It is worth noting that in a very stylized example where the only possibilities are the blobs not intersecting, or intersecting by a fixed amount, it is possible that, since this effect is captured by just three values, it can still be captured linearly, and that in such simple examples PCA could still capture such low-dimensional structure.

Pan/Zoom

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 39.

Pan/Zoom

Pan/Zoom

For the two dimensional example, we expect the structure of the underlying manifold to be a two-dimensional torus. Indeed, Figure 40 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 [90, 91] 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 [92, 93].

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 43 and respectively. Figure 43 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.

Pan/Zoom

Pan/Zoom

Figure 44 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 [94], 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)\) .

Pan/Zoom

The embedding provided by the diffusion map reveals the intrinsic two-dimensional structure of the illumination manifold, see Figure 45. 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.

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 34 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(-\frac1{2\varepsilon}u^2\right)\) , by associating each point to a vertex and, for which pair of nodes, set the edge weight as

Recall the construction of a matrix \( M = D^{-1}W\) as the transition matrix of a random walk

where \( D\) is the diagonal with \( D_{ii} = \deg(i)\) . The \( d\) -dimensional diffusion map is given by

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 2, 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^{-\frac12}WD^{-\frac12}\) as defined in (130) is related to \( \mathcal{L}_G\) by

and \( \Phi = D^{-1/2}V\) .

Pan/Zoom

Proposition 15 below establishes a connection between \( \operatorname{Ncut}\) (as described in Chapter 6) and the random walks introduced above. Let \( M\) as defined in (129) 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

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\)

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 [89], 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 [95, 96]).

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

Mean First Passage Time

The expected commute time between \( {x}_{i}\) and \( {x}_{j}\) is defined as

Next, we show that the expected commute time is a metric. Let

Notice that \( \tau_{ii} \neq \tau'_{ii}\) because \( \tau'_{ii} > 0.\) Let

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, This step needs irreducibility as otherwise some of the terms would be infinite and the cancellation in the last equality wouldn't necessarily hold

that is,

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\)

Multiply (139) by \( {\psi}_{1}^{T}\) on both sides to obtain

Thus, we have

or equivalently,

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 (139). Indeed, assume that \( T_2'\) is another solution such that

then

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,

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

can be written as

Therefore, (139) is equivalent to \( T+\mathrm{ddiag}(T')=\mathbf{1}\mathbf{1}^T+MT\) which we can rewrite as

Unfortunately, \( (I-M)\) is not invertible, so we cannot obtain \( T\) by simply inverting \( (I-M)\) .

Recall that the graph Laplacian is defined as

Multiplying (147) by \( D\) gives

We also recall that \( L\mathbf{1}=0\) and \( L\succeq 0\) .

The eigen-decomposition of \( L\) can be written as

with \( 0=\mu_{1}\leq\mu_{2}\leq…\leq\mu_{n}\) , and \( \rho_1 = \frac{1}{\sqrt{n}}{\textbf {1}}\) . Note that since \( G\) is connected then \( \mu_{2}>0\) . The pseudo-inverse of the graph Laplacian is given by

Then, \( L^\dagger\) has the properties listed in (46). Moreover,

Multiplying (149) by \( L^\dagger\) gives