Chapter 2 Learning Linear and Independent Structures

In this chapter

A Low-Dimensional Subspace
Principal Components Analysis (PCA)Pursuing Low-rank Structure via Power Iteration Probabilistic PCA Matrix Completion
A Mixture of Complete Low-Dimensional Subspaces
Mixtures of Subspaces and Sparsely-Used Dictionaries Complete Dictionary Learning Connection to ICA and Kurtosis
A Mixture of Overcomplete Low-Dimensional Subspaces
Sparse Coding with an Overcomplete Dictionary Overcomplete Dictionary Learning Learned Deep Sparse Coding
Summary and Notes
Exercises and Extensions

“The art of doing mathematics consists in finding that special case which contains all the germs of generality.”

– David Hilbert

Real data has low-dimensional structure. To see why this is true, let us consider the unassuming case of static on a TV when the satellite isn’t working. At each frame (approximately every $\frac{1}{30}$ seconds), the RGB static on a screen of size $H\times W$ is, roughly, sampled independently from a uniform distribution on $[0,1]^{3\times H\times W}$ . In theory, the static could resolve to a natural image on any given frame, but even if you spend a thousand years looking at the TV screen, it will not. This discrepancy is explained by the fact that the set of $H\times W$ natural images takes up a vanishingly small fraction of the hypercube $[0,1]^{3\times H\times W}$ . In particular, it is extremely low-dimensional, compared to the ambient space dimension. Similar phenomena occur for all other types of natural data, such as text, audio, and video. Thus, when we design systems and methods to process natural data and learn its structure or distribution, this is a central property of natural data which we need to take into account.

Therefore, our central task is to learn a distribution that has low intrinsic dimension in a high-dimensional space. In the remainder of this section, we discuss several classical methods to perform this task for several somewhat idealistic models for the distribution, namely models that are geometrically linear or statistically independent. While these models and methods are important and useful in their own right, we discuss them here as they motivate, inspire, and serve as a predecessor or analogue to more modern methods for more general distributions that involve deep (representation) learning.

Our main approach (and general problem formulation) can be summarized as:

Problem: Given one or several (noisy or incomplete) observations of a ground truth sample from the data distribution, obtain an estimate of this sample.

This approach underpins several classical methods for data processing, which we discuss in this chapter.

•

Section 2.1 — Principal Components Analysis (PCA): Given noisy samples from a distribution supported on one low-dimensional subspace, obtain an estimate of the true sample that lies on this subspace.
•

Section 2.2 — Complete Dictionary Learning and Independent Components Analysis (ICA): Given noisy samples from a distribution supported on a union (not the span) of a few low-dimensional subspaces, obtain an estimate of the true samples.
•

Section 2.3 — Sparse Coding and Overcomplete Dictionary Learning: Given noisy samples from a distribution supported on combinations of a few incoherent vectors, such as the coordinate axes, obtain an estimate of the true sample, which also has this property.

As we will soon reveal in later chapters, in the deep learning era, modern methods essentially adopt the same approach to learn.

In this chapter, as described above, we make simplifying modeling assumptions that essentially assume the data have geometrically (nearly, piece-wise) linear structures and statistically independent components. In Chapter 1, we have referred to such models of the data as “analytical models”. These modeling assumptions allow us to derive efficient algorithms with provable efficiency guarantees¹¹1Both in terms of data and computational complexity. for processing data at scale. However, they are imperfect models for often-complex real-world data distributions, and so their underlying assumptions only approximately hold. This means that the guarantees afforded by detailed analysis of these algorithms also only approximately hold in the case of real data. Nonetheless, the techniques discussed in this chapter are useful in their own right, and beyond that, serve as the “special case with the germ of generality”, so to speak, in that they present a guiding motivation and intuition for the more general paradigms within (deep) learning of more general distributions that we later address.

2.1 A Low-Dimensional Subspace

2.1.1 Principal Components Analysis (PCA)

Perhaps the simplest modeling assumption possible for low-dimensional structure is the so-called low-rank assumption. Letting $D$ be the dimension of our data space, we assume that our data belong to a low-dimensional subspace of dimension $d\ll D$ , possibly plus some small disturbances. This ends up being a nearly valid assumption for some surprisingly complex data, such as images of handwritten digits and face data [RD03] as shown in Figure 2.1, yet as we will see, it will lend itself extremely well to comprehensive analysis.

Figure 2.1 : Images of human faces and handwritten digits. Despite the seemingly large variety in their appearances, each set of these data spans (approximately) a very low-dimensional (nearly) linear subspace. — Figure 2.1: Images of human faces and handwritten digits. Despite the seemingly large variety in their appearances, each set of these data spans (approximately) a very low-dimensional (nearly) linear subspace.

Problem formulation.

To write this in mathematical notation, we represent a subspace $\mathcal{S}\subseteq\mathbb{R}^{D}$ of dimension $d$ by an orthonormal matrix $\bm{U}\in\mathsf{O}(D,d)\subseteq\mathbb{R}^{D\times d}$ such that the columns of $\bm{U}$ span $\mathcal{S}$ . Then, we say that our data $\{\bm{x}_{i}\}_{i=1}^{N}\subseteq\mathbb{R}^{D}$ have (approximate) low-rank structure if there exists an orthonormal matrix $\bm{U}\in\mathsf{O}(D,d)$ , vectors $\{\bm{z}_{i}\}_{i=1}^{N}\subseteq\mathbb{R}^{d}$ , and small vectors $\{\bm{\varepsilon}_{i}\}_{i=1}^{N}\subseteq\mathbb{R}^{D}$ such that

\bm{x}_{i}=\bm{U}\bm{z}_{i}+\bm{\varepsilon}_{i},\quad\forall i\in[N].

(2.1.1)

Here $\bm{\varepsilon}_{i}$ are disturbances that prevent the data from being perfectly low rank; their presence in our model allows us to quantify the degree to which our analysis remains relevant in the presence of deviations from our model. The true supporting subspace is $\mathcal{S}\doteq\mathop{\mathrm{col}}(\bm{U})$ , the span of the columns of $\bm{U}$ . To process all we can from these data, we need to recover $\mathcal{S}$ ; to do this it is sufficient to recover $\bm{U}$ , also called the principal components. Fortunately, this is a computationally tractable task named Principal Component Analysis, and we discuss now how to solve it.

Given data $\{\bm{x}_{i}\}_{i=1}^{N}$ distributed as in (2.1.1), we aim to recover the model $\bm{U}$ . A natural approach is to find the subspace $\bm{U}^{\star}$ and latent vectors $\{\bm{z}_{i}^{\star}\}_{i=1}^{N}$ which yield the best approximation $\bm{x}_{i}\approx\bm{U}^{\star}\bm{z}_{i}^{\star}$ . Namely, we aim to solve the problem

\min_{\tilde{\bm{U}},\{\tilde{\bm{z}}_{i}\}_{i=1}^{N}}\frac{1}{N}\sum_{i=1}^{N}\|\bm{x}_{i}-\tilde{\bm{U}}\tilde{\bm{z}}_{i}\|_{2}^{2},

(2.1.2)

where $\tilde{\bm{U}}$ is constrained to be an orthonormal matrix, as above. We will omit this constraint in similar statements below for the sake of concision.

Subspace encoding-decoding via denoising.

To simplify this problem, for a fixed $\tilde{\bm{U}}$ , we have (proof as exercise):

	$\displaystyle\min_{\{\tilde{\bm{z}}_{i}\}_{i=1}^{N}}\frac{1}{N}\sum_{i=1}^{N}\\|\bm{x}_{i}-\tilde{\bm{U}}\tilde{\bm{z}}_{i}\\|_{2}^{2}$	$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\min_{\tilde{\bm{z}}_{i}}\\|\bm{x}_{i}-\tilde{\bm{U}}\tilde{\bm{z}}_{i}\\|_{2}^{2}$		(2.1.3)
		$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\\|\bm{x}_{i}-\tilde{\bm{U}}\tilde{\bm{U}}^{\top}\bm{x}_{i}\\|_{2}^{2}.$		(2.1.4)

That is, the optimal solution $(\bm{U}^{\star},\{\bm{z}_{i}^{\star}\}_{i=1}^{N})$ to the above optimization problem has $\bm{z}_{i}^{\star}=(\bm{U}^{\star})^{\top}\bm{x}_{i}$ .

Now, we can write the original optimization problem in $\tilde{\bm{U}}^{\star}$ and $\{\tilde{\bm{z}}_{i}\}_{i=1}^{N}$ as an optimization problem just over $\tilde{\bm{U}}$ , i.e., to obtain the basis $\bm{U}^{\star}$ and compact codes $\{\bm{z}_{i}^{\star}\}_{i=1}^{N}$ it suffices to solve either of the two following equivalent problems:

\min_{\tilde{\bm{U}},\{\tilde{\bm{z}}_{i}\}_{i=1}^{N}}\frac{1}{N}\sum_{i=1}^{N}\|\bm{x}_{i}-\tilde{\bm{U}}\tilde{\bm{z}}_{i}\|_{2}^{2}=\min_{\tilde{\bm{U}}}\frac{1}{N}\sum_{i=1}^{N}\|\bm{x}_{i}-\tilde{\bm{U}}\tilde{\bm{U}}^{\top}\bm{x}_{i}\|_{2}^{2}.

(2.1.5)

Note that the problem on the right-hand side of (2.1.5) is a denoising problem: given noisy observations $\bm{x}_{i}$ of low-rank data, we aim to find the noise-free copy of $\bm{x}_{i}$ , which under the model (2.1.1) is $\bm{U}\bm{z}_{i}$ . That is, the denoised input $\hat{\bm{x}}_{i}=\bm{U}\bm{U}^{\top}\bm{x}_{i}$ . Notice that this is the point on the subspace that is closest to $\bm{x}_{i}$ , as visualized in Figure 2.2. Here by solving the equivalent problem of finding the best subspace, parameterized by the learned basis $\bm{U}^{\star}$ , we learn an approximation to the denoiser, i.e., the projection matrix $\bm{U}^{\star}(\bm{U}^{\star})^{\top}\approx\bm{U}\bm{U}^{\top}$ that projects the noisy data point $\bm{x}_{i}$ onto the subspace $\mathcal{S}$ .

Figure 2.2: Geometry of PCA. A data point

\bm{x}

(red) is projected onto the one-dimensional learned subspace spanned by the unit basis vector

\bm{u}_{1}

(blue arrow). The projection

\bm{U}\bm{U}^{\top}\bm{x}=\bm{u}_{1}\bm{u}_{1}^{\top}\bm{x}

(green) is the denoised version of

\bm{x}

using the low-dimensional structure given by

\bm{u}_{1}

, and

\bm{\varepsilon}

(brown arrow) represents the projection residual or noise.

Putting the above process together, we essentially obtain a simple encoding-decoding scheme that maps a data point $\bm{x}$ in $\mathbb{R}^{D}$ to a lower-dimensional (latent) space $\mathbb{R}^{d}$ and then back to $\mathbb{R}^{D}$ :

\bm{x}\xrightarrow{\hskip 5.69054pt\mathcal{E}=(\bm{U}^{\star})^{\top}\hskip 5.69054pt}\bm{z}\xrightarrow{\hskip 5.69054pt\mathcal{D}=\bm{U}^{\star}\hskip 5.69054pt}\hat{\bm{x}}.

(2.1.6)

Here, $\bm{z}\in\mathbb{R}^{d}$ can be viewed as the low-dimensional compact code (or a latent representation) of a data point $\bm{x}\in\mathbb{R}^{D}$ and the learned subspace basis $\bm{U}^{\star}$ as the associated codebook whose columns are the (learned) optimal code words. The process achieves the function of denoising $\bm{x}$ by projecting it onto the subspace spanned by $\bm{U}^{\star}$ .

Computing the subspace basis.

For now, we continue our calculation. Let $\bm{X}=\begin{bmatrix}\bm{x}_{1},\dots,\bm{x}_{N}\end{bmatrix}\in\mathbb{R}^{D\times N}$ be the matrix whose columns are the observations $\bm{x}_{i}$ . We have (proof as exercise)

	$\displaystyle\operatorname*{arg\ min}_{\tilde{\bm{U}}}\frac{1}{N}\sum_{i=1}^{N}\\|\bm{x}_{i}-\tilde{\bm{U}}\tilde{\bm{U}}^{\top}\bm{x}_{i}\\|_{2}^{2}$	$\displaystyle=\operatorname*{arg\ max}_{\tilde{\bm{U}}}\frac{1}{N}\sum_{i=1}^{N}\\|\tilde{\bm{U}}^{\top}\bm{x}_{i}\\|_{F}^{2}$		(2.1.7)
		$\displaystyle=\operatorname*{arg\ max}_{\tilde{\bm{U}}}\operatorname{tr}\left\{\tilde{\bm{U}}^{\top}\left(\frac{\bm{X}\bm{X}^{\top}}{N}\right)\tilde{\bm{U}}\right\}.$		(2.1.8)

Thus, to compute the principal components, we find the orthogonal matrix $\tilde{\bm{U}}$ which maximizes the term $\operatorname{tr}(\tilde{\bm{U}}^{\top}(\bm{X}\bm{X}^{\top}/N)\tilde{\bm{U}})$ . We can prove via induction that this matrix $\bm{U}^{\star}$ has columns which are the top $d$ unit eigenvectors of $\bm{X}\bm{X}^{\top}/N$ . We leave the whole proof to the reader in Exercise 2.1, but we handle the base case of the induction here. Suppose that $d=1$ . Then we only have a single unit vector $\bm{u}$ to recover, so the above problem reduces to

\max_{\tilde{\bm{u}}\colon\|\tilde{\bm{u}}\|_{2}=1}\tilde{\bm{u}}^{\top}(\bm{X}\bm{X}^{\top}/N)\tilde{\bm{u}}.

(2.1.9)

This is the so-called Rayleigh quotient of $\bm{X}\bm{X}^{\top}/N$ . By invoking the spectral theorem we diagonalize $\bm{X}\bm{X}^{\top}/N=\bm{V}\bm{\Lambda}\bm{V}^{\top}$ , where $\bm{V}$ is orthogonal and $\bm{\Lambda}$ is diagonal with non-negative entries. Hence

\tilde{\bm{u}}^{\top}(\bm{X}\bm{X}^{\top}/N)\tilde{\bm{u}}=\tilde{\bm{u}}^{\top}\bm{V}\bm{\Lambda}\bm{V}^{\top}\bm{u}=(\bm{V}^{\top}\tilde{\bm{u}})^{\top}\bm{\Lambda}(\bm{V}^{\top}\tilde{\bm{u}}).

(2.1.10)

Since $\bm{V}$ is an invertible orthogonal transformation, $\bm{V}^{\top}\tilde{\bm{u}}$ is a unit vector, and optimizing over $\tilde{\bm{u}}$ is equivalent to optimizing over $\tilde{\bm{w}}\doteq\bm{V}^{\top}\tilde{\bm{u}}$ . Hence, we can write

\tilde{\bm{u}}^{\top}(\bm{X}\bm{X}^{\top}/N)\tilde{\bm{u}}=\tilde{\bm{w}}^{\top}\bm{\Lambda}\tilde{\bm{w}},

(2.1.11)

whose optimal solutions $\bm{w}^{\star}$ among unit vectors are one-hot vectors whose only nonzero (hence unit) entry is in one of the indices corresponding to the largest eigenvalue of $\bm{X}\bm{X}^{\top}/N$ . This means that $\tilde{\bm{u}}=\bm{V}\tilde{\bm{w}}$ , the optimal solution to the original problem, corresponds to a unit eigenvector of $\bm{X}\bm{X}^{\top}/N$ (i.e., column of $\bm{V}$ ) which corresponds to the largest eigenvalue. Suitably generalizing this to the case $d>1$ , and summarizing all the previous discussion, we have the following informal Theorem.

Theorem 2.1.

Suppose that our dataset $\{\bm{x}_{i}\}_{i=1}^{N}\subseteq\mathbb{R}^{D}$ can be written as

\bm{x}_{i}=\bm{U}\bm{z}_{i}+\bm{\varepsilon}_{i},\qquad\forall i\in[N],

(2.1.12)

where $\bm{U}\in\mathsf{O}(D,d)$ captures the low-rank structure, $\{\bm{z}_{i}\}_{i=1}^{N}\subseteq\mathbb{R}^{d}$ are the compact codes of the data, and $\{\bm{\varepsilon}_{i}\}_{i=1}^{N}\subseteq\mathbb{R}^{D}$ are small vectors indicating the deviation of our data from the low-rank model. Then the principal components $\bm{U}^{\star}\in\mathsf{O}(D,d)$ of our dataset are given by the top $d$ eigenvectors of $\bm{X}\bm{X}^{\top}/N$ , where $\bm{X}=[\bm{x}_{1},\dots,\bm{x}_{N}]\in\mathbb{R}^{D\times N}$ , and approximately correspond to the optimal linear denoiser: $\bm{U}^{\star}(\bm{U}^{\star})^{\top}\approx\bm{U}\bm{U}^{\top}$ .

We do not give explicit rates of approximation here as they can become rather technical. In the special case that $\bm{\varepsilon}_{i}=\bm{0}$ for all $i$ , the learned $\bm{U}^{\star}$ spans the support of the samples $\{\bm{x}_{i}\}_{i=1}^{N}$ . If in addition the $\bm{z}_{i}$ are sufficiently diverse (say, spanning all of $\mathbb{R}^{d}$ ) then we would have perfect recovery: $\bm{U}^{\star}(\bm{U}^{\star})^{\top}=\bm{U}\bm{U}^{\top}$ .

Remark 2.1.

In some data analysis tasks, the data matrix $\bm{X}$ is formatted such that each data point is a row rather than a column as is presented here. In this case the principal components are the top $d$ eigenvectors of $\bm{X}^{\top}\bm{X}/N$ .

Remark 2.2 (Basis Selection via Denoising Eigenvalues).

In many cases, either our data will not truly be distributed according to a subspace-plus-noise model, or we will not know the true underlying dimension $d$ of the subspace. In this case, we have to choose $d$ ; this problem is called model selection. In the restricted case of PCA, one way to perform model selection is to compute $\bm{X}\bm{X}^{\top}/N$ and look for instances where adjacent eigenvalues sharply decrease; this is one indicator that the index of the larger eigenvalue is the “true dimension $d$ ”, and the rest of the eigenvalues of $\bm{X}\bm{X}^{\top}/N$ are contributed by the noise or disturbances $\bm{\varepsilon}_{i}$ . Model selection is a difficult problem and, nowadays in the era of deep learning where it is called “hyperparameter optimization”, is usually done via brute force or Bayesian optimization.

Remark 2.3 (Denoising Samples).

The expression on the right-hand side of (2.1.5), that is,

\min_{\tilde{\bm{U}}}\frac{1}{N}\sum_{i=1}^{N}\|\bm{x}_{i}-\tilde{\bm{U}}\tilde{\bm{U}}^{\top}\bm{x}_{i}\|_{2}^{2},

(2.1.13)

is what’s known as a denoising problem, thusly named because it is an optimization problem whose solution removes the noise from the samples so it fits on the subspace. Denoising—learning a map which removes noise from noisy samples so that it fits on the data structure (such as in (2.1.1), but maybe more complicated)—is a common method for learning distributions that will be discussed in the sequel and throughout the manuscript. Note that we have already discussed this notion, but it bears repeating due to its central importance in later chapters.

Remark 2.4 (Neural Network Interpretation).

If we do a PCA, we approximately recover the support of the distribution encoded by the parameter $\bm{U}^{\star}$ . The learned denoising map then takes the form $\bm{U}^{\star}(\bm{U}^{\star})^{\top}$ . On top of being a denoiser, this can be viewed as a simple two-layer weight-tied linear neural network: the first layer multiplies by $(\bm{U}^{\star})^{\top}$ , and the second layer multiplies by $\bm{U}^{\star}$ , namely

\operatorname{denoise}(\bm{x})=\underbrace{\bm{U}^{\star}\circ\underbrace{\operatorname{id}\circ\underbrace{(\bm{U}^{\star})^{\top}\bm{x}}_{\text{first ``layer''}}}_{\text{post-activation of first ``layer''}}}_{\text{output of ``NN''}}

(2.1.14)

Contrasting this to a standard two-layer neural network, we see a structural similarity:

\operatorname{NN}(\bm{x})=\underbrace{\bm{W}^{\star}\circ\underbrace{\mathrm{ReLU}\circ\underbrace{(\bm{U}^{\star})^{\top}\bm{x}}_{\text{first layer}}}_{\text{post-activation of first layer}}}_{\text{output of NN}}

(2.1.15)

In particular, PCA can be interpreted as learning a simple two-layer denoising autoencoder,²²2In fact, as we have mentioned in the previous chapter, PCA was one of the first problems that neural networks were used to solve [Oja82, BH89]. one of the simplest examples of a non-trivial neural network. In this framework, the learned representations would just be $(\bm{U}^{\star})^{\top}\bm{x}\approx\bm{z}$ . In this way, PCA serves as a model problem for (deep) representation learning, which we will draw upon further in the monograph. Notice that in this analogy, the representations reflect, or are projections of, the input data towards a learned low-dimensional structure. This property will be particularly relevant in the future.

2.1.2 Pursuing Low-rank Structure via Power Iteration

There is a computationally efficient way to estimate the top eigenvectors of $\bm{X}\bm{X}^{\top}/N$ or any symmetric positive semidefinite matrix $\bm{M}$ , called power iteration. This method is the building block of several algorithmic approaches to high-dimensional data analysis that we discuss later in the chapter, so we discuss it here.

Let $\bm{M}$ be a symmetric positive semidefinite matrix. There exists an orthonormal basis for $\mathbb{R}^{D}$ consisting of eigenvectors $(\bm{w}_{i})_{i=1}^{D}$ of $\bm{M}$ , with corresponding eigenvalues $\lambda_{1}\geq\cdots\geq\lambda_{D}\geq 0$ . By definition, any eigenvector $\bm{w}_{i}$ satisfies $\lambda_{i}\bm{w}_{i}=\bm{M}\bm{w}_{i}$ . Therefore, for any $\lambda_{i}>0$ , $\bm{w}_{i}$ is a “fixed point” to the following equation:

\bm{w}=\frac{\bm{M}\bm{w}}{\|\bm{M}\bm{w}\|_{2}}.

(2.1.16)

Theorem 2.2 (Power Iteration).

Assume that $\lambda_{1}>\lambda_{i}$ for all $i>1$ . If we compute the fixed point of (2.1.16) using the following iteration:

\bm{v}_{0}\sim\operatorname{\mathcal{N}}(\bm{0},\bm{1}),\qquad\bm{v}_{t+1}\leftarrow\frac{\bm{M}\bm{v}_{t}}{\|\bm{M}\bm{v}_{t}\|_{2}},

(2.1.17)

then, in the limit, $\bm{v}_{t}$ will converge to a top unit-norm eigenvector of $\bm{M}$ .

Proof.

First, note that for all $t$ , we have

\bm{v}_{t}=\frac{\bm{M}\bm{v}_{t-1}}{\|\bm{M}\bm{v}_{t-1}\|_{2}}=\frac{\bm{M}^{2}\bm{v}_{t-2}}{\|\bm{M}\bm{v}_{t-1}\|_{2}\|\bm{M}\bm{v}_{t-2}\|_{2}}=\cdots=\frac{\bm{M}^{t}\bm{v}_{0}}{\prod_{s=1}^{t}\|\bm{M}\bm{v}_{s}\|_{2}}.

(2.1.18)

Thus, $\bm{v}_{t}$ has the same direction as $\bm{M}^{t}\bm{v}_{0}$ and is unit norm, so we can write

\bm{v}_{t}=\frac{\bm{M}^{t}\bm{v}_{0}}{\|\bm{M}^{t}\bm{v}_{0}\|_{2}}.

(2.1.19)

Because all the eigenvectors $\bm{w}_{i}$ of $\bm{M}$ form an orthonormal basis for $\mathbb{R}^{D}$ , we can write

\bm{v}_{0}=\sum_{i=1}^{D}\alpha_{i}\bm{w}_{i},

(2.1.20)

where because $\bm{v}_{0}$ is Gaussian the $\alpha_{i}$ are all nonzero with probability $1$ . Thus, we can use our earlier expression for $\bm{v}_{t}$ to write

\bm{v}_{t}=\frac{\bm{M}^{t}\bm{v}_{0}}{\|\bm{M}^{t}\bm{v}_{0}\|_{2}}=\frac{\sum_{i=1}^{D}\lambda_{i}^{t}\alpha_{i}\bm{w}_{i}}{\|\sum_{i=1}^{D}\lambda_{i}^{t}\alpha_{i}\bm{w}_{i}\|_{2}}=\frac{\sum_{i=1}^{D}\lambda_{i}^{t}\alpha_{i}\bm{w}_{i}}{\sum_{i=1}^{D}\lambda_{i}^{t}\lvert\alpha_{i}\rvert}.

(2.1.21)

Now, let us consider the case where $\lambda_{1}>\lambda_{2}\geq\cdots\geq\lambda_{D}\geq 0$ . (The case with repeated top eigenvalues goes similarly.) Then we can write

\bm{v}_{t}=\frac{\alpha_{1}\bm{w}_{1}+\sum_{i=2}^{D}(\lambda_{i}/\lambda_{1})^{t}\alpha_{i}\bm{w}_{i}}{\lvert\alpha_{1}\rvert+\sum_{i=2}^{D}(\lambda_{i}/\lambda_{1})^{t}\lvert\alpha_{i}\rvert}.

(2.1.22)

Because $\lambda_{1}>\lambda_{i}$ for all $i>1$ , the terms inside the summation go to $0$ exponentially fast, and the remainder is the limit

\lim_{t\to\infty}\bm{v}_{t}=\frac{\alpha_{1}}{\lvert\alpha_{1}\rvert}\bm{w}_{1}=\operatorname{sign}(\alpha_{1})\bm{w}_{1},

(2.1.23)

which is a top unit eigenvector of $\bm{M}$ . The top eigenvalue $\lambda_{1}$ of $\bm{M}$ can be estimated by $\bm{v}_{t}^{\top}\bm{M}\bm{v}_{t}$ , which converges similarly rapidly to $\lambda_{1}$ . ∎

To find the second top eigenvector, we apply the power iteration algorithm to $\bm{M}-\lambda_{1}\bm{v}_{1}\bm{v}_{1}^{\top}$ , which has eigenvectors $(\bm{w}_{i})_{i=2}^{D}$ and corresponding eigenvalues $(\lambda_{i})_{i=2}^{D}$ . By repeating this procedure $d$ times in sequence, we can very efficiently estimate the top $d$ eigenvectors of $\bm{M}$ for any symmetric positive semidefinite matrix $\bm{M}$ . Thus we can also apply it to $\bm{X}\bm{X}^{\top}/N$ to recover the top $d$ principal components, which is what we were after in the first place. Notice that this approach recovers one principal component at a time; we will contrast this to other algorithmic approaches, such as gradient descent on global objective functions, in future sections.

2.1.3 Probabilistic PCA

Notice that the above formulation makes no statistical assumptions on the data-generating process. However, it is common to include statistical elements within a given data model, as it may add further enlightening interpretations about the result of the analysis. As such, we ask the natural question: what is the statistical analogue to low-dimensional structure? Our answer is that a low-dimensional distribution is one whose support is concentrated around a low-dimensional geometric structure.

To illustrate this point, we discuss probabilistic principal component analysis (PPCA). This formulation can be viewed as a statistical variant of regular PCA. Mathematically, we now consider our data as samples from a random variable $\bm{x}$ taking values in $\mathbb{R}^{D}$ (also sometimes called a random vector). We say that $\bm{x}$ has (approximate) low-rank statistical structure if and only if there exists an orthonormal matrix $\bm{U}\in\mathsf{O}(D,d)$ , a random variable $\bm{z}$ taking values in $\mathbb{R}^{d}$ , and a small random variable $\bm{\varepsilon}$ taking values in $\mathbb{R}^{D}$ such that $\bm{z}$ and $\bm{\varepsilon}$ are independent, and

\bm{x}=\bm{U}\bm{z}+\bm{\varepsilon}.

(2.1.24)

Our goal is again to recover $\bm{U}$ . Towards this end, we set up the analogous problem as in Section 2.1.1, i.e., optimizing over subspace supports $\tilde{\bm{U}}$ and random variables $\bm{z}$ to solve the following problem:

\min_{\tilde{\bm{U}},\tilde{\bm{z}}}\operatorname{\mathbb{E}}\|\bm{x}-\tilde{\bm{U}}\tilde{\bm{z}}\|_{2}^{2}.

(2.1.25)

Since we are finding the best such random variable $\bm{z}$ we can find its realization $\bm{z}(\bm{x})$ separately for each value of $\bm{x}$ . Performing the same calculations as in Section 2.1.1, we obtain

\min_{\tilde{\bm{U}},\tilde{\bm{z}}}\operatorname{\mathbb{E}}\|\bm{x}-\tilde{\bm{U}}\tilde{\bm{z}}\|_{2}^{2}=\min_{\tilde{\bm{U}}}\operatorname{\mathbb{E}}\min_{\tilde{\bm{z}}(\bm{x})}\|\bm{x}-\tilde{\bm{U}}\tilde{\bm{z}}(\bm{x})\|_{2}^{2}=\min_{\tilde{\bm{U}}}\operatorname{\mathbb{E}}\|\bm{x}-\tilde{\bm{U}}\tilde{\bm{U}}^{\top}\bm{x}\|_{2}^{2},

(2.1.26)

again re-emphasizing the fact that the estimated subspace with principal components $\bm{U}^{\star}$ corresponds to a denoiser $\bm{U}^{\star}(\bm{U}^{\star})^{\top}$ which projects onto that subspace. As before, we obtain

	$\displaystyle\operatorname*{arg\ min}_{\tilde{\bm{U}}}\operatorname{\mathbb{E}}\\|\bm{x}-\tilde{\bm{U}}\tilde{\bm{U}}^{\top}\bm{x}\\|_{2}^{2}$	$\displaystyle=\operatorname*{arg\ max}_{\tilde{\bm{U}}}\operatorname{\mathbb{E}}\\|\tilde{\bm{U}}^{\top}\bm{x}\\|_{2}^{2}$		(2.1.27)
		$\displaystyle=\operatorname*{arg\ max}_{\tilde{\bm{U}}}\operatorname{tr}(\tilde{\bm{U}}^{\top}\operatorname{\mathbb{E}}[\bm{x}\bm{x}^{\top}]\tilde{\bm{U}}),$		(2.1.28)

and the solution to the latter problem is just the top $d$ principal components of the second moment matrix $\operatorname{\mathbb{E}}[\bm{x}\bm{x}^{\top}]$ . Actually, the above problems are visually very similar to the equations for computing the principal components in the previous subsection, except with $\operatorname{\mathbb{E}}[\bm{x}\bm{x}^{\top}]$ replacing $\bm{X}\bm{X}^{\top}/N$ . In fact, the latter quantity is an estimate for the former. Both formulations effectively do the same thing, and have the same practical solution—compute the left singular vectors of the data matrix $\bm{X}$ , or equivalently the top eigenvectors of the estimated covariance matrix $\bm{X}\bm{X}^{\top}/N$ . The statistical formulation, however, has an additional interpretation. Suppose that $\operatorname{\mathbb{E}}[\bm{z}]=\bm{0}$ and $\operatorname{\mathbb{E}}[\bm{\varepsilon}]=\bm{0}$ . We have

\operatorname{\mathbb{E}}[\bm{x}]=\bm{U}\operatorname{\mathbb{E}}[\bm{z}]+\operatorname{\mathbb{E}}[\bm{\varepsilon}]=\bm{0},

(2.1.29)

so that $\operatorname{Cov}[\bm{x}]=\operatorname{\mathbb{E}}[\bm{x}\bm{x}^{\top}]$ . Now working out $\operatorname{Cov}[\bm{x}]$ we have

\operatorname{Cov}[\bm{x}]=\bm{U}\operatorname{Cov}[\bm{z}]\bm{U}^{\top}+\operatorname{Cov}[\bm{\varepsilon}]=\bm{U}\operatorname{\mathbb{E}}[\bm{z}\bm{z}^{\top}]\bm{U}^{\top}+\operatorname{Cov}[\bm{\varepsilon}].

(2.1.30)

In particular, if $\operatorname{Cov}[\bm{\varepsilon}]$ is small, it holds that $\operatorname{Cov}[\bm{x}]=\operatorname{\mathbb{E}}[\bm{x}\bm{x}^{\top}]$ is approximately a low-rank matrix, in particular rank $d$ . Thus the top $d$ eigenvectors of $\operatorname{\mathbb{E}}[\bm{x}\bm{x}^{\top}]$ essentially summarize the whole covariance matrix. But they are also the principal components, so we can interpret principal component analysis as performing a low-rank decomposition of $\operatorname{Cov}[\bm{x}]$ .

Remark 2.5.

By using the probabilistic viewpoint of PCA, we achieve a clearer and more quantitative understanding of how it relates to denoising. First, consider the denoising problem in (2.1.26), namely

\min_{\tilde{\bm{U}}}\operatorname{\mathbb{E}}\|\bm{x}-\tilde{\bm{U}}\tilde{\bm{U}}^{\top}\bm{x}\|_{2}^{2}.

(2.1.31)

It is not too hard to prove that if $\tilde{\bm{U}}$ has $d$ columns and if $\bm{\varepsilon}$ is an isotropic Gaussian random variable, i.e., with distribution $\bm{\varepsilon}\sim\operatorname{\mathcal{N}}(\bm{0},\sigma^{2}\bm{I})$ ,³³3Other distributions work so long as they support all of $\mathbb{R}^{D}$ , but the Gaussian is the easiest to work with here. then for any optimal solution $\bm{U}^{\star}$ to this problem, we have

\bm{U}^{\star}(\bm{U}^{\star})^{\top}=\bm{U}\bm{U}^{\top}

(2.1.32)

and so the true supporting subspace, say $\mathcal{S}\doteq\mathop{\mathrm{col}}(\bm{U})$ , is recovered as the span of columns of $\bm{U}^{\star}$ , since

\mathcal{S}=\mathop{\mathrm{col}}(\bm{U})=\mathop{\mathrm{col}}(\bm{U}\bm{U}^{\top})=\mathop{\mathrm{col}}(\bm{U}^{\star}(\bm{U}^{\star})^{\top})=\mathop{\mathrm{col}}(\bm{U}^{\star}).

(2.1.33)

In particular, the learned denoising map $\bm{U}^{\star}(\bm{U}^{\star})^{\top}$ is an orthogonal projection onto $\mathcal{S}$ , pushing noisy points onto the ground truth supporting subspace. We can establish a similar technical result in the case where we only have finite samples, as in Theorem 2.1, but this takes more effort and technicality. Summarizing this discussion, we have the following informal Theorem.

Theorem 2.3.

Suppose that the random variable $\bm{x}$ can be written as

\bm{x}=\bm{U}\bm{z}+\bm{\varepsilon}

(2.1.34)

where $\bm{U}\in\mathsf{O}(D,d)$ captures the low-rank structure, $\bm{z}$ is a random vector taking values in $\mathbb{R}^{d}$ , and $\bm{\varepsilon}$ is a random vector taking values in $\mathbb{R}^{D}$ such that $\bm{z}$ and $\bm{\varepsilon}$ are independent, and $\bm{\varepsilon}$ is small. Then the principal components $\bm{U}^{\star}\in\mathsf{O}(D,d)$ of our dataset are given by the top $d$ eigenvectors of $\operatorname{\mathbb{E}}[\bm{x}\bm{x}^{\top}]$ , and approximately correspond to the optimal linear denoiser: $\bm{U}^{\star}(\bm{U}^{\star})^{\top}\approx\bm{U}\bm{U}^{\top}$ .

2.1.4 Matrix Completion

In the previous subsections, we discussed the problem of learning a low-rank geometric or statistical distribution, where the data were sampled from a subspace with additive noise. But this is not the only type of disturbance from a low-dimensional distribution that is worthwhile to study. In this subsection, we introduce one more class of non-additive errors which become increasingly important in deep learning. Let us consider the case where we have some data $\{\bm{x}_{i}\}_{i=1}^{N}$ generated according to (2.1.1). Now we arrange them into a matrix $\bm{X}=\begin{bmatrix}\bm{x}_{1},\dots,\bm{x}_{N}\end{bmatrix}\in\mathbb{R}^{D\times N}$ . Unlike before, we do not observe $\bm{X}$ directly; we instead imagine that our observation was corrupted en route and we obtained

\bm{Y}=\bm{M}\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}\bm{X},

(2.1.35)

where $\bm{M}\in\{0,1\}^{D\times N}$ is a mask that is known by us, and $\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}$ is element-wise multiplication. In this case, our goal is to recover $\bm{X}$ (from which point we can use PCA to recover $\bm{U}$ , etc), given only the corrupted observation $\bm{Y}$ , the mask $\bm{M}$ , and the knowledge that $\bm{X}$ is (approximately) low-rank. This is called low-rank matrix completion.

There are many excellent resources discussing algorithms and approaches to solve this problem [WM22]. Indeed, this and similar generalizations of this low-rank structure recovery problem are resolved by “robust PCA”. We will not go into the solution method here. Instead, we will discuss under what conditions this problem is plausible to solve. On one hand, in the most absurd case, suppose that each entry of the matrix $\bm{X}$ were chosen independently from all the others. Then there would be no hope of recovering $\bm{X}$ exactly even if only one entry were missing and we had $DN-1$ entries. On the other hand, suppose that we knew that $\bm{X}$ has rank 1 exactly, which is an extremely strong condition on the low-dimensional structure of the data, and we were dealing with the mask

\bm{M}=\begin{bmatrix}\bm{1}_{(D-1)\times 1}&\bm{0}_{(D-1)\times(N-1)}\\ 1&\bm{1}_{1\times(N-1)}\end{bmatrix}.

(2.1.36)

Then we know that the data are distributed on a line, and we know a vector on that line—it is just the first column of the matrix $\bm{Y}=\bm{M}\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}\bm{X}$ . From the last coordinate of each column, also revealed to us by the mask, we can solve for each column since for each final coordinate there is only one vector on the line with this coordinate. Thus we can reconstruct $\bm{X}$ with perfect accuracy, and we only need a linear number of observations $D+N-1$ .

In the real world, actual problems are somewhere in between the two limiting cases discussed above. Yet the differences between these two extremes, as well as the earlier discussion of PCA, reveal a general kernel of truth:

The lower-dimensional and more structured the data distribution is, the easier it is to process, and the fewer observations are needed—provided that the algorithm effectively utilizes this low-dimensional structure.

As is perhaps predictable, we will encounter this motif repeatedly in the remainder of the manuscript, starting in the very next section.

2.2 A Mixture of Complete Low-Dimensional Subspaces

As we have seen, low-rank signal models are rich enough to provide a full picture of the interplay between low-dimensionality in data and efficient and scalable computational algorithms for representation and recovery under errors. These models imply a linear and symmetric representation learning pipeline (2.1.6):

\bm{z}=\mathcal{E}(\bm{x})=\bm{U}^{\top}\bm{x},\quad\hat{\bm{x}}=\mathcal{D}(\bm{z})=\bm{U}\bm{z},

which can be provably learned from finite samples of $\bm{x}$ with principal component analysis (solved efficiently, say, with the power method) whenever the distribution of $\bm{x}$ truly is linear. This is a restrictive assumption—for as Harold Hotelling, the distinguished 20th century statistician,⁴⁴4By coincidence, also famous for his contributions to the development and naming of Principal Component Analysis [Hot33]. objected following George Dantzig’s presentation of his theory of linear programming for the first time [Dan02],

“…we all know the world is nonlinear.”

Figure 2.3 : Left: features tracked on independently moving objects in a scene. Right: image patches associated with different regions of an image. — Figure 2.3: Left: features tracked on independently moving objects in a scene. Right: image patches associated with different regions of an image.

Even accounting for its elegance and simplicity, the low-rank assumption is too restrictive to be broadly applicable to modeling real-world data. A key limitation is the assumption of a single linear subspace that is responsible for generating the structured observations. In many practical applications, structure generated by a mixture of distinct low-dimensional subspaces provides a more realistic model. For example, consider a video sequence that captures the motion of several distinct objects, each subject to its own independent displacement (Figure 2.3 left). Under suitable assumptions on the individual motions, each object becomes responsible for an independent low-dimensional subspace in the concatenated sequence of video frames [VM04]. As another example, consider modeling natural images via learning a model for the distribution of patches, spatially-contiguous collections of pixels, within an image (Figure 2.3 right). Unlike in the Eigenface example we saw previously, where images of faces with matched poses can be well-approximated by a single low-dimensional subspace, the patch at a specific location in a natural image can correspond to objects with very different properties—for example, distinct color or shape due to occlusion boundaries. Therefore, modeling the distribution of patches with a single subspace is futile, but a mixture of subspaces, one for each region, performs surprisingly well in practice, say for segmenting or compressing purposes [MRY+11].⁵⁵5We will return to this observation in Chapter 4, where we will show it can be significantly generalized to learn representations for large-scale modern datasets. We will see a concrete example in the next chapter (Example 3.12).

Figure 2.4 : Data on a mixture of low-dimensional subspaces, say 𝒮 j = col ( 𝑼 j ) \mathcal{S}_{j}=\mathop{\mathrm{col}}(\bm{U}_{j}) caligraphic_S start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT = roman_col ( bold_italic_U start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ) . — Figure 2.4: Data on a mixture of low-dimensional subspaces, say $\mathcal{S}_{j}=\mathop{\mathrm{col}}(\bm{U}_{j})$ .

In this section, we will begin by discussing the conceptual and algorithmic foundations for structured representation learning when the data distribution can be modeled by a mixture of low-dimensional subspaces, as illustrated in Figure 2.4. In this setting, the decoder mapping will be almost as simple as the case of a single subspace: we simply represent via

\hat{\bm{x}}=\mathcal{D}(\bm{z})=\left(\sum_{k=1}^{K}\pi_{k}(\bm{z})\bm{U}_{k}\right)\bm{z},

(2.2.1)

where $\pi_{k}:\mathbb{R}^{d}\to\{0,1\}$ are a set of sparse weighting random variables, such that only a single subspace $\mathcal{S}_{k}=\mathop{\mathrm{col}}(\bm{U}_{k})$ is selected in the sum. However, the task of encoding such data $\bm{x}$ to suitable representations $\bm{z}$ , and learning such an encoder-decoder pair from data, will prove to be more involved.

We will see how ideas from the rich literature on sparse representation and independent component analysis lead to a natural reformulation of the above decoder architecture through the lens of sparsity, the corresponding encoder architecture (obtained through a power-method-like algorithm analogous to that of principal component analysis), and strong guarantees of correctness and efficiency for learning such encoder-decoder pairs from data. In this sense, the case of mixed linear low-dimensional structure already leads to many of the key conceptual components of structured representation learning that we will develop in far greater generality in this book.

2.2.1 Mixtures of Subspaces and Sparsely-Used Dictionaries

Let $\bm{U}_{1},\dots,\bm{U}_{K}$ , each of size $D\times d$ , denote a collection of orthonormal bases for $K$ subspaces of dimension $d$ in $\mathbb{R}^{D}$ . To say that $\bm{x}$ follows a mixture-of-subspaces distribution parameterized by $\bm{U}_{1},\dots,\bm{U}_{K}$ means, geometrically speaking, that

\bm{x}=\bm{U}_{k}\bm{z}\quad\text{for some}\enspace k\in[K],\enspace\bm{z}\in\mathbb{R}^{d}.

(2.2.2)

The statistical analogue of this geometric model, as we saw for the case of PCA and linear structure, is that $\bm{x}$ follows a mixture of Gaussians distribution: that is,

\bm{x}\sim\sum_{k=1}^{K}\pi_{k}\mathcal{N}(\mathbf{0},\bm{U}_{k}\bm{U}_{k}^{\top}),\quad\text{for some}\enspace\pi_{k}\geq 0,\enspace\sum_{k=1}^{K}\pi_{k}=1.

(2.2.3)

In other words, for each $k\in[K]$ , $\bm{x}$ is Gaussian on the low-dimensional subspace $\mathop{\mathrm{col}}(\bm{U}_{k})$ with probability $\pi_{k}$ .

Remark 2.6 (A Mixture of Gaussians v.s. A Superposition of Gaussians).

One should be aware that the above model (2.2.3) is a mixture of Gaussian distributions, not to be confused with a mixture of Gaussian variables by superposition, say

\bm{x}=\sum_{i=1}^{n}w_{i}\bm{x}_{i},\quad\bm{x}_{i}\sim\mathcal{N}(\mathbf{0},\bm{U}_{i}\bm{U}_{i}^{\top}),

(2.2.4)

where $\bm{x}_{i}$ are independent random Gaussian vectors and $w_{i}$ are a set of fixed weights. As we know from the properties of Gaussian vectors, such a superposition $\bm{x}$ will remain a Gaussian distribution.

For now, we focus on the geometric perspective offered by (2.2.2). There is an algebraically convenient alternative to this conditional representation. Consider a lifted representation vector $\bm{z}=[\bm{z}_{1}^{\top},\dots,\bm{z}_{K}^{\top}]^{\top}\in\mathbb{R}^{dK}$ , such that $\bm{z}$ is $d$ -sparse with support on one of the $K$ consecutive non-overlapping blocks of $d$ coordinates out of $dK$ . Then (2.2.2) can be written equivalently as

\bm{x}=\underbrace{\begin{bmatrix}|&\ldots&|\\ \bm{U}_{1}&\ldots&\bm{U}_{K}\\ |&\ldots&|\end{bmatrix}}_{\bm{U}}\underbrace{\begin{bmatrix}\bm{z}_{1}\\ \vdots\\ \bm{z}_{K}\end{bmatrix}}_{\bm{z}},\quad\left\|\begin{bmatrix}\left\|\bm{z}_{1}\right\|_{2}\\ \vdots\\ \left\|\bm{z}_{K}\right\|_{2}\end{bmatrix}\right\|_{0}=1.

(2.2.5)

Here, the $\ell^{0}$ “norm” $\|\,\cdot\,\|_{0}$ measures sparsity by counting the number of nonzero entries:

\|\bm{z}\|_{0}=\left\lvert\{i\mid z_{i}\neq 0\}\right\rvert,

(2.2.6)

and the matrix $\bm{U}\in\mathbb{R}^{D\times Kd}$ is called a dictionary with all the $\{\bm{U}_{i}\}_{i=1}^{K}$ as code words. In general, if the number of subspaces in the mixture $K$ is large enough, there is no bound on the number of columns contained in the dictionary $\bm{U}$ . In the case where $Kd<D$ , $\bm{U}$ is called undercomplete; when $Kd=D$ , it is called complete; and when $Kd>D$ , it is called overcomplete.

Now, (2.2.5) suggests a convenient relaxation for tractability of analysis: rather than modeling $\bm{x}$ as coming from a mixture of $K$ specific subspaces $\bm{U}_{1},\dots,\bm{U}_{K}$ , we may instead start with a dictionary $\bm{U}\in\mathbb{R}^{D\times m}$ , where $m$ may be smaller or larger than $D$ , and simply seek to represent $\bm{x}=\bm{U}\bm{z}$ with the sparsity $\|\bm{z}\|_{0}$ sufficiently small. This leads to the sparse dictionary model for $\bm{x}$ :

\bm{x}=\bm{U}\bm{z}+\bm{\varepsilon},\quad\|\bm{z}\|_{0}\ll d,

(2.2.7)

where $\bm{\varepsilon}$ represents an unknown noise vector. Geometrically, this implies that $\bm{x}$ lies close to the span of a subset of $\|\bm{z}\|_{0}$ columns of $\bm{U}$ , making this an instantiation of the mixture-of-subspaces model (2.2.2) with a very large value of $K$ , and specific correlations between the subspaces $\bm{U}_{k}$ .

Orthogonal dictionary for sparse coding.

Now we can formulate the structured representation learning problem for mixtures of low-dimensional subspaces that we will study in this section. We assume that we have samples $\bm{X}=[\bm{x}_{1},\dots,\bm{x}_{N}]$ from an unknown sparse dictionary model (2.2.7), possibly with added noises $\bm{\varepsilon}_{i}$ . Let us begin from the assumption that the dictionary $\bm{U}$ in the sparse dictionary model (2.2.7) is complete and orthogonal,⁶⁶6It can be shown that for the complete case, we do not lose any generality by making the orthogonal assumption (Exercise 2.4). and that each coefficient vector $\bm{z}$ is $d$ -sparse, with $d\ll D$ . In this setting, representation learning amounts to correctly learning the orthogonal dictionary $\bm{U}$ via optimization: we can then take $\mathcal{E}(\bm{x})=\bm{U}^{\top}\bm{x}$ as the encoder and $\mathcal{D}(\bm{z})=\bm{U}\bm{z}$ as the decoder, and $\mathcal{D}=\mathcal{E}^{-1}$ . In diagram form:

\bm{x}\xrightarrow{\hskip 5.69054pt\mathcal{E}=\bm{U}^{\top}\hskip 5.69054pt}\bm{z}\xrightarrow{\hskip 5.69054pt\mathcal{D}=\bm{U}\hskip 5.69054pt}\hat{\bm{x}}.

(2.2.8)

We see that the autoencoding pair $(\mathcal{E},\mathcal{D})$ for complete dictionary learning is symmetric, as in the case of a single linear subspace, making the computational task of encoding and decoding no more difficult than in the linear case. On the other hand, the task of learning the dictionary $\bm{U}$ is strictly more difficult than learning a single linear subspace by PCA. To see why we cannot simply use PCA to learn the orthogonal dictionary $\bm{U}$ correctly, note that the loss function that gave rise to PCA, namely (2.1.5), is completely invariant to rotations of the rows of the matrix $\bm{U}$ : that is, if $\bm{Q}$ is any $d\times d$ orthogonal matrix, then $\bm{U}$ and $\bm{U}\bm{Q}$ are both feasible and have an identical loss for (2.1.5). The sparse dictionary model is decidedly not invariant to such transformations: if we replaced $\bm{U}$ by $\bm{U}\bm{Q}$ and made a corresponding rotation $\bm{Q}^{\top}\bm{z}$ of the representation coefficients $\bm{z}$ , we would in general destroy the sparsity structure of $\bm{z}$ , violating the modeling assumption. Thus, we need new algorithms for learning orthogonal dictionaries.

2.2.2 Complete Dictionary Learning

In this section, we will derive algorithms for solving the orthogonal dictionary learning problem. To be more precise, we assume that the observed vector $\bm{x}\in\mathbb{R}^{D}$ follows a statistical model

\bm{x}=\bm{U}\bm{z}+\bm{\varepsilon},

(2.2.9)

where $\bm{U}\in\mathbb{R}^{D\times D}$ is an unknown orthogonal dictionary, $\bm{z}$ is a random vector with statistically independent components $z_{i}$ , each with zero mean, and $\bm{\varepsilon}\in\mathbb{R}^{D}$ is an independent random vector of small (Gaussian) noises. The goal is to recover $\bm{U}$ (and hence $\bm{z}$ ) from samples of $\bm{x}$ .

Here we assume that each independent component $z_{i}$ is distributed as

z_{i}\sim\mathrm{Bern}(\theta)\cdot\mathcal{N}(0,1/\theta).

That is, it is the product of a Bernoulli random variable with probability $\theta$ of being $1$ and $1-\theta$ of being $0$ , and an independent Gaussian random variable with variance $1/\theta$ . This distribution is formally known as the Bernoulli-Gaussian distribution. The normalization is chosen so that $\operatorname{Var}(z_{i})=1$ and hence $\mathbb{E}[\|\bm{z}\|_{2}^{2}]=d$ . This modeling assumption implies that the vector of independent components $\bm{z}$ is typically very sparse: we calculate $\mathbb{E}\left[\|\bm{z}\|_{0}\right]=d\theta$ , which is small when $\theta$ is inversely proportional to $d$ .

Remark 2.7 (The Orthogonal Assumption).

At first sight, the assumption that the dictionary $\bm{U}$ is orthogonal might seem to be somewhat restrictive. But there is actually no loss of generality. One may consider a complete dictionary to be any square invertible matrix $\bm{U}$ . With samples generated from this dictionary: $\bm{X}=\bm{U}\bm{Z}\in\mathbb{R}^{D\times N}$ , it is easy to show that with some preconditioning of the data matrix $\bm{X}$ :

\bar{\bm{X}}=\Big{(}\frac{1}{N\theta}\bm{X}\bm{X}^{\top}\Big{)}^{-\frac{1}{2}}\bm{X},

(2.2.10)

then there exists an orthogonal matrix $\bm{U}_{o}\in\mathsf{O}(D)$ such that

\bar{\bm{X}}=\bm{U}_{o}\bm{Z}.

(2.2.11)

See Exercise 2.4 or [SQW17a] for more details.

Figure 2.5 : Maximizing ℓ 4 \ell^{4} roman_ℓ start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT norm or minimizing ℓ 1 \ell^{1} roman_ℓ start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT norm promotes sparsity (for vectors on the sphere). — Figure 2.5: Maximizing $\ell^{4}$ norm or minimizing $\ell^{1}$ norm promotes sparsity (for vectors on the sphere).

Dictionary learning via the MSP algorithm.

Now suppose that we are given a set of observations:

\bm{x}_{i}=\bm{U}\bm{z}_{i}+\bm{\varepsilon}_{i},\ \forall i\in[N].

(2.2.12)

Let $\bm{X}=[\bm{x}_{1},\dots,\bm{x}_{N}]$ and $\bm{Z}=[\bm{z}_{1},\dots,\bm{z}_{N}]$ . The goal is to recover $\bm{U}$ from the data $\bm{X}$ . Therefore, given any orthogonal matrix $\bm{A}\in\mathsf{O}(D)$ ,

\bm{A}\bm{x}_{i}=\bm{A}\bm{U}\bm{z}_{i}+\bm{A}\bm{\varepsilon}_{i}

(2.2.13)

would be nearly sparse if $\bm{A}=\bm{U}^{T}$ (as by assumption the noise $\bm{\varepsilon}_{i}$ is of small magnitude).

Also, given $\bm{U}$ is orthogonal and the fact $\bm{\varepsilon}$ is small, the vector $\bm{x}$ has a predictable expected norm, i.e., $\mathbb{E}[\|\bm{x}\|_{2}^{2}]\approx\mathbb{E}[\|\bm{z}\|_{2}^{2}]=d$ . It is a known fact that for vectors on a sphere, maximizing the $\ell^{4}$ norm is equivalent to minimizing the $\ell^{0}$ norm (for promoting sparsity),

\operatorname*{arg\ max}_{\bm{z}\in\mathbb{S}^{n}}\|\bm{z}\|_{4}\quad\Leftrightarrow\quad\operatorname*{arg\ min}_{\bm{z}\in\mathbb{S}^{n}}\|\bm{z}\|_{0}.

(2.2.14)

This is illustrated in Figure 2.5.

An orthogonal matrix $\bm{A}$ preserves the Euclidean ( $\ell^{2}$ ) norm: $\|\bm{A}\bm{x}\|_{2}^{2}=\|\bm{x}\|_{2}^{2}$ . Therefore, to find the correct orthogonal dictionary $\bm{U}$ from $\bm{X}$ , we may try to solve the following optimization program:

\max_{\tilde{\bm{A}}\in\mathsf{O}(D)}\,\frac{1}{4}\left\|\tilde{\bm{A}}\bm{X}\right\|_{4}^{4}=\frac{1}{4}\sum_{i=1}^{N}\left\|\tilde{\bm{A}}\bm{x}_{i}\right\|_{4}^{4}

(2.2.15)

This is known as the $\ell^{4}$ maximization problem [ZMZ+20]. After we find the solution $\bm{A}^{\star}$ , we can take the transpose $\bm{U}^{\star}=(\bm{A}^{\star})^{\top}$ .

Remark 2.8.

It is also known that for vectors on a sphere, minimizing the $\ell^{1}$ norm is equivalent to minimizing the $\ell^{0}$ norm (for promoting sparsity),

\operatorname*{arg\ min}_{\bm{z}\in\mathbb{S}^{n}}\|\bm{z}\|_{1}\quad\Leftrightarrow\quad\operatorname*{arg\ min}_{\bm{z}\in\mathbb{S}^{n}}\|\bm{z}\|_{0},

which is also illustrated in Figure 2.5. This fact can also be exploited to learn the dictionary $\bm{A}$ effectively and efficiently. This was actually explored earlier than the $\ell^{4}$ norm used here. Interested readers may refer to the work of [QZL+20a].

Note that the above problem is equivalent to the following constrained optimization problem:

\min\,-\frac{1}{4}\left\|\tilde{\bm{A}}\bm{X}\right\|_{4}^{4}\quad\mbox{subject to}\quad\tilde{\bm{A}}^{\top}\tilde{\bm{A}}=\bm{I}.

(2.2.16)

As shown in [WM22], using the Lagrange multiplier method, one can derive that the optimal solution to the problem should satisfy the following “fixed point” condition:

\bm{A}^{\star}=\mathcal{P}_{\mathrm{O}(D)}[({\bm{A}^{\star}\bm{X}})^{\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}3}\bm{X}^{\top}],

(2.2.17)

where $\mathcal{P}_{\mathrm{O}(D)}[\,\cdot\,]$ is a projection onto the space of orthogonal matrices $\mathrm{O}(D)$ .⁷⁷7For any matrix $\bm{M}$ with SVD $\bm{M}=\bm{U}\bm{\Sigma}\bm{V}^{\top}$ , $\mathcal{P}_{\mathrm{O}(D)}[\bm{M}]=\bm{U}\bm{V}^{\top}$ . We leave this as an exercise for the reader.

To compute the fixed point for the above equation, similar to how we computed eigenvectors for PCA (2.1.16), we may take the following power iteration:

\bm{A}_{t+1}=\mathcal{P}_{\mathrm{O}(D)}[({\bm{A}_{t}\bm{X}})^{\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}3}\bm{X}^{\top}].

(2.2.18)

This is known as the matching, stretching, and projection (MSP) algorithm proposed by [ZMZ+20]. It was shown that under broad conditions such a greedy algorithm indeed converges to the correct solution at a superlinear rate.

Remark 2.9 (Global Optimality of $\ell^{4}$ Maximization).

The constrained $\ell^{4}$ maximization problem is a nonconvex program. In general one should not expect that any greedy (say gradient-descent type) algorithms would converge to the globally optimal solution. Surprisingly, one can show that, unlike general nonconvex programs, the landscape of $\ell^{4}$ maximization over a sphere

\min\,-\frac{1}{4}\left\|\bm{q}^{\top}\bm{X}\right\|_{4}^{4}\quad\mbox{subject to}\quad\bm{q}^{\top}\bm{q}=1.

(2.2.19)

is very benign: All local minima are close to the global optima and all critical points are saddle points with a direction of negative curvature. Hence, any descent method with the ability of escaping strict saddle points provably finds global optimal solutions. For more precise statements, interested readers may refer to [QZL+20].

Remark 2.10 (Stable Deep Linear Network).

The above iterative process of computing the dictionary has a natural incremental “deep learning” interpretation. Let us define $\delta\bm{A}_{t+1}=\bm{A}_{t+1}\bm{A}_{t}^{\top}$ and $\bm{Z}_{t}=\bm{A}_{t}\bm{X}$ , then it is easy to show that

\delta\bm{A}_{t+1}=\mathcal{P}_{\mathrm{O}(D)}[(\bm{Z}_{t})^{\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}3}\bm{Z}_{t}^{\top}].

If $\bm{A}_{t}$ converges to the correct dictionary $\bm{D}_{o}$ , then the above iterative encoding process is essentially equivalent to a “deep linear network”:

\bm{Z}\;\longleftarrow\;\bm{Z}_{t+1}=\underbrace{\delta\bm{A}_{t+1}\delta\bm{A}_{t}\ldots\delta\bm{A}_{1}}_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\text{forward constructed layers}}\bm{X}.

Note that the computation of the increment transforms $\delta\bm{A}_{t+1}$ at each layer depends only on the feature output from the previous layer $\bm{Z}_{t}$ . The network is naturally stable as each layer is a norm-preserving orthogonal transform. Despite its resemblance to a linear deep network, backpropagation is unnecessary to learn each layer. All layers are learned in one forward pass!

2.2.3 Connection to ICA and Kurtosis

With the Bernoulli-Gaussian model, the variables $z_{i}$ are independent and non-Gaussian. Then, there is a clear correspondence between the dictionary learning and the classic independent component analysis (ICA), to the extent that algorithms to solve one problem can be used to solve the other.⁸⁸8We explore this issue in more depth in Exercise 2.3, where a connection between non-Gaussianity of the independent components and the purely geometric notion of symmetry is made. This issue is related to our observation above that PCA does not work for recovering sparsely-used orthogonal dictionaries: in the statistical setting, it can be related to rotational invariance of the Gaussian distribution (Exercise 2.2).

Towards deriving an algorithm based on ICA, we focus on an objective function known as kurtosis, which is used in ICA as a direct consequence of the non-Gaussianity of the components. The kurtosis, or fourth-order cumulant, of a zero-mean random variable $X$ is defined as

\mathop{\mathrm{kurt}}(X)=\operatorname{\mathbb{E}}{X^{4}}-3(\operatorname{\mathbb{E}}{X^{2}})^{2}.

(2.2.20)

If we have only finite samples from the random variable $X$ arranged into a vector $\bm{x}=[x_{1},\dots,x_{N}]$ , we define kurtosis through their empirical average, which yields

\mathop{\mathrm{kurt}}(\bm{x})=\frac{1}{N}\|\bm{x}\|_{4}^{4}-\frac{3}{N^{2}}\|\bm{x}\|_{2}^{4}.

(2.2.21)

Finally, for random vectors, we define their kurtosis as the sum of each component’s scalar kurtosis. Kurtosis is a natural loss function for ICA because for Gaussian $X$ , kurtosis is zero; the reader can verify further that the Bernoulli-Gaussian distribution has positive kurtosis. Thus a natural procedure for seeking non-Gaussian independent components is to search for a set of mutually-orthogonal directions $\bm{V}\in\mathbb{R}^{d\times k}$ such that $\bm{V}^{\top}\bm{X}$ has maximal kurtosis, where $\bm{X}=\bm{U}\bm{Z}\in\mathbb{R}^{D\times N}$ is the Bernoulli-Gaussian ICA data matrix. Formally, we seek to solve the problem

\max_{\bm{V}^{\top}\bm{V}=\bm{I}}\mathop{\mathrm{kurt}}(\bm{V}^{\top}\bm{X}).

(2.2.22)

At one extreme, we may set $k=D$ and seek to recover the entire dictionary $\bm{U}$ in a single shot. It can be shown that this problem can be solved with the MSP algorithm we have seen previously. At the other extreme, we may set $k=1$ and seek to recover a single direction (column of $\bm{U}$ ) at a time, performing deflation, i.e., replacing the data matrix $\bm{X}$ by $(\bm{I}-\bm{u}\bm{u}^{\top})\bm{X}$ , after each step before finding another direction. There is a natural tradeoff between the scalability of the $k=1$ incremental approach and the efficiency and robustness of the $k=D$ approach.

Incremental ICA: correctness and FastICA algorithm.

The FastICA algorithm, advanced by Hyvärinen and Oja [HO97], is a fast fixed-point algorithm for solving the $k=1$ kurtosis maximization scheme for ICA. The problem at hand is

\max_{\|\bm{v}\|_{2}^{2}=1}\,\mathop{\mathrm{kurt}}(\bm{X}^{\top}\bm{v}).

(2.2.23)

First, we will perform some very basic analysis of this objective to verify its correctness. Notice by the change of variables $\bm{w}=\bm{U}^{\top}\bm{v}$ that this problem is equivalent to

\max_{\|\bm{w}\|_{2}^{2}=1}\,\mathrm{kurt}(\bm{Z}^{\top}\bm{w}).

This objective is simple enough that we can make strong statements about its correctness as a formulation for recovering the dictionary $\bm{U}$ . For example, in the population setting where $N\to\infty$ , we may use additivity properties of the kurtosis (Exercise 2.5) and our assumed normalization on the independent components to write the previous problem equivalently as

\max_{\|\bm{w}\|_{2}^{2}=1}\,\sum_{i=1}^{d}\mathrm{kurt}(z_{i})w_{i}^{4}.

(2.2.24)

It can be shown that under the Bernoulli-Gaussian assumption, the optimization landscape of this problem is “benign” (Exercise 2.7)—meaning that all local maxima of the objective function correspond to the recovery of one of the independent components. One efficient and scalable way to compute one of these maxima is via first-order optimization algorithms, which iteratively follow the gradient of the objective function and project onto the constraint set $\{\bm{w}\mid\|\bm{w}\|_{2}^{2}=1\}$ . Since we have assumed that each $z_{i}$ satisfies $\operatorname{Var}(z_{i})=1$ , we have for large $N$

\mathop{\mathrm{kurt}}(\bm{X}^{\top}\bm{u})\approx\tfrac{1}{N}\|\bm{X}^{\top}\bm{u}\|_{4}^{4}-3\|\bm{u}\|_{2}^{4}.

(2.2.25)

We can then derive a corresponding approximation to the gradient:

\nabla_{\bm{u}}\mathop{\mathrm{kurt}}(\bm{X}^{\top}\bm{u})\approx\tfrac{4}{N}\bm{X}(\bm{X}^{\top}\bm{u})^{\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}3}-12\|\bm{u}\|_{2}^{2}\bm{u}.

The FastICA algorithm uses a fixed-point method to compute directions of maximum kurtosis. It starts from the first-order optimality conditions for the kurtosis maximization problem, given the preceding gradient approximation and the constraint set, which read

\displaystyle\bm{X}(\bm{X}^{\top}\bm{u})^{\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}3}=\underbrace{\left\langle\bm{u},\bm{X}(\bm{X}^{\top}\bm{u})^{\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}3}\right\rangle}_{\lambda}\bm{u},

(2.2.26)

where the specific value of $\lambda$ is determined using the unit norm constraint on $\bm{u}$ . Exercise 2.6 describes the mathematical background necessary to derive these optimality conditions from first principles. Equation (2.2.26) is satisfied by any critical point of the kurtosis maximization problem; we want to derive an equation satisfied by only the maximizers. After noticing that $\lambda=\|\bm{X}^{\top}\bm{u}\|_{4}^{4}$ , we equivalently re-express (2.2.26) as the modified equation

\displaystyle\frac{1}{N}\bm{X}(\bm{X}^{\top}\bm{u})^{\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}3}-3\bm{u}=\left(\frac{\lambda}{N}-3\right)\bm{u},

(2.2.27)

and realize that any maximizer of (2.2.23) must satisfy $\lambda/N-3>0$ , assuming that $N$ is sufficiently large. Hence, we may normalize both sides of (2.2.27), giving the following fixed-point equation satisfied by every maximizer of (2.2.23):

\displaystyle\frac{\frac{1}{N}\bm{X}(\bm{X}^{\top}\bm{u})^{\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}3}-3\bm{u}}{\left\|\frac{1}{N}\bm{X}(\bm{X}^{\top}\bm{u})^{\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}3}-3\bm{u}\right\|_{2}}=\bm{u}.

(2.2.28)

Iterating the mapping defined by the lefthand side of this fixed point expression then yields the FastICA algorithm of Hyvärinen and Oja [HO97]:

\begin{split}\bm{v}^{+}&=\tfrac{1}{N}\bm{X}(\bm{X}^{\top}\bm{u})^{\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}3}-3\bm{u},\\ \bm{u}^{+}&=\bm{v}^{+}/\left\|\bm{v}^{+}\right\|_{2}.\end{split}

(2.2.29)

It turns out that the FastICA algorithm converges extremely rapidly (actually at a cubic rate) to columns of the dictionary $\bm{U}$ ; interested readers may consult [HO97] for details. Comparing the FastICA algorithm (2.2.29) to the power method studied in Section 2.1.1 for the PCA problem and the MSP algorithm (2.2.18), we notice a striking similarity. Indeed, FastICA is essentially a modified power method, involving the gradient of the empirical kurtosis rather than the simpler linear gradient of the PCA objective.

2.3 A Mixture of Overcomplete Low-Dimensional Subspaces

As we have seen, complete dictionary learning enjoys an elegant computational theory in which we maintain a symmetric autoencoding structure $\mathcal{E}(\bm{x})=\bm{U}^{\top}\bm{x}$ , $\mathcal{D}(\bm{z})=\bm{U}\bm{z}$ , with a scalable power-method-like algorithm (the MSP algorithm) for learning an orthogonal dictionary/codebook $\bm{U}$ from data. Nevertheless, for learning representations of general high-dimensional data distributions, one must expand the size of the codebook beyond the orthogonality requirement—meaning that we must have $\bm{A}\in\mathbb{R}^{D\times m}$ , with $m\gg D$ , corresponding to the case of an overcomplete dictionary/codebook,⁹⁹9We change the notation here from $\bm{U}$ to $\bm{A}$ in order to emphasize the non-orthogonality and non-square-shape of the overcomplete dictionary $\bm{A}$ . and the signal model

\bm{x}=\bm{A}\bm{z}+\bm{\varepsilon},\quad\|\bm{z}\|_{0}=d\ll m.

(2.3.1)

There are both geometric and physical/modeling motivations for passing to the overcomplete case. Geometrically, recall that in our original reduction from the mixture of subspace data model to the sparse dictionary model, a mixture of $K$ subspaces in $\mathbb{R}^{D}$ , each of dimension $d$ , led to a dictionary of shape $\bm{A}\in\mathbb{R}^{D\times Kd}$ . In other words, overcomplete dictionaries correspond to richer mixtures of subspaces, with more distinct modes of variability for modeling the high-dimensional data distribution. On the modeling side, we may run a computational experiment on real-world data that reveals the additional modeling power conferred by an overcomplete representation.

Example 2.1.

Given sampled images of hand-written digits, Figure 2.6(a) shows the result of fitting an orthogonal dictionary to the dataset. In contrast, Figure 2.6(b) shows the result of running an optimization algorithm for learning overcomplete dictionaries (which we will present in detail later in the Chapter) on these samples. Notice that the representations become far sparser and the codebooks far more interpretable—they consist of fundamental primitives for the strokes composing the digits, i.e. oriented edges. $\blacksquare$

In fact, overcomplete dictionary learning was originally proposed as a biologically plausible algorithm for image representation based on empirical evidence of how early stages of the visual cortex represent visual stimuli [OF96, OF97].

In the remainder of this section, we will overview the conceptual and computational foundations of overcomplete dictionary learning. Supposing that the model (2.3.1) is satisfied with sparse codes $\bm{z}$ , overcomplete dictionary $\bm{A}$ , and sparsity level $d$ , and given samples $\bm{X}=[\bm{x}_{1},\dots,\bm{x}_{N}]$ of $\bm{x}$ , we want to learn an encoder $\mathcal{E}:\mathbb{R}^{D}\to\mathbb{R}^{m}$ mapping each $\bm{x}$ to its sparse code $\bm{z}$ , and a decoder $\mathcal{D}(\bm{z})=\bm{A}\bm{z}$ reconstructing each $\bm{x}$ from its sparse code. In diagram form:

\bm{x}\xrightarrow{\hskip 11.38109pt\mathcal{E}\hskip 11.38109pt}\bm{z}\xrightarrow{\hskip 5.69054pt\mathcal{D}=\bm{A}\hskip 5.69054pt}\hat{\bm{x}}.

(2.3.2)

We will start from the construction of the encoder $\mathcal{E}$ . We will work incrementally: first, given the true dictionary $\bm{A}$ , we will show how the problem of sparse coding gives an elegant, scalable, and provably-correct algorithm for recovering the sparse code $\bm{z}$ of $\bm{x}$ . Although this problem is NP-hard in the worst case, it can be solved efficiently and scalably for dictionaries $\bm{A}$ which are incoherent, i.e. having columns that are not too correlated. The encoder architecture encompassed by this solution will no longer be symmetric: we will see it has the form of a primitive deep network, which depends on the dictionary $\bm{A}$ .

Then we will proceed to the task of learning the decoder $\mathcal{D}$ , or equivalently the overcomplete dictionary $\bm{A}$ . We will derive an algorithm for overcomplete dictionary learning that allows us to simultaneously learn the codebook $\bm{A}$ and the sparse codes $\bm{z}$ , using ideas from sparse coding. Finally, we will discuss a more modern perspective on learnable sparse coding that leads us to a fully asymmetric encoder-decoder structure, as an alternative to (2.3.2). Here, the decoder will correspond to an incremental solution to the sparse dictionary learning problem, and yield a pair of deep network-like encoder decoders for sparse dictionary learning. This structure will foreshadow many developments to come in the remainder of the monograph, as we progress from analytical models to modern neural networks.

2.3.1 Sparse Coding with an Overcomplete Dictionary

In this section, we will consider the data model (2.3.1), which accommodates sparse linear combinations of many motifs, or atoms. Given data $\{\bm{x}_{i}\}_{i=1}^{N}\subseteq\mathbb{R}^{D}$ satisfying this model, i.e. expressible as

\bm{x}_{i}=\bm{A}\bm{z}_{i}+\bm{\varepsilon}_{i},\qquad\forall i\in[N]

(2.3.3)

for some dictionary $\bm{A}\in\mathbb{R}^{D\times m}$ with $m$ atoms, sparse codes $\bm{z}_{i}$ such that $\|\bm{z}_{i}\|_{0}\leq d$ , and small errors $\bm{\varepsilon}_{i}$ , the sparse coding problem is to recover the codes $\bm{z}_{i}$ as accurately as possible from the data $\bm{x}_{i}$ , given the dictionary $\bm{A}$ . Efficient algorithms to solve this problem succeed when the dictionary $\bm{A}$ is incoherent in the sense that the inner products $\bm{a}_{i}^{\top}\bm{a}_{j}$ are uniformly small, hence the atoms are nearly orthogonal.¹⁰¹⁰10As it turns out, in a high-dimensional space, it is rather easy to pack a number of nearly orthogonal vectors that is much larger than the ambient dimension [WM22].

Note that we can collect the $\bm{x}_{i}$ into $\bm{X}=\begin{bmatrix}\bm{x}_{1},\dots,\bm{x}_{N}\end{bmatrix}\in\mathbb{R}^{D\times N}$ , collect the $\bm{z}_{i}$ into $\bm{Z}=\begin{bmatrix}\bm{z}_{1},\dots,\bm{z}_{N}\end{bmatrix}\in\mathbb{R}^{d\times N}$ , and collect the $\bm{\varepsilon}_{i}$ into $\bm{E}=\begin{bmatrix}\bm{\varepsilon}_{1},\dots,\bm{\varepsilon}_{N}\end{bmatrix}\in\mathbb{R}^{D\times N}$ , to rewrite (2.3.3) as

\bm{X}=\bm{A}\bm{Z}+\bm{E}.

(2.3.4)

A natural approach to solving the sparse coding problem is to seek the sparsest signals that are consistent with our observations, and this naturally leads to the following optimization problem:

\min_{\bm{Z}\in\mathbb{R}^{d\times N}}\left\{\|\bm{X}-\bm{A}\bm{Z}\|_{F}^{2}+\lambda\|\bm{Z}\|_{1}\right\},

(2.3.5)

where the $\ell^{1}$ norm $\|\bm{Z}\|_{1}$ , taken elementwise, is known to promote sparsity of the solution [WM22]. This optimization problem is known as the LASSO problem.

However, unlike PCA or the complete dictionary learning problem, there is no clear power iteration-type algorithm to recover $\bm{Z}^{\star}$ . A natural alternative is to consider solving the above optimization problem with gradient descent type algorithms. Let $f(\bm{Z})=\|\bm{X}-\bm{A}\bm{Z}\|_{2}^{2}+\lambda\|\bm{Z}\|_{1}$ . Conceptually, we could try to find $\bm{Z}^{\star}$ with the following iterations:

\bm{Z}_{t+1}\leftarrow\bm{Z}_{t}+\eta\nabla f(\bm{Z}_{t}).

(2.3.6)

However, because the term associated with the $\ell^{1}$ norm $\|\bm{Z}\|_{1}$ is non-smooth, we cannot just run gradient descent. For this type of function, we need to replace the gradient $\nabla f(\bm{Z})$ with something that generalizes the notion of gradient, known as the subgradient $\partial f(\bm{Z})$ :

\bm{Z}_{t+1}\leftarrow\bm{Z}_{t}+\eta\partial f(\bm{Z}_{t}).

(2.3.7)

However, it is known that the convergence of subgradient descent is usually very slow. Hence, for this type of optimization problems, it is conventional to adopt a so-called proximal gradient descent-type algorithm. We give a technical overview of this method in Section A.1.3.

Applying proximal gradient to the LASSO objective function (2.3.5) leads to the classic iterative shrinkage-thresholding algorithm (ISTA), which implements the iteration

	$\displaystyle\bm{Z}_{1}$	$\displaystyle\sim$	$\displaystyle\operatorname{\mathcal{N}}(\bm{0},\bm{I}),$		(2.3.8)
	$\displaystyle\bm{Z}_{t+1}$	$\displaystyle=$	$\displaystyle S_{\eta\lambda}\left(\bm{Z}_{t}-2\eta\bm{A}^{\top}(\bm{A}\bm{Z}_{t}-\bm{X})\right),\quad\forall t\geq 1,$		(2.3.9)

with step size $\eta\geq 0$ , and the soft-thresholding operator $S_{\alpha}$ defined on scalars by

	$\displaystyle S_{\alpha}(x)$	$\displaystyle\doteq\begin{cases}x-\alpha,&x\geq\alpha,\\ 0,&-\alpha<x<\alpha,\\ x+\alpha,&x\leq-\alpha\end{cases}$		(2.3.10)
		$\displaystyle=\operatorname{sign}(x)\max\{\lvert x\rvert-\alpha,0\},$		(2.3.11)

and applied element-wise to the input matrix. As a proximal gradient descent algorithm applied to a convex problem, convergence to a global optimum is assured, and a precise convergence rate can be derived straightforwardly [WM22].

We now take a moment to remark on the functional form of the update operator in (2.3.9). It takes the form

\bm{Z}_{t+1}=\texttt{nonlinearity}(\bm{Z}_{t}+\texttt{linear}^{\top}(\texttt{linear}(\bm{Z}_{t})+\texttt{bias})).

(2.3.12)

This functional form is very similar to that of a residual network layer, namely,

\bm{Z}_{t+1}=\bm{Z}_{t}+\texttt{linear}_{1}^{\top}(\texttt{nonlinearity}(\texttt{linear}_{2}(\bm{Z}_{t})+\texttt{bias}_{1})+\texttt{bias}_{2},

(2.3.13)

only decoupling the linear maps (i.e., matrix multiplications), adding a bias, and moving the nonlinearity. The nonlinearity in (2.3.9) is notably similar to the commonly-used ReLU activation function $x\mapsto\max\{x,0\}$ in deep learning—in particular, the soft-thresholding operator is like a “signed” ReLU activation, which is also shifted by a bias. The ISTA, then, can be viewed as a forward pass of a primitive (recurrent one-layer) neural network. We argue in Chapter 4 that such operations are essential to deep representation learning.

2.3.2 Overcomplete Dictionary Learning

Recall that we have the data model

\bm{X}=\bm{A}\bm{Z}+\bm{E},

(2.3.14)

where $\bm{Z}$ is sparse, and our goal previously was to estimate $\bm{Z}$ given knowledge of the data $\bm{X}$ and the dictionary atoms $\bm{A}$ . Now we turn to the more practical and more difficult case where we do not know either $\bm{A}$ or $\bm{Z}$ and seek to learn them from a large dataset.

A direct generalization of Equation 2.3.5 suggests solving the problem

\min_{\tilde{\bm{A}},\tilde{\bm{Z}}}\left\{\|\bm{X}-\tilde{\bm{A}}\tilde{\bm{Z}}\|_{F}^{2}+\lambda\|\tilde{\bm{Z}}\|_{1}\right\}.

(2.3.15)

However, the bilinear term in Equation 2.3.15 introduces a scale ambiguity that hinders convergence: given any point $(\tilde{\bm{A}},\tilde{\bm{Z}})$ and any constant $c>0$ , the substitution $(c^{-1}\tilde{\bm{A}},c\tilde{\bm{Z}})$ gives loss value

\|\bm{X}-\tilde{\bm{A}}\tilde{\bm{Z}}\|_{F}^{2}+c\lambda\|\tilde{\bm{Z}}\|_{1}.

(2.3.16)

This loss is evidently minimized over $c$ by taking $c\to 0$ , which corresponds to making the rescaled dictionary $c^{-1}\tilde{\bm{A}}$ go ‘to infinity’. In particular, the iterates of any optimization algorithm solving (2.3.15) will not converge.

This issue with (2.3.15) is dealt with by adding a constraint on the scale of the columns of the dictionary $\tilde{\bm{A}}$ . For example, it is common to assume that each column $\bm{A}_{j}$ of the dictionary $\bm{A}$ in (2.3.14) has bounded $\ell^{2}$ norm—say, without loss of generality, by $1$ . We then enforce this as a constraint, giving the objective

\min_{\tilde{\bm{Z}},\tilde{\bm{A}}\,:\,\|\tilde{\bm{A}}_{j}\|_{2}\leq 1}\left\{\|\bm{X}-\tilde{\bm{A}}\tilde{\bm{Z}}\|_{F}^{2}+\lambda\|\tilde{\bm{Z}}\|_{1}\right\}.

(2.3.17)

Constrained optimization problems such as (2.3.17) can be solved by a host of sophisticated algorithms [NW06]. However, a simple and scalable one is actually furnished by the same proximal gradient descent algorithm that we used to solve the sparse coding problem in the previous section. We can encode each constraint as additional regularization term, via the characteristic function for the constraint set—details are given in Example A.1. Applying proximal gradient descent to the resulting regularized problem is equivalent to projected gradient descent, in which, at each iteration, the iterates after taking a gradient descent step are projected onto the constraint set.

Remark 2.11 ( $\ell^{4}$ maximization versus $\ell^{1}$ minimization).

Note that the above problem formulation follows naturally from the LASSO formulation (2.3.5) for sparse coding. We promote the sparsity of the solution via the $\ell^{1}$ norm. Nevertheless, if we are only interested in recovering the over-complete dictionary $\bm{A}$ , the $\ell^{4}$ maximization scheme introduced in Section 2.2.2 also generalizes to the over-complete case, without any significant modification. Interested readers may refer to the work of [QZL+20].

The above problem (2.3.17), which we call overcomplete dictionary learning, is nonconvex as here both $\bm{A}$ and $\bm{Z}$ are unknown. It cannot be solved easily by the standard convex optimization toolkit. Nevertheless, because it is interesting, simple to state, and practically important, there have been many important works dedicated to different algorithms and theoretical analysis for this problem. Here, for the interest of this manuscript, we present an idiomatic approach to solve this problem which is closer to the spirit of deep learning.

From our experience with the LASSO problem above, it is easy to see that, for the two unknowns $\bm{A}$ and $\bm{Z}$ , if we fix one and optimize the other, each subproblem is in fact convex and easy to solve. This naturally suggests that we could attempt to solve the above program (2.3.17) by minimizing against $\bm{Z}$ or $\bm{A}$ alternatively, say using gradient descent. Coupled with a natural choice of initialization, this leads to the following iterative scheme:

	$\displaystyle\bm{Z}^{\ell+1}=S_{\eta\lambda}\left(\bm{Z}^{\ell}-2\eta\bm{A}_{+}^{\top}(\bm{A}_{+}\bm{Z}^{\ell}-\bm{X})\right),\quad\bm{Z}^{1}=\mathbf{0},\quad\forall\ell\in[L]$		(2.3.18)
	$\displaystyle\bm{Z}^{+}=\bm{Z}^{L},$		(2.3.19)
	$\displaystyle\bm{A}_{t+1}=\operatorname{proj}_{\\|(\,\cdot\,)_{j}\\|_{2}\leq 1,\,\forall j}\left(\bm{A}_{t}-2\nu(\bm{A}_{t}\bm{Z}^{+}-\bm{X})(\bm{Z}^{+})^{\top}\right),\quad(\bm{A}_{1})_{j}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}\operatorname{\mathcal{N}}(\bm{0},\tfrac{1}{D}\bm{I}),\enspace\forall j\in[m],\quad\;\;\forall t\in[T],$		(2.3.20)
	$\displaystyle\bm{A}_{+}=\bm{A}_{T},$		(2.3.21)

where the projection operation in the update for $\bm{A}$ ensures each column has at most unit $\ell^{2}$ norm, via $\bm{A}_{j}\mapsto\bm{A}_{j}/\max\{\|\bm{A}_{j}\|_{2},1\}$ , and where $\bm{A}_{+}$ is initialized with each column i.i.d. $\mathcal{N}(\mathbf{0},\tfrac{1}{D}\bm{I})$ . The above consists of one ‘block’ of alternating minimization, and we repeatedly perform such blocks, each with independent initializations, until convergence. Above, we have used two separate indices $\{t\}$ and $\{\ell\}$ to indicate the iterations. As we will see later, this allows us to interpret the two updates separately in the context of deep learning.

Despite the dictionary learning problem being a nonconvex problem, it has been shown that alternating minimization type algorithms indeed converge to the correct solution, at least locally. See, for example, [AAJ+16]. As a practical demonstration, the above algorithm (with $L=T=1$ ) was used to generate the results for overcomplete dictionary learning in Figure 2.6.

2.3.3 Learned Deep Sparse Coding

The main insight from the alternating minimization algorithm for overcomplete dictionary learning in the previous section (Equations 2.3.18 and 2.3.20) is to notice that when we fix $\bm{A}$ , the ISTA update for $\bm{Z}^{\ell}$ (2.3.18) looks like the forward pass of a deep neural network with weights given by $\bm{A}$ (and $\bm{A}^{\top}$ ). But in general, we do not know the true $\bm{A}$ , and the current estimate $\bm{A}_{+}$ could be erroneous. Hence it needs to be further updated using (2.3.20) based on the residual of using the current estimate of the sparse codes $\bm{Z}^{+}$ to reconstruct $\bm{X}$ . The alternating minimization algorithm iterates these two procedures until convergence. But we can instead extrapolate, and design other learning procedures by combining these insights with techniques from deep learning. This leads to more interpretable network architectures, which will be a recurring theme throughout this manuscript.

Learned ISTA.

The above deep-network interpretation of the alternating minimization is more conceptual than practical, as the process could be rather inefficient and take many layers or iterations to converge. But this is mainly because we try to infer both $\bm{Z}$ and $\bm{A}$ from $\bm{X}$ . The problem can be significantly simplified and the above iterations can be made much more efficient in the supervised setting, where we have a dataset of input and output pairs $(\bm{X},\bm{Z})$ distributed according to (2.3.14) and we only seek to learn $\bm{A}^{\ell}$ for the layerwise learnable sparse coding iterations (2.3.29):

\displaystyle\bm{Z}^{\ell+1}=S_{\eta\lambda}\left(\bm{Z}^{\ell}-2\eta(\bm{A}^{\ell})^{\top}(\bm{A}^{\ell}\bm{Z}^{\ell}-\bm{X})\right),\quad\forall\ell\in[L].

(2.3.22)

If we denote the operator for each iteration as $\bm{Z}^{\ell+1}=f(\bm{A}^{\ell},\bm{Z}^{\ell})$ , the above iteration can be illustrated in terms of a diagram:

\bm{X},\bm{Z}^{1}\xrightarrow{\hskip 2.84526ptf(\bm{A}^{1},\,\cdot\,)\hskip 2.84526pt}\bm{Z}^{2}\xrightarrow{\hskip 2.84526ptf(\bm{A}^{2},\,\cdot\,)\hskip 2.84526pt}\bm{Z}^{3}\xrightarrow{\hskip 2.84526ptf(\bm{A}^{3},\,\cdot\,)\hskip 2.84526pt}\cdots\bm{Z}^{L}\xrightarrow{\hskip 2.84526ptf(\bm{A}^{L},\,\cdot\,)\hskip 2.84526pt}\bm{Z}^{L+1}\approx\bm{Z}.

Thus, given the sequential architecture, to learn the operator $\bm{A}^{\ell}$ at each layer, it is completely natural to learn it, say via back propagation (BP),¹¹¹¹11See Section A.2 for a brief description of BP. by minimizing the error between the final code $\bm{Z}^{L}$ and the ground truth $\bm{Z}$ :

\min_{\{\bm{A}^{\ell}\}}\big{\|}\bm{Z}^{L}(\bm{A}^{1},\ldots,\bm{A}^{L})-\bm{Z}\big{\|}_{2}^{2}.

(2.3.23)

This is the basis of the Learned ISTA (LISTA) algorithm [GL10], which can be viewed as the learning algorithm for a deep neural network, which tries to emulate the sparse encoding process from $\bm{X}$ to $\bm{Z}$ . In particular, it can be viewed as a simple representation learning algorithm. In fact, this same methodology can be used as a basis to understand the representations computed in more powerful network architectures, such as transformers. We develop these implications in detail in Chapter 4.

Sparse Autoencoders.

The original motivation for overcomplete dictionary learning was to provide a simple generative model for high-dimensional data. We have seen with LISTA that, in addition, iterative algorithms for learning sparsely-used overcomplete dictionaries provide an interpretation for ReLU-like deep networks, which we will generalize in later chapters to more complex data distributions than (2.3.14). But it is also worth noting that even in the modern era of large models, the data generating model (2.3.14) provides a useful practical basis for interpreting features in pretrained large-scale deep networks, such as transformers, following the hypothesis that the (non-interpretable, a priori) features in these networks consist of sparse “superpositions” of underlying features, which are themselves interpretable [EHO+22a]. These unsupervised learning paradigms are generally more data friendly than LISTA, as well, which requires large amounts of labeled $(\bm{X},\bm{Z})$ pairs for supervised training.

We can use our development of the LISTA algorithm above to understand common practices in this field of research. In the most straightforward instantiation (see [HCS+24, GTT+25]), a large number of features from a pretrained deep network $h$ are collected from different inputs $\bm{x}_{i}$ , which themselves are chosen based on a desired interpretation task.¹²¹²12For example, the inputs $\bm{x}_{i}$ could correspond to texts containing samples of computer code in different programming languages, with our task being to try to identify interpretable features in a transformer feature map $h$ corresponding to different salient aspects of the input, such as the specific programming language (distinct across input “classes”) or the need to insert a matching parenthesis at the current position (common across input “classes”). We discuss the use of deep networks, and in particular transformers, for text representation learning in greater detail in Chapter 7. For simplicity, we will use $h$ to denote the pre-selected feature map in question, with $D$ -dimensional features; given $N$ sample inputs, let $\bm{H}\in\mathbb{R}^{D\times N}$ denote the full matrix of features of $h$ . Then a so-called sparse autoencoder $f:\mathbb{R}^{D}\to\mathbb{R}^{d}$ , with decoder $g:\mathbb{R}^{d}\to\mathbb{R}^{D}$ , is trained via the LASSO loss (2.3.5):

\min_{f,g}\|\bm{H}-g(f(\bm{H})))\|_{F}^{2}+\lambda\|f(\bm{H})\|_{1},

(2.3.24)

where the sparse autoencoder $f$ takes the form of a one-layer neural network, i.e. $f(\bm{h}_{i})=\sigma(\bm{W}_{\mathrm{enc}}(\bm{h}_{i}-\bm{b})+\bm{b}_{\mathrm{enc}})$ , where $\sigma(x)=\max\{x,0\}$ is the ReLU activation function, and the decoder $g$ is linear, so that $g(\bm{z}_{i})=\bm{W}_{\mathrm{dec}}\bm{z}+\bm{b}$ .

The parameterization and training procedure (2.3.24) may initially seem to be an arbitrary application of deep learning to the sparse coding problem, but it is actually highly aligned with the algorithms we have studied above for layerwise sparse coding with a learned dictionary. In particular, recall the LISTA architecture $\bm{Z}^{L}=f(\bm{A}^{L},f(\bm{A}^{L-1},\cdots,f(\bm{A}^{1},\bm{X})\cdots))$ . In the special case $L=2$ , we have

\bm{Z}^{2}=f(\bm{A}^{1},\bm{X})=S_{\eta\lambda}\left(\bm{Z}^{1}-2\eta(\bm{A}^{1})^{\top}(\bm{A}^{1}\bm{Z}^{1}-\bm{X})\right).

(2.3.25)

Let us assume that the sparse codes $\bm{Z}$ in question are nonnegative, i.e., that $\bm{Z}\geq\mathbf{0}$ .¹³¹³13In the data generating model (2.3.14), an arbitrary dictionary-and-sparse-code pair $(\bm{A},\bm{Z})$ can be replaced by one in which $\bm{Z}\geq\mathbf{0}$ simply by doubling the number of columns in $\bm{A}$ , so from a modeling perspective, this is a very mild assumption. Then (see Example A.3), we can consider an equivalent LISTA architecture obtained from the sparse coding objective with an additional nonnegativity constraint on $\bm{Z}$ as

\bm{Z}^{2}=f(\bm{A}^{1},\bm{X})=\max\left\{\bm{Z}^{1}-2\eta(\bm{A}^{1})^{\top}(\bm{A}^{1}\bm{Z}^{1}-\bm{X})-\lambda\eta\mathbf{1},0\right\},

(2.3.26)

and after some algebra, express this as

\bm{Z}^{2}=f(\bm{A}^{1},\bm{X})=\max\left\{2\eta(\bm{A}^{1})^{\top}+\left(\bm{Z}^{1}-2\eta(\bm{A}^{1})^{\top}\bm{A}^{1}\bm{Z}^{1}-\lambda\eta\mathbf{1}\right),0\right\}.

(2.3.27)

Given the ability to change the sparse code initialization $\bm{Z}^{1}$ as a learnable parameter (which, in the current framework, must have all columns equal to the same learnable vector), this has the form of a ReLU neural network with learnable bias—identical to the sparse autoencoder $f$ ! Moreover, to decode the learned sparse codes $\bm{Z}^{2}$ , it is natural to apply the learned dictionary $\bm{Z}^{2}\mapsto\bm{A}^{1}\bm{Z}^{2}$ . Then the only difference between this and the SAE decoder $g$ is the additional bias $\bm{b}$ , which can technically be absorbed into $\bm{H}$ and $f$ in the training objective (2.3.24).

Thus, the SAE parameterization and training procedure coincides with LISTA training with $L=1$ , and a modified training objective—using the LASSO objective (2.3.5), which remains unsupervised, instead of the supervised reconstruction loss (2.3.23) used in vanilla LISTA. In particular, we can understand the SAE architecture in terms of our interpretation of the LISTA architecture in terms of layerwise sparse coding in (2.3.29). This connection is suggestive of a host of new design strategies for improving practical interpretability methodology, many of which remain tantalizingly unexplored. We begin to lay out some connections to broader autoencoding methodology in Chapter 5.

Layerwise learned sparse coding?

In the supervised setting, LISTA provides a deep neural network analogue of the sparse coding iteration, with layerwise-learned dictionaries, inspired by alternating minimization; even in the unsupervised setting, the same methodology can be applied to learning, as with sparse autoencoders. But the connection between low-dimensional-structure-seeking optimization algorithms and deep network architectures goes much deeper than this, and suggests an array of scalable and natural neural learning architectures which may even be usable without backpropagation.

As a simple illustration, we return the alternating minimization iterations (2.3.18) and (2.3.20). This scheme randomly re-initializes the dictionary $\bm{A}_{1}$ on every such update. An improvement uses instead warm starting, where the residual is generated using the previous estimate $\bm{A}_{+}$ for the dictionary. If we then view each ISTA update (2.3.18) as a layer and allow the associated dictionary, now coupled with the sparse code updates as $\bm{A}_{\ell}$ , to update in time, this leads to a “layerwise-learnable” sparse coding scheme:

	$\displaystyle\bm{Z}^{1}=\mathbf{0},\quad(\bm{A}_{1})_{j}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}\operatorname{\mathcal{N}}(\bm{0},\tfrac{1}{D}\bm{I}),\enspace\forall j\in[m],$		(2.3.28)
	$\displaystyle\bm{Z}^{\ell+1}=S_{\eta\lambda}\left(\bm{Z}^{\ell}-2\eta(\bm{A}_{\ell})^{\top}(\bm{A}_{\ell}\bm{Z}^{\ell}-\bm{X})\right),$		(2.3.29)
	$\displaystyle\bm{A}_{\ell+1}=\bm{A}_{\ell}-2\nu(\bm{A}_{\ell}\bm{Z}^{\ell+1}-\bm{X})(\bm{Z}^{\ell+1})^{\top}.$		(2.3.30)

Note that this iteration corresponds to a relabeling of (2.3.18) and (2.3.20) for $T=L=1$ , over infinitely many blocks. Each of the ‘inner’ steps updating $\bm{Z}$ can be considered as a one-layer forward pass, while each of the ‘outer’ steps updating $\bm{A}$ can be considered as a one-layer backward pass, of a primitive deep neural network. In particular, this algorithm is the simplest case in which a clear divide between forward optimization and backward learning manifests. This divide is still observed in current neural networks and autoencoders—we will have much more to say about it in Chapter 4 and in Chapter 5.

Notice that the above layer-wise scheme also suggests a plausible alternative to the current end-to-end optimization strategy that primarily relies on back propagation (BP) detailed in Section A.2.3. Freeing training large networks from BP would be one of the biggest challenges and opportunities in the future, as we will discuss more at the end of the book in Chapter 8.

2.4 Summary and Notes

The idealistic models we have presented in this chapter—PCA, ICA, and dictionary learning—were developed over the course of the twentieth century. Many books have been written solely about each method, so we will only attempt here to give a broad overview of the key works and history.

Jolliffe [Jol86] attributes principal component analysis to Pearson [Pea01], and independently Hotelling [Hot33]. In mathematics, the main result on the related problem of low-rank approximation in unitarily invariant norms is attributed to Eckart and Young [EY36], and to Mirsky for full generality [Mir60]. PCA continues to play an important role in research as perhaps the simplest model problem for unsupervised representation learning: as early as the 1980s, works such as [Oja82] and [BH89] used the problem to understand learning in primitive neural networks, and more recently, it has served as a tool for understanding more complex representation learning frameworks, such as diffusion models [WZZ+24].

Independent component analysis was proposed by [BJC85] and pioneered by Aapo Hyvärinen in the 1990s and early 2000s in a series of influential works: see [HO00a] for a summary. As a simple model for structure that arises in practical data, it initially saw significant use in applications such as blind source separation, where each independent component $z_{i}$ represents an independent source (such as sound associated to a distinct instrument in a musical recording) that is superimposed to produce the observation $\bm{x}=\bm{U}\bm{z}$ .

The problem of dictionary learning can, in the complete or orthogonal case, be seen as one of the foundational problems of twentieth-century signal processing, particularly in linear systems theory, where the Fourier basis plays the key role; from the 1980s onward, the field of computational harmonic analysis crystallized around the study of alternate such dictionaries for classes of signals in which optimal approximation could only be realized in a basis other than Fourier (e.g., wavelets) [DVD+98]. However, the importance of the case of redundant bases, or overcomplete dictionaries, was only highlighted following the pioneering work of Olshausen and Field [OF96, OF97]. Early subsequent work established the conceptual and algorithmic foundations for learning sparsely-used overcomplete dictionaries, often aimed at representing natural images [Don01, DM03, EA06, MK07, AEB06, MBP14, GJB15]. Later, a significant amount of theoretical interest in the problem, as an important and nontrivial model problem for unsupervised representation learning, led to its study by the signal processing, theoretical machine learning, and theoretical computer science communities, in particular focused on conditions under which the problem could be provably and efficiently solved. A non-exhaustive list of notable works in this line include those of [SWW12] and [SQW17] on the complete case; [AGM+15]; [BKS15]; and [QZL+20]. Many deep theoretical questions about this simple-to-state problem remain open, perhaps in part due to a tension with the problem’s worst-case NP-hardness (e.g., see [Til15]).

Table 2.1: Summary of (generalized) power methods presented in the Chapter

Problem	Algorithm	Iteration	Type of Structure Enforced
PCA	Power Method	$\bm{u}_{t}=\frac{\bm{X}\bm{X}^{\top}\bm{u}_{t}}{\\|\bm{X}\bm{X}^{\top}\bm{u}_{t}\\|_{2}}$	$1$ -dim. subspace (unit vector)
ICA	FastICA	$\bm{u}_{t+1}=\frac{\tfrac{1}{N}\bm{X}(\bm{X}^{\top}\bm{u}_{t})^{\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}3}-3\bm{u}_{t}}{\left\\|\tfrac{1}{N}\bm{X}(\bm{X}^{\top}\bm{u}_{t})^{\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}3}-3\bm{u}_{t}\right\\|_{2}}$	$1$ -dim. subspace (unit vector)
Complete DL	MSP Algorithm	$\bm{U}_{t+1}=\mathcal{P}_{\mathrm{O}(D)}[({\bm{U}_{t}\bm{X}})^{\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}3}\bm{X}^{\top}]$	$D$ -dim. subspace (orthogonal matrix)

One point that we wish to highlight from the study of these classical analytical models for low-dimensional structure is the common role played by various generalized power methods—algorithms that very rapidly converge, at least locally, to various types of low-dimensional structures. The terminology for this class of algorithms follows the work of [MYP+10]. At a high level, modeled on the classical power iteration for computation of the top eigenvector of a semidefinite matrix $\bm{A}\in\mathbb{R}^{n\times n}$ , that is

\bm{u}_{t+1}=\frac{\bm{A}\bm{u}_{t}}{\|\bm{A}\bm{u}_{t}\|_{2}},

(2.4.1)

this class of algorithms consists of a “powering” operation involving a matrix $\bm{A}$ associated to the data, along with a “projection” operation that enforces a desired type of structure. Table 2.1 presents a summary of the algorithms we have studied in this chapter that follow this structure. The reader may appreciate the applicability of this methodology to different types of low-dimensional structure, and different losses (i.e., both the quadratic loss from PCA, and the kurtosis-type losses from ICA), as well as the lack of such an algorithm for overcomplete dictionary learning, despite the breadth of the literature on these algorithms. We see the development of power methods for further families of low-dimensional structures, particularly those relevant to applications where deep learning is prevalent, as one of the more important (and open) research questions suggested by this chapter.

The connection we make in Section 2.2.1 between the geometric mixture-of-subspaces distributional assumption and the more analytically-convenient sparse dictionary assumption has been mentioned in prior work, especially by those focused on generalized principal component analysis and applications such as subspace clustering, e.g., work of [VMS16]. The mixture of subspaces assumption will continue to play a significant role throughout this manuscript, both as an analytical test case for different algorithmic paradigms, and as a foundation for deriving different deep network architectures, as with LISTA in Section 2.3.3, but which can scale to more complex data distributions.

2.5 Exercises and Extensions

Exercise 2.1.

Prove that, for any symmetric matrix $\bm{A}$ , the solution to the problem $\max_{\bm{U}\in\mathsf{O}(D,d)}\operatorname{tr}\left(\bm{U}^{\top}\bm{A}\bm{U}\right)$ is the matrix $\bm{U}^{\star}$ whose columns are the top $d$ unit eigenvectors of $\bm{A}$ .

Exercise 2.2.

Let $\bm{z}\sim\mathcal{N}(\mathbf{0},\sigma^{2}\bm{I})$ be a Gaussian random variable with independent components, each with variance $\sigma^{2}$ . Prove that for any orthogonal matrix $\bm{Q}$ (i.e., $\bm{Q}^{\top}\bm{Q}=\bm{I}$ ), the random variable $\bm{Q}\bm{z}$ is distributed identically to $\bm{z}$ . (Hint: recall the formula for the Gaussian probability density function, and the formula for the density of a linear function of a random variable.)

Exercise 2.3.

The notion of statistical identifiability discussed above can be related to symmetries of the model class, allowing estimation to be understood in a purely deterministic fashion without any statistical assumptions.

Consider the model $\bm{X}=\bm{U}\bm{Z}$ for matrices $\bm{X},\bm{U},\bm{Z}$ of compatible sizes.

1.

Show that if $\bm{A}$ is any square invertible matrix of compatible size, then the pair $(\bm{U}\bm{A},\bm{A}^{-1}\bm{Z})$ also equals $\bm{X}$ under the model. We call this a $\mathsf{GL}(d)$ symmetry.
2.

Suppose each column of $\bm{Z}$ is an independent and identically distributed observation from a common statistical model $\bm{z}$ , which moreover has zero mean and independent components $z_{i}$ with positive variance. Show that for any square invertible matrix $\bm{A}$ , if $\bm{A}\bm{z}$ has uncorrelated components, then $\bm{A}$ can be written as $\bm{D}_{1}\bm{Q}\bm{D}_{2}$ , where $\bm{Q}$ is an orthogonal matrix and $\bm{D}_{1}$ , $\bm{D}_{2}$ are diagonal matrices. This links the “independence” assumption in ICA to a “symmetry breaking” effect, which only allows scale and rotational symmetries.

Exercise 2.4.

Consider the model $\bm{x}=\bm{U}\bm{z}$ , where $\bm{U}\in\mathbb{R}^{D\times d}$ with $D\geq d$ is fixed and has rank $d$ , and $\bm{z}$ is a zero-mean random variable. Let $\bm{x}_{1},\dots,\bm{x}_{N}$ denote i.i.d. observations from this model.

1.

Show that the matrix $\bm{X}=[\bm{x}_{1},\dots,\bm{x}_{N}]$ has rank no larger than $d$ , and therefore there is an orthonormal matrix $\bm{V}\in\mathbb{R}^{D\times d}$ so that $\bm{X}=\bm{V}\bm{Y}$ , where $\bm{Y}\in\mathbb{R}^{d\times N}$ . (Hint: use PCA.)
2.

Show that the whitened matrix $(\bm{Y}\bm{Y}^{\top})^{-1/2}\bm{Y}$ exists in expectation whenever $\operatorname{Cov}(\bm{z})$ is nonsingular, and that it has identity empirical covariance.¹⁴¹⁴14In particular, it can be proved mathematically that this is enough to guarantee that the whitened matrix exists with high probability whenever $\bm{z}$ satisfies a suitable concentration inequality and $N$ is sufficiently large.
3.

Show, by using the singular value decomposition of $\bm{U}$ , that the matrix $\bm{V}$ can be chosen so that the whitened matrix satisfies $(\bm{Y}\bm{Y}^{\top})^{-1/2}\bm{Y}=\bm{W}[\bm{z}_{1},\dots,\bm{z}_{N}]$ , where $\bm{W}$ is an orthonormal matrix.

Exercise 2.5.

Let $X$ and $Y$ be zero-mean independent random variables.

1.

Show that $\mathop{\mathrm{kurt}}(X+Y)=\mathop{\mathrm{kurt}}(X)+\mathop{\mathrm{kurt}}(Y)$ .
2.

For any $\alpha\in\mathbb{R}$ , show that $\mathop{\mathrm{kurt}}(\alpha X)=\alpha^{4}\mathop{\mathrm{kurt}}(X)$ .

Exercise 2.6.

Let $f:\mathbb{R}^{d}\to\mathbb{R}$ be a given twice-continuously-differentiable objective function. Consider the spherically-constrained optimization problem

\max_{\|\bm{u}\|_{2}^{2}=1}\,f(\bm{u}).

(2.5.1)

In this exercise, we will derive the expressions we gave in the FastICA derivation for maximizing kurtosis over the sphere via a gradient ascent algorithm. These expressions are special cases of a rich theory of calculus and optimization on manifolds, of which the sphere is a particular example. A deep technical study of this field is out-of-scope for our purposes, so we only mention two key references for the interested reader: the pioneering textbook by Absil, Mahony, and Sepulchre [AMS09], and a more recent introductory treatise by Boumal [Bou23].

For any constraint set $\mathcal{M}$ that is a differentiable submanifold of $\mathbb{R}^{d}$ , the tangent space at a point $\bm{u}\in\mathcal{M}$ is, informally, the best local linear approximation to the manifold $\mathcal{M}$ at the point $\bm{u}$ . In the important special case where $\mathcal{M}$ is defined locally at $\bm{u}$ as a level set of a function $F:\mathbb{R}^{d}\to\mathbb{R}$ , that is

U\cap\mathcal{M}=F^{-1}(\{0\})

for some open set $U\subset\mathcal{M}$ with $\bm{u}\in U$ , the tangent space to $\mathcal{M}$ at $\bm{u}$ can be calculated via differentiation:

T_{\bm{u}}\mathcal{M}=\operatorname{Ker}(DF_{\bm{u}}).

It is easily seen that the sphere has the defining equation $F(\bm{u})=\|\bm{u}\|_{2}^{2}-1$ . Show, using these facts, that the tangent space to the sphere at $\bm{u}$ is given by

T_{\bm{u}}\mathbb{S}^{d-1}=\{\bm{v}\in\mathbb{R}^{d}\mid\langle\bm{v},\bm{u}\rangle=0\},

and that the orthogonal projection onto this subspace is $\bm{P}_{\bm{u}}^{\perp}=\bm{I}-\bm{u}\bm{u}^{\top}$ .

The vector field

\mathrm{grad}\,f(\bm{u})=\bm{P}_{\bm{u}}^{\perp}\nabla f

(2.5.2)

is known as the Riemannian gradient of the function $f$ restricted to the sphere. The first order optimality conditions for the optimization problem (2.5.1) can be expressed in terms of the Riemann gradient:

\mathrm{grad}\,f(\bm{u})=\mathbf{0}.

Geometrically, this says that the Euclidean gradient of $f$ at $\bm{u}$ must be orthogonal to the tangent space to the sphere at $\bm{u}$ . Now suppose $\bm{v}\in\mathbb{R}^{d}$ is nonzero. Show that

\mathrm{proj}_{\mathbb{S}^{d-1}}(\bm{v})\doteq\min_{\|\bm{u}\|_{2}^{2}=1}\,\|\bm{u}-\bm{v}\|_{2}=\frac{\bm{v}}{\|\bm{v}\|_{2}},

using the first-order optimality conditions.

In optimization over $\mathbb{R}^{d}$ , one checks the second-order optimality conditions (to determine whether a critical point is a maximizer, a minimizer, or a saddle point) using the Hessian matrix $\nabla^{2}f(\bm{u})$ . Show, by differentiating the Riemann gradient $\mathrm{grad}\,f(\bm{u})$ for the sphere with respect to $\bm{u}$ as in the first part of this exercise, that the corresponding object for determining second-order optimality conditions for sphere-constrained optimization is the Riemannian Hessian, defined as

\mathrm{Hess}\,f(\bm{u})=\bm{P}_{\bm{u}}^{\perp}\left(\nabla^{2}f(\bm{u})-\langle\nabla f(\bm{u}),\bm{u}\rangle\bm{I}\right)\bm{P}_{\bm{u}}^{\perp}.

(2.5.3)

Exercise 2.7.

In this exercise, we sketch an argument referred to in the literature as a landscape analysis for the spherically-constrained population kurtosis maximization problem (2.2.24). We will show that when there is at least one independent component with positive kurtosis, its global maximizers indeed lead to the recovery of one column of the dictionary $\bm{U}$ . For simplicity, we will assume that $\mathop{\mathrm{kurt}}(z_{i})\neq 0$ for each $i=1,\dots,d$ .

Using the results of Part 1 of Exercise 2.6, show that the first-order optimality condition for (2.2.24) is

\left(\sum_{i=1}^{d}\mathop{\mathrm{kurt}}(z_{i})w_{i}^{4}\right)\bm{w}=\mathop{\mathrm{kurt}}(\bm{z})\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}\bm{w}^{\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}3},

(2.5.4)

where the kurtosis is calculated elementwise, $\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}$ denotes elementwise multiplication of vectors and $\bm{w}^{\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\odot$}}{\scalebox{0.8}{$\textstyle\odot$}}{\scalebox{0.8}{$\scriptstyle\odot$}}{\scalebox{0.8}{$\scriptscriptstyle\odot$}}$}}}3}$ denotes the elementwise cube of its argument.

Show that the vectors $\bm{w}$ with unit norm that also satisfy (2.5.4) all take the following form. Let $S^{+}=\{i\in[d]\mid\mathop{\mathrm{kurt}}(z_{i})>0\}$ , and $S^{-}=\{i\in[d]\mid\mathop{\mathrm{kurt}}(z_{i})<0\}$ . Let $S$ be a subset either of $S^{+}$ or $S^{-}$ . Then

\bm{w}_{S}=\sum_{i\in S}\pm\sqrt{\frac{1}{\mathop{\mathrm{kurt}}(z_{i})\sum_{j\in S}\frac{1}{\mathop{\mathrm{kurt}}(z_{j})}}}\bm{e}_{i}

(2.5.5)

satisfies (2.5.4), where $\bm{e}_{i}$ is the vector with a $1$ in the $i$ -th position and $0$ s elsewhere, and the $\pm$ sign denotes the choice of either a positive or negative sign.

3.

Assume that there is at least one $i$ such that $\mathop{\mathrm{kurt}}(z_{i})>0$ . Using the results of Part 2 of Exercise 2.6, show that the only local maxima of the objective of (2.2.24) are the signed one-sparse vectors $\pm\bm{e}_{i}$ with $i\in S^{+}$ . Conclude that the global maximizers of (2.2.24) are the signed one-sparse vectors corresponding to components with maximum kurtosis. (Hint: count the number of positive and negative eigenvalues of the Riemannian Hessian (2.5.3) at each critical point.)
4.

Now assume that $\mathop{\mathrm{kurt}}(z_{j})<0$ for every $j=1,\dots,d$ . This corresponds to an “over-deflated” instantiation of the kurtosis maximization problem. Using again the results of Part 2 of Exercise 2.6, show that the only local maxima of the objective of (2.2.24) are the signed dense vectors $\sum_{i=1}^{d}\pm\bm{e}_{i}$ . This shows that the optimization formulation (2.2.24) cannot be applied naively.

Exercise 2.8.

This exercise follows the structure and formalism introduced in Exercise 2.6, but applies it instead to the orthogonal group $\mathsf{O}(d)=\{\bm{U}\in\mathbb{R}^{d\times d}\mid\bm{U}^{\top}\bm{U}=\bm{I}\}$ . Consult the description of Exercise 2.6 for the necessary conceptual background; the formalism applies identically to the case where the ambient space is the set of $d\times d$ matrices as long as one recalls that the relevant inner product on matrices is $\langle\bm{X},\bm{Y}\rangle=\operatorname{tr}(\bm{X}^{\top}\bm{Y})$ . An excellent general reference for facts about optimization on the orthogonal group is Edelman, Arias, and Smith [EAS98].

Let $f:\mathbb{R}^{d\times d}\to\mathbb{R}$ be a given twice-continuously-differentiable objective function. Consider the orthogonally-constrained optimization problem

\max_{\bm{Q}^{\top}\bm{Q}=\bm{I}}\,f(\bm{Q}).

(2.5.6)

It is easily seen that the orthogonal group has the defining equation $F(\bm{Q})=\bm{Q}^{\top}\bm{Q}=\bm{I}$ . Show, using this fact, that the tangent space to the orthogonal group at $\bm{Q}$ is given by

T_{\bm{Q}}\mathsf{O}(d)=\{\bm{Q}\bm{\Omega}\in\mathbb{R}^{d\times d}\mid\bm{\Omega}^{\top}=-\bm{\Omega}\},

and that the orthogonal projection onto this subspace is

\mathcal{P}_{T_{\bm{Q}}\mathsf{O}(d)}(\bm{\Delta})=\bm{Q}\operatorname{skew}(\bm{Q}^{\top}\bm{\Delta}),

where $\operatorname{skew}(\bm{\Delta})=\tfrac{1}{2}(\bm{\Delta}-\bm{\Delta}^{\top})$ is the orthogonal projection onto the set of skew-symmetric matrices. The vector field

\mathrm{grad}\,f(\bm{Q})=\mathcal{P}_{T_{\bm{Q}}\mathsf{O}(d)}\left(\nabla f(\bm{Q})\right)

(2.5.7)

is known as the Riemannian gradient of the function $f$ restricted to the orthogonal group. The first order optimality conditions for the optimization problem (2.5.6) can be expressed in terms of the Riemann gradient:

\mathrm{grad}\,f(\bm{Q})=\mathbf{0}.

Show, by differentiating the Riemann gradient $\mathrm{grad}\,f(\bm{Q})$ for the orthogonal group with respect to $\bm{Q}$ as in the first part of this exercise, that the Riemannian Hessian is given by

\mathrm{Hess}\,f(\bm{Q})=\mathcal{P}_{T_{\bm{Q}}\mathsf{O}(d)}\left(\nabla^{2}f(\bm{Q})-\operatorname{sym}(\bm{Q}^{\top}\nabla f(\bm{Q}))\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\otimes$}}{\scalebox{0.8}{$\textstyle\otimes$}}{\scalebox{0.8}{$\scriptstyle\otimes$}}{\scalebox{0.8}{$\scriptscriptstyle\otimes$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\otimes$}}{\scalebox{0.8}{$\textstyle\otimes$}}{\scalebox{0.8}{$\scriptstyle\otimes$}}{\scalebox{0.8}{$\scriptscriptstyle\otimes$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\otimes$}}{\scalebox{0.8}{$\textstyle\otimes$}}{\scalebox{0.8}{$\scriptstyle\otimes$}}{\scalebox{0.8}{$\scriptscriptstyle\otimes$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\otimes$}}{\scalebox{0.8}{$\textstyle\otimes$}}{\scalebox{0.8}{$\scriptstyle\otimes$}}{\scalebox{0.8}{$\scriptscriptstyle\otimes$}}$}}}\bm{I}\right)\mathcal{P}_{T_{\bm{Q}}\mathsf{O}(d)},

(2.5.8)

where $\operatorname{sym}(\bm{\Delta})=\tfrac{1}{2}(\bm{\Delta}+\bm{\Delta}^{\top})$ denotes the orthogonal projection onto the set of symmetric matrices, and $\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\otimes$}}{\scalebox{0.8}{$\textstyle\otimes$}}{\scalebox{0.8}{$\scriptstyle\otimes$}}{\scalebox{0.8}{$\scriptscriptstyle\otimes$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\otimes$}}{\scalebox{0.8}{$\textstyle\otimes$}}{\scalebox{0.8}{$\scriptstyle\otimes$}}{\scalebox{0.8}{$\scriptscriptstyle\otimes$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\otimes$}}{\scalebox{0.8}{$\textstyle\otimes$}}{\scalebox{0.8}{$\scriptstyle\otimes$}}{\scalebox{0.8}{$\scriptscriptstyle\otimes$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\otimes$}}{\scalebox{0.8}{$\textstyle\otimes$}}{\scalebox{0.8}{$\scriptstyle\otimes$}}{\scalebox{0.8}{$\scriptscriptstyle\otimes$}}$}}}$ denotes the Kronecker product of matrices. Take care to interpret the operators appearing in the previous expression as linear transformations on ${d\times d}$ matrices, not as $d\times d$ matrices themselves. The second-order optimality conditions for the optimization problem (2.5.6) can be expressed in terms of the Riemann Hessian:

\mathrm{Hess}\,f(\bm{Q})\preceq\mathbf{0}.

For a minimization problem, the sign is reversed.

(Hint: The key is to manipulate one’s calculations to obtain the form (2.5.8), which is as compact as possible. To this end, make use of the following isomorphism of the Kronecker product: if $\bm{A}$ , $\bm{X}$ , and $\bm{B}$ are matrices of compatible sizes, then one has

(\bm{B}^{\top}\mathbin{\mathchoice{\raisebox{1.3pt}{$\displaystyle\mathchoice{\scalebox{0.8}{$\displaystyle\otimes$}}{\scalebox{0.8}{$\textstyle\otimes$}}{\scalebox{0.8}{$\scriptstyle\otimes$}}{\scalebox{0.8}{$\scriptscriptstyle\otimes$}}$}}{\raisebox{1.3pt}{$\mathchoice{\scalebox{0.8}{$\displaystyle\otimes$}}{\scalebox{0.8}{$\textstyle\otimes$}}{\scalebox{0.8}{$\scriptstyle\otimes$}}{\scalebox{0.8}{$\scriptscriptstyle\otimes$}}$}}{\raisebox{0.75pt}{$\scriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\otimes$}}{\scalebox{0.8}{$\textstyle\otimes$}}{\scalebox{0.8}{$\scriptstyle\otimes$}}{\scalebox{0.8}{$\scriptscriptstyle\otimes$}}$}}{\raisebox{0.6pt}{$\scriptscriptstyle\mathchoice{\scalebox{0.8}{$\displaystyle\otimes$}}{\scalebox{0.8}{$\textstyle\otimes$}}{\scalebox{0.8}{$\scriptstyle\otimes$}}{\scalebox{0.8}{$\scriptscriptstyle\otimes$}}$}}}\bm{A})\operatorname{\mathrm{vec}}(\bm{X})=\operatorname{\mathrm{vec}}(\bm{A}\bm{X}\bm{B}),

where $\operatorname{\mathrm{vec}}(\bm{X})$ denotes the “left-to-right” stacking of the columns of the matrix argument into a vector. We use this isomorphism in (2.5.8) in order to define the Kronecker product of two matrices as an operator on matrices in a canonical way.)

Now suppose $\bm{X}\in\mathbb{R}^{d\times d}$ is full-rank. In this and the next part of the exercise, we consider the projection onto the orthogonal group of $\bm{X}$ :

\mathrm{proj}_{\mathsf{O}(d)}(\bm{X})\doteq\min_{\bm{Q}\in\mathsf{O}(d)}\,\|\bm{Q}-\bm{X}\|_{F}^{2}.

(2.5.9)

We will prove that the solution to this problem is given by

\mathrm{proj}_{\mathsf{O}(d)}(\bm{X})=\bm{U}\bm{V}^{\top},

where $\bm{X}=\bm{U}\bm{S}\bm{V}^{\top}$ is a singular value decomposition of $\bm{X}$ .

(a)

Using the first and second-order optimality conditions, show that every local minimizer $\bm{Q}$ of (2.5.9) satisfies

	$\displaystyle\left(\bm{Q}^{\top}\bm{X}\right)^{\top}$	$\displaystyle=\bm{Q}^{\top}\bm{X},$
	$\displaystyle\bm{Q}^{\top}\bm{X}$	$\displaystyle\succeq\mathbf{0}.$

(Hint: use linearity of the Kronecker product in either of its two arguments when the other is fixed.)

(b)

Using these conditions, argue that at every local minimizer $\bm{Q}$ of (2.5.9), one has $\bm{Q}^{\top}\bm{X}=(\bm{X}^{\top}\bm{X})^{1/2}$ . (Hint: Use the fact from linear algebra that if $\bm{S}\succeq\mathbf{0}$ is a symmetric positive semidefinite matrix, then $(\bm{S}^{\top}\bm{S})^{1/2}=\bm{S}$ .)

(c)

Using the singular value decomposition $\bm{X}=\bm{U}\bm{S}\bm{V}^{\top}$ , conclude that

\bm{U}\bm{V}^{\top}=\mathrm{proj}_{\mathsf{O}(d)}(\bm{X}).

	$\displaystyle\min_{\{\tilde{\bm{z}}_{i}\}_{i=1}^{N}}\frac{1}{N}\sum_{i=1}^{N}\\|\bm{x}_{i}-\tilde{\bm{U}}\tilde{\bm{z}}_{i}\\|_{2}^{2}$	$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\min_{\tilde{\bm{z}}_{i}}\\|\bm{x}_{i}-\tilde{\bm{U}}\tilde{\bm{z}}_{i}\\|_{2}^{2}$		(2.1.3)
		$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\\|\bm{x}_{i}-\tilde{\bm{U}}\tilde{\bm{U}}^{\top}\bm{x}_{i}\\|_{2}^{2}.$		(2.1.4)