Canonical Correlation Analysis: review

Canonical Correlation Analysis: review

\address[Anna Bykhovskaya]{Duke University} \email{[email removed]}

\address[Vadim Gorin]{University of California at Berkeley} \email{[email removed]}

Preface

For over a century canonical correlations, variables, and related concepts have been studied across various fields, with contributions dating back to jordan1875essai and harold1936relations. This text surveys the evolution of canonical correlation analysis, a fundamental statistical tool, beginning with its foundational theorems and progressing to recent developments and open research problems. Along the way we introduce and review methods, notions, and fundamental concepts from linear algebra, random matrix theory, and high-dimensional statistics, placing particular emphasis on rigorous mathematical treatment.

The survey is intended for technically proficient graduate students and other researchers with an interest in this area. The content is organized into five chapters, supplemented by six sets of exercises found in Chapter (ref). These exercises introduce additional material, reinforce key concepts, and serve to bridge ideas across chapters. We recommend the following sequence: first, solve Problem Set 0, then proceed with Chapter 1, solve Problem Set 1, and so on through the text.

The review grew out of lectures delivered by Gorin at the 2024 Random Matrix Theory Summer School at the University of Michigan, and from the topics class taught by Bykhovskaya at Duke in the Fall 2024. We are grateful to the organizers and participants of these events for their valuable contributions. Gorin's work was partially supported by NSF grant DMS - 2246449.

Introduction and basic definitions

Opening remarks

Broadly speaking, Canonical Correlation Analysis (CCA) is a method for identifying common factors between two multi-dimensional objects. It can be seen as a companion to the Principal Component Analysis (PCA):

Empirically, CCA is widely applied in testing and inferring relationships between data sets: given two data sets $A$ and $B$, one can test for independence by examining canonical correlations between the row spaces of their respective matrices. If independence is rejected, CCA can then be used to identify the most interdependent components, or linear combinations, of $A$ and $B$. Additionally, applying CCA to spaces derived from the same data sets but through slightly varied procedures can uncover structural and temporal properties within the data.

Examples include: in genomics, finding shared patterns across assays from the same individuals; in neuroscience, linking brain activity with behavioral data; in ecology, correlating species characteristics with habitats; in econometrics, selecting common factors and performing cointegration analysis; and in finance, optimizing portfolio allocations. We will not cover all possible applications but provide references as starting points: gittins1985canonical,johansen_book, simon1998assessing, sherry2005conducting,witten2009extensions,breitung2013canonical,andreou2019inference,wang2020finding,zhuang2020technical, choi2021canonical,franchi2023estimating,firoozye2023canonical.

Three points of view on CCA

CCA can be approached from the probabilistic, statistical, and geometrical perspectives. The probabilistic (or population) framework deals with two families of random variables and measures their dependence, while the statistics (or sample) framework assumes that instead of observing the actual distributions (i.e., knowing the mean, variance, and so on) one observes realizations or samples from the distributions. This reflects real-world data scenarios, such as when researchers do not know the true distribution of stock returns, but have access to daily observations. Finally, the geometrical framework unifies both the probabilistic and statistical settings by interpreting random variables or their samples as vectors in an appropriately defined vector space.

Before introducing general definitions of CCA, it is helpful to first examine a one-dimensional case, where everything reduces to well-known notions of (sample) correlations and angles.

Depending on the reader's preference, one may choose to work with one of the above frameworks. Yet, it is helpful to remember that (ref), (ref), and (ref) are three faces of the same mathematical object: both (ref) and (ref) can be turned into the form (ref) by properly specifying the vector space and the scalar product.

Canonical correlations provide a multi-dimensional extension of (ref), (ref), and (ref). Instead of two random variables we consider two groups of random variables; rather than using two vectors of data, we work with two rectangular matrices of data; and instead of lines we examine linear subspaces.

The concepts behind CCA were originally introduced in jordan1875essai,harold1936relations and remain widely used today. While we will cover all the key definitions, for additional reading, we recommend statistics textbooks such as thompson1984canonical,gittins1985canonical,anderson1958introduction, and muirhead2009aspects.

Definitions: canonical correlations and variables

In this subsection we provide formal definitions of canonical correlations and canonical variables and introduce different algorithms, which can be used to calculate them. We start with geometric framework via Theorem (ref) and then explain how to specialize it to probabilistic and statistics frameworks. We then present two alternative ways to calculate canonical correlations and canonical variables, which are commonly used perspectives in data science when thinking about CCA. The first approach can be viewed as a form of dimension reduction, where we iteratively capture as much of the correlation structure as possible by working with linear combinations of the data. It is also used to construct the bases of Theorem (ref). The second approach relates CCA to the eigenvalues and eigenvectors of specific matrices. The proof of Theorem (ref) as well as some additional results are presented in Section (ref).

Two bases

The numbers $c_1\ge c_2 \ge \dots \ge c_{K}$ are called canonical correlations between subspaces $\mathbf U$ and $\mathbf V$; they are also cosines of the canonical angles between the subspaces. The vectors $\mathbf u_i$, $1\le i \le {K}$, and $\mathbf v_j$, $1\le j\le {M}$, are called canonical variables; they split into ${K}$ pairs $(\mathbf u_i,\mathbf v_i)$ and $M-K$ singletons $\mathbf v_j$, $j\ge K$.

Theorem (ref) constructs an orthonormal basis in the ambient space, which encompasses $\mathbf U$ and $\mathbf V$. It can be treated as an analogue of the diagonalization procedure for symmetric matrices (i.e., the theorem which says that any symmetric matrix can be made diagonal by a choice of an appropriate basis of the space). As will be seen in Corollary (ref), canonical correlations $c_i$ are invariants of spaces $\mathbf U$ and $\mathbf V$ under orthogonal transformations of the space: they are present, no matter what orthogonal basis one chooses. The point of view of $c_i$ being invariants is explored, e.g., in halmos1969two and bottcher2010gentle.

Let us adapt Theorem (ref) to probability and statistics frameworks.

Theorem (ref) in the probability framework. Suppose we are given a collection of mean $0$ random variables $\{\mathsf x^i\}_{i=1}^K$ and another collection $\{\mathsf y^j\}_{j=1}^M$. Then we take as $\mathbf U$ the space of all linear combinations $\alpha_1 \mathsf x^1+\alpha_2 \mathsf x^2+\dots+\alpha_K\mathsf x^K$ with deterministic real coefficients $\alpha_i$. Similarly, $\mathbf V$ is the space of all linear combinations of $\mathsf y^j$ with deterministic real coefficients. The scalar product is given by expectation: $\langle \mathsf u,\mathsf v\rangle=\mathbb E \mathsf u \mathsf v$ for mean $0$ random variables $\mathsf u$ and $\mathsf v$. In this situation, the canonical variables are also random variables: $\mathbf u_i$ are linear combinations of $\mathsf x^1,\dots,\mathsf x^K$ and $\mathbf v_j$ are linear combinations of $\mathsf y^1,\dots,\mathsf y^M$ with deterministic coefficients. On the other hand, the canonical correlations $c_1,\dots,c_K$ are deterministic --- they describe the correlation structure between families $\{\mathsf x^i\}_{i=1}^K$ and $\{\mathsf y^j\}_{j=1}^M$. In particular, if $K=M=1$, then we are back to (ref) with $\rho=c_1$: $$ \mathbf u_1=\pm\frac{\mathsf x^1}{\sqrt{\mathbb E (\mathsf x^1)^2}}, \qquad \mathbf v_1=\pm\frac{\mathsf y^1}{\sqrt{\mathbb E (\mathsf y^1)^2}}, \qquad c_1= \frac{\left|\mathbb E \mathsf x^1\mathsf y^1\right|}{\sqrt{\mathbb E (\mathsf x^1)^2 \mathbb E (\mathsf y^1)^2}}, $$\ where the signs are chosen so that $\langle \mathbf u_1,\mathbf v_1\rangle \ge 0$.

For general $M\ge K$ Theorem (ref) gives the simplest possible form to which a correlation structure between two families of random variables $\{\mathsf x^i\}_{i=1}^K$ and $\{\mathsf y^j\}_{j=1}^M$ can be reduced by linear transformations. A useful analogy can be drawn with the theorem for Gaussian vectors, which states that any zero-mean Gaussian vector can be obtained by applying a linear transformation to a Gaussian vector with identity covariance matrix. Similarly, Theorem (ref) asserts that the $(K+M)\times (K+M)$ covariance matrix of the vector $(\mathsf x^1,\dots, \mathsf x^K, \mathsf y^1, \dots, \mathsf y^M)$ can be reduced by two linear transformations (one applied to $\{\mathsf x^i\}_{i=1}^K$ and another to $\{\mathsf y^j\}_{j=1}^M$) to a simple block form:

{ $$\left(

\right), $$} where the block sizes are $K\times K$, $K\times M$ in the top and $M\times K$, $M\times M$ in the bottom.

Theorem (ref) in the statistics framework. Next, suppose we are given two rectangular matrices of data: $K\times S$ matrix $X=[X_{ij}]$ and $M\times S$ matrix $Y=[Y_{ij}]$. Then we take as $\mathbf U\subset \mathbb R^S$ the space of all linear combinations of $K$ rows of $X$ and $\mathbf V\subset\mathbb R^S$ is the space of all linear combinations of $M$ rows of $Y$. The scalar product is the standard scalar product in $\mathbb R^S$: $\langle (x_1,\dots,x_S), (y_1,\dots,y_S)\rangle = \sum_{i=1}^S x_i y_i$. In this situation the sample\footnote{In this context we add the term “sample” to indicate the statistical, or sampled, nature of the data.} canonical variables are $S$-dimensional vectors: $\mathbf u_i$ are linear combinations of rows of $X$ and $\mathbf v_j$ are linear combinations of rows of $Y$. The sample canonical correlations $c_1,\dots,c_K$ are numbers describing the cross-correlation structure between the rows of $X$ and $Y$. In particular, if $K=M=1$, then we are back to (ref) with $\hat \rho=c_1$: $$ \mathbf u_1=\pm \frac{ (X_{11},X_{12},\dots,X_{1S})}{\sqrt{\sum\limits_{i=1}^S X_{1i}^2}},\qquad \mathbf v_1=\pm \frac{(Y_{11},Y_{12},\dots,Y_{1S})}{\sqrt{\sum\limits_{i=1}^S Y_{1i}^2}},\qquad c_1=\frac{\left|\sum\limits_{i=1}^S X_{1i}Y_{1i} \right| }{\sqrt{\sum\limits_{i=1}^S X_{1i}^2\sum\limits_{i=1}^S Y_{1i}^2}}. $$

Maximization problem

The first approach for finding the canonical correlations and variables of Theorem (ref) is to consider a function $f(\mathbf u,\mathbf v)=\langle \mathbf u, \mathbf v\rangle$, in which $\mathbf u$ varies over all unit vectors in $\mathbf U$ and $\mathbf v$ varies over all unit vectors in $\mathbf V$. This function represents (sample) correlation\footnote{Since correlation between two variables is unaffected by rescaling of the variables, we can restrict the space to be the unit sphere.} between $\mathbf u$ and $\mathbf v$. Then $c_1$ is the maximum of $f$, which is achieved at $(\mathbf u_1,\mathbf v_1)$. After $(\mathbf u_1,\mathbf v_1,c_1)$ is found, we consider restriction of the same function $f(\mathbf u,\mathbf v)$ on the vectors $\mathbf u \in \mathbf U$ orthogonal to $\mathbf u_1$ and vectors $\mathbf v\in \mathbf V$ orthogonal to $\mathbf v_1$. Maximizing the restricted function, we get $c_2$, which is achieved at $(\mathbf u_2,\mathbf v_2)$. Next, we restrict $f$ to $\mathbf u\in \mathbf U$ orthogonal to both $\mathbf u_1$ and $\mathbf u_2$ and to $\mathbf v\in \mathbf V$ orthogonal to both $\mathbf v_1$ and $\mathbf v_2$, etc. Formally, step $i$ of the sequential maximization problem can be written as

Once all $(\mathbf u_i,\mathbf v_i,c_i)$, $1\le i \le {K}$, are identified, the remaining vectors $\mathbf v_j$, $j>{K}$, are an arbitrary orthogonal basis in the part of $\mathbf V$ orthogonal to $\mathbf v_1,\dots,\mathbf v_{K}$. This approach reduces the task of finding canonical correlations and variables to a series of iterative maximization problems, which can be efficiently solved using numerous numerical algorithms and also admit various generalizations.

If we wish to solve maximization problems (ref) simultaneously for all $1\le i \le K$, then we can view $(\mathbf u_i,\mathbf v_i,c_i)$, $1\le i \le {K}$, as critical points and corresponding values of $f$. Given the constraints $\langle\mathbf u,\mathbf u\rangle=\langle\mathbf v,\mathbf v\rangle=1$, these critical points can be found through the Lagrangian function: $$ \mathcal L(\mathbf u,\mathbf v)=\langle \mathbf u, \mathbf v\rangle+ a \langle\mathbf u,\mathbf u\rangle + b \langle\mathbf v,\mathbf v\rangle, \qquad \mathbf u\in\mathbf U, \quad \mathbf v\in \mathbf V. $$ where $a$ and $b$ are Lagrange multipliers. Expanding $\mathbf u$ and $\mathbf v$ in arbitrary bases of $\mathbf U$ and $\mathbf V$, respectively (in the probabilistic setting it is natural to expand in $\{\mathsf x^i\}_{i=1}^K$ and $\{\mathsf y^j\}_{j=1}^M$; while in the statistical setting it is natural to expand in rows of $X$ and $Y$), we differentiate $\mathcal{L}(\mathbf u, \mathbf v)$ with respect to all coefficients of the expansions and set the partial derivatives to zero. This yields a system of linear equations, whose solutions provide the pairs $(\pm \mathbf u_i, \pm \mathbf v_i)$.

Eigenvalues and eigenvectors

Here is an alternative approach for finding the canonical correlations and variables of Theorem (ref). For a subspace $\mathbf Q\subset \mathbf W$ let us denote through $P_{\mathbf Q}$ the orthogonal projector onto $\mathbf Q$.

In order to use Proposition (ref) as a numerical algorithm, one needs to calculate the matrices of the projectors. If we use the statistics point of view where $\mathbf U$ and $\mathbf V$ are rows of rectangular matrices $X$ and $Y$, then

where $(\cdot)^{\mathsf T}$ here and below denotes matrix transposition. Note that if the rows of $X$ and $Y$ are orthogonal, then $X X^\mathsf T$ and $Y Y^\mathsf T$ are identity matrices, whose inversion is straightforward. This feature is sometimes used for developing more efficient numerical methods.

The product of projectors $P_{\mathbf U} P_{\mathbf V}$ in (ref) is a $\dim(\mathbf W)\times \dim(\mathbf W)$ dimensional matrix, which can be hard to operate with if $\dim(\mathbf W)$ is large. In particular, if we are in the probability setting with random variables $\{\mathsf x^i\}_{i=1}^K$ and $\{\mathsf y^j\}_{j=1}^M$, it is natural to take as $\mathbf W$ the space of all square-integrable random variables (on some probability space), which makes $\mathbf W$ infinite-dimensional and, therefore, directly evaluating $P_{\mathbf U} P_{\mathbf V}$ becomes computationally challenging. However, as formally proved later in Proposition (ref), we can simplify the algorithm and deal with $K\times K$ and $M\times M$ matrices instead:

$\bullet$ In the statistics framework the squared sample canonical correlations can be equivalently found as eigenvalues of either of the matrices

where we note that the first matrix is $K\times K$ and the second one is $M\times M$. If $\boldsymbol\alpha$ is an eigenvector of the first matrix and $\boldsymbol\beta$ is an eigenvector of the second matrix with the same eigenvalue $c$, then, up to normalization of vectors, the triplets $(\mathbf U^\mathsf T\boldsymbol\alpha, \mathbf V^\mathsf T\boldsymbol\beta, c)$ represent sample canonical variables and correlations between $\mathbf U$ and $\mathbf V$. Sometimes, $\boldsymbol\alpha$ and $\boldsymbol\beta$ are referred to as canonical vectors, which should not be confused with canonical variables $\mathbf U^\mathsf T\boldsymbol\alpha$ and $\mathbf V^\mathsf T\boldsymbol\beta$.

$\bullet$ Similarly, in the probability framework the squared (populational) canonical correlations can be found as eigenvalues of either of the following matrices, in which we used the vector notations $\mathbf u=(\mathsf x^1,\mathsf x^2,\dots,\mathsf x^K)^{\mathsf T}$, $\mathbf v=(\mathsf y^1,\mathsf y^2,\dots,\mathsf y^M)^{\mathsf T}$,

If $\boldsymbol\alpha$ is an eigenvector of the first matrix and $\boldsymbol\beta$ is an eigenvector of the second matrix with the same eigenvalue $c$, then, up to normalization of vectors, the triplets $(\mathbf u^\mathsf T\boldsymbol\alpha, \mathbf v^\mathsf T\boldsymbol\beta, c)$ represent canonical variables and correlations between the families $\{\mathsf x^i\}_{i=1}^K$ and $\{\mathsf y^j\}_{j=1}^M$.

Proof of Theorem (ref) and corollaries

This section presents proofs of all the previously mentioned facts, along with some additional results. We start by showing that the iterative maximization algorithm (ref) produces the desired bases $\mathbf u_1, \mathbf u_2, \dots, \mathbf u_{K}$ and $\mathbf v_1,\dots, \mathbf v_{M}$.

Let $U^1=\{\mathbf u\in\mathbf U\mid \langle \mathbf u,\mathbf u\rangle =1\}$, $V^1=\{\mathbf v\in\mathbf V\mid \langle \mathbf v,\mathbf v\rangle=1\}$ be two unit spheres in $K$ and $M$ dimensional spaces, respectively. Consider the function $f(\mathbf u,\mathbf v)=\langle \mathbf u,\mathbf v\rangle$ restricted to these spheres. If $f$ is identical zero, then we are done: the canonical correlations $c_i$ are $0$ and we can choose an arbitrary basis of $\mathbf U$ as $\mathbf u_1,\dots,\mathbf u_K$ and an arbitrary basis of $\mathbf V$ ad $\mathbf v_1,\dots,\mathbf v_M$.

If $f(\mathbf u,\mathbf v)$ is not identical zero, then $f(-\mathbf u,\mathbf v)=-f(\mathbf u,\mathbf v)$ implies that $f$ takes positive values at some points. Since $f$ is a continuous function and spheres are compact, $f$ achieves its maximum at a point. Let $0<c_1\le 1$ denote the maximal value of $f$, and let $(\mathbf u_1,\mathbf v_1)$ denote the point where it is achieved. By definition $\langle \mathbf u_1,\mathbf u_1\rangle=\langle \mathbf v_1,\mathbf v_1\rangle=1$ and $\langle \mathbf u_1,\mathbf v_1\rangle=c_1$, fitting into (ref).

Further, set $U^2=\{\mathbf u\in U^1\mid \langle \mathbf u,\mathbf u_1\rangle =0\}$, $V^2=\{v\in V^1\mid \langle \mathbf v,\mathbf v_1\rangle=0\}$. These are $K-1$ and $M-1$ dimensional spheres, respectively, with an important property:

Next, we consider the restriction of $f$ onto $U^2\times V^2$. If it is identical zero, then $c_2=c_3=\dots=0$, and we can complement $\mathbf u_1$ to an orthonormal basis of $\mathbf U$ and complement $\mathbf v_1$ to an orthonormal basis of $\mathbf V$, getting the desired (ref).

Otherwise, let $c_2$ denote the maximum of $f$ on $U^2\times V^2$, achieved at $(\mathbf u_2,\mathbf v_2)$. By definitions and Lemma (ref), we have: $$ \langle \mathbf u_2,\mathbf u_2\rangle =\langle \mathbf v_2,\mathbf v_2\rangle=1,\quad \langle \mathbf u_2,\mathbf u_1\rangle=\langle \mathbf u_2,\mathbf v_1\rangle=\langle \mathbf v_2,\mathbf u_1\rangle= \langle \mathbf v_2,\mathbf v_1\rangle =0,\quad \langle \mathbf u_2,\mathbf v_2\rangle = c_2, $$ fitting into (ref). We can continue in the same way by defining $K-2$ and $M-2$ dimensional spheres $ U^3=\{\mathbf u\in U^2\mid \langle \mathbf u,\mathbf u_2\rangle =0\}$, $V^3=\{v\in V^2\mid \langle \mathbf v,\mathbf v_2\rangle=0\}$. Then we use the same argument as in Lemma (ref) to show that $U^3$ is orthogonal to $\mathbf v_1$ and $\mathbf v_2$ and $V^3$ is orthogonal to $\mathbf u_1$ and $\mathbf u_2$, maximize $f$ over $U^3\times V^3$, etc.

Repeating the argument $K$ times, we exhaust $\mathbf U$ and construct its orthonormal basis $\mathbf u_1,\dots,\mathbf u_K$ and corresponding partial basis $\mathbf v_1,\dots,\mathbf v_K$ of $\mathbf V$. Yet again by a version of Lemma (ref), the complement of $\mathbf v_1,\dots,\mathbf v_K$ in $\mathbf V$ is orthogonal to all of $\mathbf U$. Hence, we can extend $\mathbf v_1,\dots,\mathbf v_K$ to an orthonormal basis of $\mathbf V$ in an arbitrary way, and the scalar products (ref) will hold. This finishes the proof of Theorem (ref).

Let us now explain the validity of the second algorithm outlined in Section (ref).

We end this section with a discussion of uniqueness of the bases satisfying (ref) in Theorem (ref). If all $c_i$ are distinct, then Proposition (ref) implies that the only freedom in choosing $\mathbf u_1,\dots,\mathbf u_K$ and $\mathbf v_1,\dots,\mathbf v_K$ is in simultaneous multiplication of $\mathbf u_i$ and $\mathbf v_i$ by $-1$. For $\mathbf v_{K+1},\dots,\mathbf v_M$ we always have more options: any $M-K$ orthonormal vectors complementing $\mathbf v_1,\dots,\mathbf v_K$ to a basis of $\mathbf V$ would work.

If some $c_i$ coincide, then the situation is more delicate: eigenvectors in Proposition (ref) are no longer unique, only the linear spaces spanned by eigenvectors with the same eigenvalue are uniquely determined. In particular, imagine that $c_1$ has multiplicity two, so that for four unit vectors $\mathbf u_1,\mathbf u_2\in \mathbf U$ ,$\mathbf v_1,\mathbf v_2\in \mathbf V$ we have: $$ \langle \mathbf u_1,\mathbf v_1\rangle = \langle \mathbf u_2,\mathbf v_2\rangle = c_1, \qquad \langle \mathbf u_1,\mathbf u_2\rangle = \langle \mathbf u_1,\mathbf v_2\rangle = \langle \mathbf u_2,\mathbf v_1\rangle = \langle \mathbf v_1,\mathbf v_2\rangle = 0. $$ Then for any $-1\le \alpha\le 1$, we also have $$ \langle \alpha \mathbf u_1 + \sqrt{1-\alpha^2} \mathbf u_2, \alpha \mathbf v_1 + \sqrt{1-\alpha^2} \mathbf v_2 \rangle = \alpha^2 c_1 +(1-\alpha^2)c_1=c_1. $$ Hence the vectors $(\alpha \mathbf u_1 + \sqrt{1-\alpha^2} \mathbf u_2, \alpha \mathbf v_1 + \sqrt{1-\alpha^2} \mathbf v_2)$ also deliver the maximum of the function $f(\mathbf u,\mathbf v)$ and (by the algorithm in subsection (ref)) can be included as the first vectors in a system satisfying (ref).

Two extreme examples are when the subspaces $\mathbf U$ and $\mathbf V$ are either orthogonal or coincide. In both situations all $c_i$ are equal.

Questions and applications

There are many open questions that continue to drive research and applications of CCA. On the theoretical side, CCA-related questions in probability and statistics often center on the distributional properties -- both exact and asymptotic -- of canonical correlations and variables. In this review we focus on some of these topics. Specifically, we highlight that:

Returning to empirical applications, these two features allow us to test mathematical theorems about CCA on real-world data sets, yielding several insights. First, we often observe a remarkable alignment between theoretical predictions and empirical results in various data sets, highlighting the applicability of our modeling approach to real-world data. Second, we can leverage these theorems to derive structural insights from the data sets.

As a preview, we present three plots generated from financial data sets. Figure (ref) shows the histogram of sample canonical correlations for weekly data from S&P 100 stocks over a ten-year period, correlating the logarithms of stock prices with their time increments. Figure (ref) provides a similar plot but for daily data on 25 cryptocurrencies over two years. Figure (ref) presents squared sample canonical correlations for weekly returns over ten years for two groups of 80 stocks: largest “cyclical” versus largest “non-cyclical (defensive)” stocks. In each case the empirical histograms are accompanied by (closely matching) theoretical curves, that represent the Wachter distribution. The relevance of the Wachter distribution to CCA and these data sets in particular is discussed in Chapters 3 and 5. The canonical correlations appearing to the right of the support of Wachter distribution (the flat region of the theoretical curves) can be interpreted as signals in the data, which we elaborate on in Chapter 4. For further details on these figures and the corresponding theoretical results we refer to BG1, BG2, and BG3, respectively.

IID setting

CCA and Maximum Likelihood Estimation

This chapter develops a connection between statistical and probability frameworks. We start from the population model with two random vectors: $\mathbf u=(u^1,\dots,u^K)^{\mathsf T}$ and $\mathbf v=(v^1,\dots,v^M)^\mathsf T$, $M\ge K$. We do not have direct access to these vectors and we do not know anything about them, except that they have mean $0$. Instead, we are observing $S$ i.i.d.\ samples of these vectors. Formally, this means that we have two matrices: $\mathbf U$ of size $K\times S$ and $\mathbf V$ of size $M\times S$. The combined $(K+M)\times S$ matrix $\mathbf W=

$ has i.i.d.\ columns, with each column following the distribution of $\mathbf w=

$.

Let us clarify the meaning of the words “correlation structure”. There are three covariance matrices in play: $K\times K$ matrix of auto-covariances of $\mathbf u$, $\Lambda_{uu}[i,j]=\mathbb E u^i u^j$, and $M\times M$ matrix of auto-covariances of $\mathbf v$, $\Lambda_{vv}[i,j]=\mathbb E v^i v^j$, and $K\times M$ matrix of cross-covariances $\Lambda_{uv}[i,j]=\mathbb E u^i v^j$; we also denote $\Lambda_{vu}=\Lambda_{uv}^\mathsf T$. Theorem (ref) implies that there exists $K\times K$ matrix $A$ and $M\times M$ matrix $B$, such that

where $I_K$ and $I_M$ are identity matrices of $K\times K$ and $M\times M$ dimensions, respectively, and $\mathrm{diag}(c_1,\dots,c_K)$ is a rectangular matrix with $c_1,c_2,\dots,c_K$ on the main diagonal and $0$ everywhere else. Here $1\ge c_1\ge c_2\ge\dots\ge c_K\ge 0$ are canonical correlations between the spaces spanned by coordinates of $\mathbf u$ and coordinates of $\mathbf v$, $A \mathbf u$ is the random vector of $K$ canonical variables $(\mathbf u_1,\dots,\mathbf u_K)^{\mathsf T}$ and $B\mathbf v$ is the vector of $M$ canonical variables $(\mathbf v_1,\dots,\mathbf v_M)^{\mathsf T}$.

From (ref) we conclude that knowing the covariance matrices $\Lambda_{uu}$, $\Lambda_{vv}$, $\Lambda_{uv}$ is equivalent to knowing the canonical correlations and variables between the coordinates of $\mathbf u$ and $\mathbf v$. The advantage of (ref) is its clear separation of the internal correlation structures of $\mathbf u$ and $\mathbf v$, captured by $A$ and $B$, from the strength of the correlations between $\mathbf u$ and $\mathbf v$, represented $c_1,\dots,c_K$. Hence, our task can be reformulated as:

There are various statistical methods for estimating parameters, with one of the most classical being the Maximum Likelihood Estimator (MLE). The MLE procedure begins by assuming that the random variables of interest are sampled from a specific, predefined family of probability distributions (population model), that depend on our parameter of interest. Next, we compute the likelihood function -- a function of both the parameters and the observed samples -- that quantifies the probability, or probability density, of observing the given samples. The parameters that maximize this function, treating the observed samples as fixed, provide the desired estimator. The estimator is then treated as a function of the observed sample.

Theorem (ref) says that sample covariances are MLE for their population counterparts. Thus, in view of (ref), the same should hold for canonical correlations and variables and Task (ref) can be approached by working with the sample CCA:\footnote{For a different interpretation of canonical correlations as MLE estimators (for another probability model) see bach2005probabilistic.}

Several remarks are in order. First, the canonical variables for $\mathbf u$ and $\mathbf v$ are random variables (mathematically, these are functions on an abstract probability space), so we cannot estimate them directly, instead we are estimating the coefficients $A_{ij}$ and $B_{ij}$ of their decomposition in coordinates $u^1,\dots,u^K$ and $v^1,\dots,v^M$. On the other hand, canonical correlations are deterministic numbers, so they are being estimated directly in the third line. Second, the eigenvectors in the first two lines of (ref) are defined only up to multiplication by an arbitrary real number. If we want to have some normalization, we need to additionally impose it. Note that in Theorem (ref) there was also freedom in multiplying $\mathbf u_i$, $\mathbf v_i$ by $-1$, so one natural approach is to avoid normalizations at all, thinking about $\mathbf u_i$ and $\mathbf v_i$ as directions (or one-dimensional lines), rather than vectors.

The preceding statistical reasoning suggests to use Procedure (ref), however, we do not yet have any guarantees on the quality of estimation. Is there some sense, in which the estimates $\hat A, \hat B$, and $\hat c_1,\dots,\hat c_K$ are close to the true values of $A$, $B$ and $c_1,\dots,c_K$? We start addressing this question below and continue in the next two chapters. We will see that the answer crucially depends on how large $K$, $M$, and $S$ are.

Two independent data sets

Our next task is to understand the probability distribution of (ref), given some assumptions on vectors $\mathbf u$ and $\mathbf v$. The simplest possible setting is to assume that $\mathbf u$ and $\mathbf v$ are independent mean $0$ Gaussian vectors with arbitrary non-degenerate covariance matrices $\Lambda_{uu}$ and $\Lambda_{vv}$. That is, $\Lambda_{uv}$ is a zero matrix.

From the applied point of view, this setting can be used to construct tests for independence of two data sets. From the theoretical point of view, this setting reveals fruitful connections to classical ensembles of the random matrix theory and serves as a foundation for further studies of more complicated cases.

Jacobi ensemble

Miraculously, in the i.i.d.\ Gaussian setting for $\mathbf U$ and $\mathbf V$, the distribution of the sample canonical correlations is explicit, as was first noticed by hsu1939distribution.

For $N=1$, the matrix reduces to its single matrix element $\mathcal M_{11}\in (0,1)$ and Definition (ref) becomes the density of Beta distribution: $$ \frac{1}{\mathfrak B(p,q)} x^{p-1} (1-x)^{q-1},\quad 0<x<1,\qquad\qquad \mathfrak B(p,q)=\int_0^1 x^{p-1}(1-x)^{q-1}\mathrm d x=\frac{\Gamma(p)\Gamma(q)}{\Gamma(p+q)}. $$

Both (ref) and (ref) can be viewed as multidimensional generalizations of the Beta distribution. The two generalizations are closely related, as we explain in Theorem (ref). Their connection to the canonical correlation analysis is stated in Theorem (ref).

In $K=M=1$ case Theorem (ref) states that the sample correlation coefficient between two Gaussian vectors has a Beta distribution. The $K\le M$ restriction is just for convenience and can be removed by symmetry between $\mathbf U$ and $\mathbf V$. The $K+M\le S$ restriction is more interesting: note that if $K+M>S$, then any $K$--dimensional linear subspace in $S$--dimensional space will have a $(K+M-S)$--dimensional intersection with any $M$--dimensional linear subspace in $S$--dimensional space. Hence, in this situation there are $K+M-S$ deterministic canonical correlations equal to $1$. One can still analyze the joint distribution of the remaining correlations (it will be again given by an instance of the Jacobi ensemble with appropriate parameters), yet we will not address it here: from the point of view of statistical applications the case of large $S$ is the most natural one, as we will see later.

Theorem (ref) can be generalized to cases in which the columns of $\mathbf U$ and $\mathbf V$ are i.i.d.\ samples of Gaussian $\mathbf u$ and $\mathbf v$, where $\mathbf u$ and $\mathbf v$ allowed to be correlated, see constantine1963some. The distribution becomes much more complicated than the Jacobi ensemble --- it involves so-called matrix hypergeometric functions.

The proofs of Theorems (ref) and (ref) rely on the techniques for computing matrix integrals and Jacobians for transformations between various matrices and their eigenvalues. Three encyclopedic sources for such techniques are Hua, muirhead2009aspects, and Forrester10. We present the general schemes of the proofs in the next section omitting some of the details.

Scheme of the proof of Theorem (ref)

This is a very general computation which works for many ensembles of random matrices. The key point is to figure out the Jacobian of the map between a symmetric matrix and its eigenvalues.

Consider the map $$ \Phi: (O; \lambda_1,\dots,\lambda_N) \mapsto \mathbf X = O\, \mathrm{diag}(\lambda_1,\dots,\lambda_N)\, O^{\mathsf T}, $$ which takes a real orthogonal matrix $O$ and an $N$-tuple of real numbers $\lambda_1\ge\lambda_2\ge\dots\ge\lambda_N$ as an input and produces a real symmetric matrix as an output.

Note that for a generic $\mathbf X$, $\Phi^{-1}(\mathbf X)$ has $2^N$ elements: $\lambda_1,\dots,\lambda_N$ are eigenvalues of $\mathbf X$ and columns of $O$ are normalized eigenvectors of $\mathbf X$. The latter are defined up to multiplication by $-1$, hence, $2^N$ options. Strictly speaking, this is only true when all $\lambda_i$ are distinct, as otherwise there is more freedom in choosing eigenvectors; however, the $\lambda_i$ coincide on a lower-dimensional set, which has probability $0$ with respect to distributions of interest, and, therefore, is not important for us.

We would like to compute the Jacobian of the transformation $\Phi$. In order to make this precise, note that the set of all $N\times N$ orthogonal matrices is a $N(N-1)/2$ dimensional manifold in $\mathbb R^{N^2}$: it is being cut by $N(N+1)/2$ quadratic equations expressing that the columns of the matrix $O$ are orthonormal. Hence, the space of all orthogonal matrices has a natural measure $\mathrm d O$ inherited from the Lebesgue measure on $\mathbb R^{N^2}$. In parallel, $(\lambda_1,\dots,\lambda_N)$ are equipped with the Lebesgue measure $\mathrm d \lambda_1\cdots \mathrm d \lambda_N$ inherited from $\mathbb R^{N}$ and the symmetric matrices $\mathbf X$ are equipped with the Lebesgue measure $\mathrm d X$ inherited from $\mathbb R^{N(N+1)/2}$, identified with the space of matrix elements on and above the diagonal. The following lemma explains how these measures are related under $\Phi$:

Scheme of the proof of Theorem (ref)

An important part of (ref) is the matrices $\mathbf U \mathbf U^\mathsf T$ and $\mathbf V\mathbf V^\mathsf T$. The first step is to find their distribution.

The matrix $\mathbf Z \mathbf Z^\mathsf T$ is often called the sample covariance matrix in the literature.

Two independent Wishart matrices give rise to an instance of the Jacobi ensemble.

The matrix in (ref) is commonly referred to in the literature as the MANOVA ensemble, with the name derived from Multivariate Analysis of Variance.

Limits of the Jacobi ensemble

In the last chapter we have shown that the distribution of the $K$ squared (sample) canonical correlations between two independent spaces spanned by i.i.d.\ Gaussian matrices of $K\times S$ and $M\times S$ sizes has density proportional to

The goal for this chapter is to analyze the behavior of (ref) in various asymptotic regimes.

Finite $K$ and $M$

We begin by assuming that $K$ and $M$ are fixed and small, while $S$ is large. Based on the interpretations in Section (ref), this means that we have many i.i.d.\ samples of low-dimensional vectors $\mathbf u$ and $\mathbf v$. The following result goes back to hsu1941limiting.

Theorem (ref) has two conceptual consequences. First, it shows that the squared sample canonical correlations converge as $S\to\infty$ to the true squared canonical correlations, which are all zeros in this case, as two data-sets are independent; the speed of convergence is $\frac{1}{S}$. Second, the theorem can be used to construct a statistical test for the hypothesis “two data-sets are independent”:

Theorem (ref) coupled with Theorem (ref) implies that for $\mathbf U$ and $\mathbf V$ with i.i.d.\ Gaussian columns, the independence of $\mathbf U$ and $\mathbf V$ implies that $Sc_1^2>q_\alpha$ with probability close to $1-\alpha$ for large $S$. Hence, for $\alpha$ close to $1$, this is a very unlikely event, which is the reason for the rejection of the hypothesis. In fact, one can relax the Gaussianity assumption in Theorem (ref) (see Problem set 3 in Section (ref)), which implies the asymptotic validity of Procedure (ref) for non-Gaussian $\mathbf U$ and $\mathbf V$. On the other hand, a limitation of the procedure is that the assumption of i.i.d.\ columns can only be slightly relaxed --- uncorrelated columns with some mild technical conditions are permissible --- but cannot be entirely removed.

Large $K$ and $M$

The next step is to investigate the case of large $K$ and large $M$. Figure (ref) shows a histogram of all squared canonical correlations obtained from a single simulation for i.i.d. $\mathcal{N}(0,1)$ matrix elements, $K=100$, $M=150$, $S=500$. The histogram is constructed by splitting the $[0,1]$ interval into bins, counting the squared correlations within each bin, and normalizing so that the area under the histogram equals $1$. While the histogram is random, the simulation reveals a notable feature: Even though the true canonical correlations are all zero, the sample correlations are far from being zero. We will quantify this observation in two theorems below, first by examining all canonical correlations together and then by focusing on the largest ones.

For the first theorem, we represent $K$ squared sample correlations $\hat c_i^2$ by an atomic probability measure, called the empirical measure: $$ \mu_{K,M,S}=\frac{1}{K}\sum_{i=1}^K \delta_{\hat c_i^2}. $$

The asymptotic density of $\mu_{K,M,S}$ is precisely what we see in orange in Figure (ref). It is typical to think of a measure through its pairings with test-functions $f$: $$ \mu_{K,M,S}[f]=\frac{1}{K} \sum_{i=1}^k f(\hat c_i^2). $$ Note that since $\hat c_i^2$ are random, $\mu_{K,M,S}[f]$ is also a random variable. In particular, when $f$ is an indicator function of the interval $(-\infty,t]$, we get the empirical distribution function: $$ \mu_{K,M,S}[ \mathbf 1_{(-\infty,t]}]=\frac{1}{K} \sum_{i=1}^K \mathbf 1_{c_i^2\le t} = \frac{\#\{c_i^2 \text{ which are smaller or equal than } t\}}{K}. $$

Let us introduce the limiting object from Figure (ref), which is called the Wachter distribution, paying tribute to wachter1980limiting.

One can check directly that $0<\lambda_-<\lambda_+<1$ for every $\tau_K>\tau_M>1$ with $\tau_K^{-1}+\tau_M^{-1}<1$ and that (ref) is a probability density, i.e., it integrates to $1$. For any integrable function $f(x)$ we denote $$ \omega_{\tau_K,\tau_M}[f]=\int_0^1 f(x) \omega_{\tau_K,\tau_M}(x) \mathrm d x. $$

There are many approaches to proving Theorem (ref) and analyzing fluctuations (ref), (ref):

We sketch a proof via the last method in the next section. Our approach will crucially rely on the Gaussianity assumption, but other tools allow one to extend Theorem (ref) to non-Gaussian matrix elements of $\mathbf U$ and $\mathbf V$.

Note that while Theorem (ref) describes all canonical correlations together, it does not say anything about their individual limits: for instance, the limits (ref), (ref) will not change if we shift $\hat{c}_1^2$ by an arbitrary number. Hence, we complement Theorem (ref) with the following theorem describing the largest and smallest canonical correlations.\footnote{While the correlations in the middle (or “in the bulk”) are also well-understood, we will not cover them in these lectures, as they are less relevant for the statistical applications we focus on.} We define:

and note that $$ \omega_{\tau_K,\tau_M}(x)\approx \frac{c_\pm}{\pi} \sqrt{|x-\lambda_{\pm}|}, \text{ as } x\to\lambda_{\pm}\quad \text{ inside }\quad [\lambda_-,\lambda_+]. $$

The limiting point process $\{\mathfrak a_i\}_{i=1}^{\infty}$ frequently shows up in random matrix theory: largest (or smallest) eigenvalues in many ensembles of real symmetric random matrices converge to it. The law of the largest coordinate $\mathfrak a_1$ is called the Tracy--Widom distribution $F_1$, and it is even more widespread, appearing also in the limit theorems for interacting particles systems and stochastic PDES, such as the celebrated Kardar-Parisi-Zhang (KPZ) equation. All existing formulas for the finite-dimensional distributions of $\{\mathfrak a_i\}_{i=1}^{\infty}$ are quite complicated. Yet, the distribution can be easily sampled and tabulated, see bornemann2009numerical,vignette_largevars.

We do not provide a proof of Theorem (ref) here. Possible approaches include:

Theorems (ref) and (ref) have two conceptual consequences.

First, they serve as a warning to practitioners applying CCA to real-world data sets: In high-dimensional setting the presence of bounded away from $0$ squared sample canonical correlations does not imply that two data sets are dependent. Care is needed when interpreting CCA results to avoid drawing spurious conclusions.

Second, in parallel with Procedure (ref), Theorems (ref) and (ref) can serve as a foundation for constructing statistical tests for the hypothesis of independence\footnote{For recent work on testing joint independence of more than two datasets see bao2024spectral.} of $\mathbf U$ and $\mathbf V$: one can either use all canonical correlations, relying on (ref), (ref) or use the largest canonical correlations relying on (ref).

Existence of the limit of empirical measure in CCA

We begin by outlining the main ideas behind one particular proof of Theorem (ref). Many technical details in this approach can be reconstructed from guionnet2019asymptotics and references therein.

The first step is to develop the variational principle, which is based on the following rewriting of the Jacobi density $\mathcal J(K,p,q)$:

Identification of the limit

Our next step is to solve the variational problem ${\tilde{\mathfrak E}(\mu)\to \max}$ and show that the maximizer $\tilde \omega_{\tau_K,\tau_M}$ in Lemma (ref) is the Wachter distribution $\omega_{\tau_K,\tau_M}$. We do this by writing down an algebraic equation which uniquely determins $\tilde \omega_{\tau_K,\tau_M}$ and then solve it. It is possible to deduce this equation by directly manipulating the functional $\tilde{\mathfrak E}(\mu)$ (such computations are often done in the logarithmic potential theory, cf.\ saff2013logarithmic). We proceed in a slightly different way and instead derive exact equations for the prelimit object --- Jacobi ensemble $\mathcal J(K,p,q)$. The equation for $\tilde{\mathfrak E}(\mu)$ is then obtained by taking a limit.

Equation (ref) is an instance of the Dyson--Schwinger or loop equation. It can be made a part of a larger hierarchy of similar equations, which uniquely determine the distribution of $\mathcal J(K,p,q)$. Replacing a complicated distribution with a simpler equations for expectations with respect to this distribution is a powerful approach, as we demonstrate on the example of the Gaussian distribution in Problem set 3 in Chapter (ref).

One can find an explicit formula for $\tilde \omega_{\tau_K,\tau_M}(\mathrm d x)$ from (ref) by using the Stieltes tranform of the measure, defined as $$ G(z)=\int_0^1 \frac{1}{z-x} \tilde \omega_{\tau_K,\tau_M}(\mathrm d x),\qquad z\in \mathbb C\setminus \text{support}(\tilde \omega_{\tau_K,\tau_M}). $$

Directly solving equation (ref) for $G(z)$, we get the explicit formula:

We are now ready to finish the proof of Theorem (ref).

The case $\tau_K=\tau_M$ in Theorem (ref) can be treated similarly, but one should be slightly more careful when handling the $x=0$ point: in this case $\lambda_-=0$ and also we might have $p=1$ in Theorem (ref), so that a new boundary term appears in integration by parts.

Signal plus noise

In the previous chapter we studied the asymptotic properties of the sample canonical correlations in the situation when no signal exists: two data sets are independent, so all canonical correlations are zeros. I.e., there was no information to be inferred about the correlation structure. Nevertheless, the asymptotic behavior turned out to be quite rich. In this chapter we make the next step and allow for some signal in the data. That is, we discuss the situation when some of the canonical correlations $c_i$ are not zeros.

Finite $K$ and $M$

As previously, we start by assuming that $K$ and $M$ are fixed and small, while $S$ is large. I.e., we work with matrices $\mathbf U$ and $\mathbf V$ that are composed of many i.i.d.\ samples of low-dimensional vectors $\mathbf u$ and $\mathbf v$. We start by analyzing canonical correlations and then move to canonical variables.

From the statistics perspective, Theorem (ref) tells that $\hat c_i(S)$ are consistent estimates of $c_i$: they approximate the true (population) values of the parameters of the model as $S\to\infty$.

We now turn to canonical variables and seek an appropriate metric to assess how closely the variables derived from the “sample” $\mathbf U$ and $\mathbf V$ approximate those from the “population” $\mathbf u$ and $\mathbf v$. In the finite $K$ and $M$ setting, all proximity metrics are somewhat equivalent; however, selecting the right metric will become crucial in the next subsection, where we allow for large $K$ and $M$.

Let $c_i$ be the $i$--th largest canonical correlation between $\mathbf u$ and $\mathbf v$ and let $\boldsymbol\alpha=(\alpha_1,\alpha_2,\dots,\alpha_K)^\mathsf T$ and $\boldsymbol\beta=(\beta_1,\beta_2,\dots,\beta_M)^\mathsf T$ be the corresponding canonical vectors, i.e., eigenvectors of the matrices $$ (\mathbb E \mathbf u \mathbf u^\mathsf T)^{-1} \mathbb E \mathbf u \mathbf v^\mathsf T (\mathbb E \mathbf v\mathbf v^\mathsf T)^{-1} \mathbb E \mathbf v \mathbf u^{\mathsf T}\quad \text{ and }\quad (\mathbb E \mathbf v \mathbf v^\mathsf T)^{-1} \mathbb E \mathbf v \mathbf u^\mathsf T (\mathbb E \mathbf u\mathbf u^\mathsf T)^{-1} \mathbb E \mathbf u \mathbf v^{\mathsf T}, \qquad \text{respectively.} $$ As shown in Proposition (ref), the random variables $\boldsymbol\alpha^{\mathsf T}\mathbf u$ and $\boldsymbol\beta^\mathsf T\mathbf v$ are canonical variables, corresponding to $c_i$. An advantage of looking at canonical vectors $\boldsymbol\alpha$ and $\boldsymbol\beta$ instead of the canonical variables is that $\boldsymbol\alpha$ and $\boldsymbol\beta$ are non-random and, therefore, can be directly compared with their sample versions $\hat \boldsymbol\alpha (S)$ and $\hat \boldsymbol\beta(S)$, where $\hat \boldsymbol\alpha (S)$ and $\hat \boldsymbol\beta(S)$ are the eigenvectors of $$ (\mathbf U \mathbf U^\mathsf T)^{-1} \mathbf U \mathbf V^\mathsf T (\mathbf V \mathbf V^\mathsf T)^{-1} \mathbf V \mathbf U^\mathsf T \quad \text{ and }\quad (\mathbf V \mathbf V^\mathsf T)^{-1} \mathbf V \mathbf U^\mathsf T (\mathbf U \mathbf U^\mathsf T)^{-1} \mathbf U \mathbf V^\mathsf T, \qquad \text{respectively,} $$ corresponding to the $i$th largest eigenvalue or squared sample canonical correlation $\hat c_i^2(S)$.

Let us emphasize two features of Definition (ref). First, the eigenvectors are only defined up to multiplication by an arbitrary real constant (or up to multiplication by $-1$ if we normalize them), hence, it is tricky to compare them directly; the use of angles avoids this problem.\footnote{The limits of individual coordinates of $\boldsymbol\alpha$ and $\boldsymbol\beta$ are also discussed in the literature, see anderson1999asymptotic and references therein.} Second, instead of measuring the angles between $\widehat \boldsymbol\alpha$ and $\boldsymbol\alpha$ (or $\widehat \boldsymbol\beta$ and $\boldsymbol\beta$), we use the angles between $\mathbf U^\mathsf T \boldsymbol\alpha$ and $\mathbf U^\mathsf T \widehat \boldsymbol\alpha$ (or between $\mathbf V^{\mathsf T}\boldsymbol\beta$ and $\mathbf V^\mathsf T\widehat \boldsymbol\beta$). The reason behind this is that the former angles depend on the choice of units of measurement: the angle changes if we multiply one of the coordinates of $\boldsymbol\alpha$ by a constant (equivalently, divide a component of $\mathbf u$ or a row in $\mathbf U$ by the same constant). In practice, there is generally no default normalization of real data, and this dependence on measurement units implies that, asymptotically, the angles are influenced not only by canonical correlations but also by the covariances of $\mathbf u$ and $\mathbf v$. By focusing on $\theta^\mathbf u_i$ and $ \theta^\mathbf v_i$, we eliminate the issue of normalization, ensuring that the results in the theorems of the next subsection are independent of the covariance matrix of the vector $\mathbf u$ or of the covariance matrix of the vector $\mathbf v$.

Linearly growing $K$ and $M$

For the case of growing $K$ and $M$ we introduce the “signal plus noise” or “spiked random matrix” model, with general terminology going back to johnstone2001distribution. In this model, independent vectors $\mathbf u$ and $\mathbf v$ (or their sample version, arrays $\mathbf U$ and $\mathbf V$) with all canonical correlations $c_i=0$ represent pure noise, meaning no signal is present in the data. Conversely, when some $c_i$ are positive, they and their corresponding canonical variables are treated as signals. Our goal is to separate signal from noise, reconstructing signals from the observed noisy data.

The natural first step toward this goal is to assume that only few $c_i$ are positive, possibly just one. The results from Chapter 3 and, in particular, Figure (ref) indicate that inferring a non-zero correlation from observed data is challenging: the presence of non-zero sample correlations is not sufficient to conclude that even a single $c_i$ is non-zero.

Let us assume that there is only one signal: $c_1=r>0$ and all other $c_i$ equal zero. Thus, $r^2$ equals the single nonzero eigenvalue of the $K\times K$ matrix $(\mathbb E \mathbf u \mathbf u^\mathsf T)^{-1} (\mathbb E\mathbf u \mathbf v^\mathsf T) (\mathbb E\mathbf v \mathbf v^\mathsf T)^{-1} (\mathbb E\mathbf v \mathbf u^\mathsf T)$ and the single nonzero eigenvalue of the ${M\times M}$ matrix $(\mathbb E\mathbf v \mathbf v^\mathsf T)^{-1} (\mathbb E\mathbf v \mathbf u^\mathsf T) (\mathbb E\mathbf u \mathbf u^\mathsf T)^{-1} (\mathbb E\mathbf u \mathbf v^\mathsf T)$. We let $\boldsymbol\alpha$ and $\boldsymbol\beta$ be the corresponding eigenvectors of the former and the latter matrices, respectively. Note that $(\mathbf u^\mathsf T\boldsymbol\alpha, \mathbf v^\mathsf T\boldsymbol\beta, r)$ is the single non-trivial triplet of canonical variables and correlations for $\mathbf u$ and $\mathbf v$.

We further assume that as $S\to\infty$, $r=r(S)$ tends to a constant $\rho>0$. The reason for distinguishing between $r$ and $\rho$ is that the vectors $\mathbf u$ and $\mathbf v$ are now necessarily $S$--dependent, because their dimensions $K$ and $M$ grow with $S$, thus, their correlation structure, captured by $r$, is allowed to change with $S$.

We recall the definition of $\lambda_+$ from (ref) and introduce three additional functions of $\rho$:

The part of Theorem (ref) about the limit of $\hat{c}_i^2$ is proved in bao2019canonical and is further extended to non-Gaussian $\mathbf U$ and $\mathbf V$ in yang2022limiting; another approach to the Gaussian case can be found in hou2023spiked. These articles also show that in case $\mathbf{1}$ the rescaled difference $S^{1/2}(\hat c_1^2-\lambda_+)$ is asymptotically Gaussian with explicit variance. The part about the limit of the angles $\theta^\mathbf u,\,\theta^\mathbf v$ is proved in BG3.

Figure (ref) plots the dependence of $z_\rho$, $\mathfrak s_u$, and $\mathfrak s_\mathbf v$ on $\rho^2$ for a particular choice of $\tau_K$ and $\tau_M$. Analyzing the formulas (ref), (ref), and (ref) one can deduce many features of this dependence. In particular, for $\rho^2<1$, we get $z_\rho>\rho^2$, so that the largest sample canonical correlation always overestimates the true correlation. In addition, $\mathfrak s_\mathbf u>0$ and $\mathfrak s_\mathbf v>0$, which means that the estimates for the canonical variables are not consistent, but rather inclined by certain angles toward the desired true direction, as illustrated in Figure (ref). However, whenever either $\tau_K,\tau_M\to\infty$ or $\rho^2\to 1$, both $z_\rho$ approaches $\rho^2$ and the angles tend to $0$. Therefore, the consistency is restored in such a limit, matching the result of Theorems (ref) and (ref).

When we are working with real data and have no knowledge of the true value of $\rho$, the results of Theorem (ref) should be applied in the following way: if the model matches the data, then most of the squared canonical correlations should belong to the $[\lambda_-,\lambda_+]$ interval. Further, if there is a gap between the largest sample canonical correlation $\lambda_1$ and $\lambda_+$, as in Figure (ref), then using Eq. (ref), we take $\lambda_1$ as an approximation of $z_\rho$. Treating $z_\rho$ as known and approximating $\tau_K$, $\tau_M$ with $S/K$, $S/M$, Eq. (ref) becomes a quadratic equation in $\rho^2$. Solving it and using $ \tfrac{1}{\sqrt{(\tau_M-1)(\tau_K-1)}}\le \rho^2\le 1$ to choose the correct root out of the two, we obtain an estimate for $\rho^2$. Further plugging into (ref), (ref) and using (ref), (ref), we obtain an estimate for the angles between sample and population canonical variables.

Some ingredients of the proof of Theorem (ref)

We briefly outline an approach towards proving Theorem (ref) from BG3. The idea is to develop CCA perturbation theory and then to produce an equation explaining how the canonical correlations and variables change when we add a vector to each of the subspaces.

Suppose that we are given two subspaces, $\widetilde \mathbf U$ and $\widetilde \mathbf V$ in $S$--dimensional space and let $\dim(\widetilde \mathbf V)={M-1}\ge {K-1}=\dim(\widetilde \mathbf U)$. In addition, we are given two vectors $\mathbf u^*$ and $\mathbf v^*$, that we add to spaces $\widetilde \mathbf U$ and $\widetilde \mathbf V$. Define $$ \mathbf U=\mathrm{span}(\mathbf u^*, \widetilde \mathbf U), \qquad \mathbf V=\mathrm{span}(\mathbf v^*, \widetilde \mathbf V). $$ Both spaces $\widetilde \mathbf U$, $\widetilde \mathbf V$ and vectors $\mathbf u^*$, $\mathbf v^*$ are arbitrary, either deterministic or random, and no assumptions are made about them. Our task is to connect the canonical correlations and variables between $\mathbf U$ and $\mathbf V$ to those between $\widetilde \mathbf U$ and $\widetilde \mathbf V$.

Let $\{\widetilde \mathbf u_i\}_{i=1}^{K-1}$, $\{\widetilde \mathbf v_j\}_{j=1}^{M-1}$ and $\{\widetilde c_i\}_{i=1}^{K-1}$ be the canonical correlations and variables between $\widetilde \mathbf U$ and $\widetilde \mathbf V$, so that: $$ \langle \widetilde \mathbf u_i, \widetilde \mathbf u_j\rangle = \delta_{i=j}, \qquad \langle \widetilde \mathbf v_i, \widetilde \mathbf v_j\rangle = \delta_{i=j},\qquad \langle \widetilde \mathbf u_i, \widetilde \mathbf v_j\rangle = \widetilde c_i \delta_{i=j}. $$ We also set $\widetilde c_j=0$ for $K\le j \le M-1$ for convenience of the notations.

Now take two vectors $\widehat \boldsymbol\alpha=(\alpha_0,\alpha_1,\dots,\alpha_{K-1})$ and $\widehat \boldsymbol\beta=(\beta_0,\beta_1,\dots,\beta_{M-1})$, such that $ \alpha_0\mathbf u^*+\sum_{i=1}^{K-1} \alpha_i \widetilde \mathbf u_i$ and $\beta_0 \mathbf v^*+\sum_{j=1}^{M-1} \beta_j \widetilde \mathbf v_j$ is a pair of canonical variables for $\mathbf U$ and $\mathbf V$ corresponding to a squared canonical correlation denoted $z$. We normalize the vectors so that

Then $z=\left\langle \alpha_0 \mathbf u^*+ \sum_{i=1}^{K-1} \alpha_i \widetilde \mathbf u_i,\, \beta_0 \mathbf v^*+ \sum_{j=1}^{M-1} \beta_j \widetilde \mathbf v_j\right\rangle^2$.

Note that previously, in Sections (ref) and (ref), we wrote canonical variables in the bases of rows of the matrices $\mathbf U$ and $\mathbf V$; this is how canonical vectors were defined there. This time, the basis is slightly different: for $\mathbf U$ we are expanding the canonical variables as combinations of $\mathbf u^*$ and $\widetilde \mathbf u_1,\dots,\widetilde \mathbf u_{K-1}$; for $\mathbf V$ we are expanding the canonical variables as combinations of $\mathbf v^*$ and $\widetilde \mathbf v_1,\dots,\widetilde \mathbf v_{M-1}$. Up to this detail, the heuristic meaning of $\widehat \boldsymbol\alpha$ and $\widehat \boldsymbol\beta$ is the same as before: these are still coefficients or coordinates of the expansion of canonical variables.

Multiple signals

Theorem (ref) covers a situation with a single non-zero canonical correlation $c_i$ in the population and one can ask what happens when there are more. The case of finitely many non-zero $c_i$ is well understood, whereas for the case of infinitely many non-zero $c_i$ many questions remain open. In the finite case, when some of the non-zero canonical correlations $c_i$ coincide, the non-uniqueness of canonical variables is again relevant, see, e.g., Remark (ref). In addition, the centered distributional limits of $\hat c_i$ become non-Gaussian, see bao2019canonical. If the canonical correlations are all distinct then essentially the same Theorem (ref) continues to hold for each $c_i$ individually and independently, see bao2019canonical and BG3. The histogram in Figure (ref), referenced in Chapter 1, illustrates multiple signals --- three, to be specific--- using real data. Let us explain that picture in more detail.

We calculate canonical correlations between consumer cyclical and non-cyclical (consumer defensive) stocks. The former are known to follow the state of the economy, while the latter are not related to business cycles and are useful during economic slowdowns. Knowing their correlations can be helpful for portfolio allocation. Hence, we apply CCA and search for maximally correlated combinations of cyclical and non-cyclical stocks. We use weekly returns (which are widely believed to be uncorrelated across time) for the $80$ largest cyclical and $80$ largest non-cyclical stocks over ten years ($01.01.2010-01.01.2020$), which gives us $521$ observations across time; the exact list of stocks is reported in BG3. We can see that the overall shape of the histogram of the correlations matches the Wachter distribution. In view of Theorem (ref) from Chapter 3 and Remark (ref), this indicates that the data matches our modelling assumptions and we can use the formulas of Theorem (ref). At the same time, three largest correlations are clearly separated from the rest. Interestingly, the presence of three correlations is resonant with the celebrated three factors of fama1992cross.

When a linear proportion of the population canonical correlations $c_i$ are nonzero, several new questions arise. In this regime, the Wachter distribution of Theorem (ref) no longer describes the limiting empirical distribution of the squared sample canonical correlations. A more intricate limiting law, which depends on the empirical distribution of the population canonical correlations, has recently been derived in bai2022limiting, hou2023spiked, zhang2025limiting. These works provide explicit formulas for the limiting empirical measure of the sample canonical correlations as a functional of the population counterpart.

From a statistical perspective, the corresponding inverse problem is also of interest: given the observed sample canonical correlations and variables, one aims to recover their population analogues in the setting where many $c_i$ are nonzero. Theorem (ref) suggests that consistent recovery is not possible, yet optimal estimation procedures can still be developed. This question has been extensively studied in the context of Principal Component Analysis (PCA), where such procedures are known as cleaning or shrinkage; see ledoit2022power for a review. Closer to the CCA setting, optimal cleaning methods for cross-covariance matrices were developed in benaych2023optimal. For CCA itself, however, we are not aware of any works establishing analogous optimal shrinkage or cleaning schemes.

Subcritical signals and regularizations

Theorem (ref) requires $\rho^2> \tfrac{1}{\sqrt{(\tau_M-1)(\tau_K-1)}}$, i.e., the correlation of the signal to be large enough, in order to be useful. When $\rho^2$ is smaller, by (ref) the histogram of canonical correlations does not have spikes anymore and the theorem is of no use. The situation $\rho^2<\tfrac{1}{\sqrt{(\tau_M-1)(\tau_K-1)}}$ is called subcritical in the literature, and one can be wondering: is there any way to infer asymptotic information on the signal, i.e., the non-zero canonical correlation and corresponding variables in this case? Is it even possible to detect the presence of a subcrticial signal?

If one sticks to spectral methods, then the most likely answer is “no”. In this direction, assuming the data to be Gaussian, johnstone2020testing computed the asymptotics of the ratio of the probability densities of the canonical correlations under the alternative of a single subcritical signal and under the null of no signal. They show that the ratio neither grows to $\infty$, nor decays to $0$, but rather converges to a non-trivial Gaussian process, parameterized by the value of subcritical $\rho^2$. This is an indication of impossibility of the detection with probability close to $1$ of subcritical $\rho^2$. The transition to subcritical regime at $\rho^2 \approx \tfrac{1}{\sqrt{(\tau_M-1)(\tau_K-1)}}$ for the largest eigenvalues is explored in detail in bykhovskaya2025weak.

For more general class of methods it might be possible to improve the $\tfrac{1}{\sqrt{(\tau_M-1)(\tau_K-1)}}$ threshold and detect some subcritical $\rho^2$ when the noise is non-Gaussian. The details have not yet been figured out in the CCA setting, however, in the related PCA setup improvement can be achieved by applying a certain function (which depends on the distribution on the noise and is the identical function for the Gaussian noise --- leading to no improvements in the latter case) to each matrix element before using spectral methods, see, e.g., perry2018optimality and references therein.

An alternative approach to overcoming the problem of the deterioration of the quality of CCA estimates as the ratio $S/K$ and $S/M$ become smaller is by imposing additional restrictions. If we have extra information about the nature of canonical variables, then we can use it by introducing regularizations into CCA. Some of the popular approaches require the canonical vectors to be sparse (i.e., only a few coordinates of canonical variables are non-zero) or smooth (i.e., the coordinates depend on their labels in a continuous or smooth way). A way to achieve these regularizing requirements is by using the variational formula of Section (ref) for the canonical correlations and variables and then adding extra terms penalizing undesired features to the functions which we maximize. This is a very flexible approach, as various penalties can be used to guarantee various features of regularized canonical variables. However, a drawback is that one needs to rely on the assumptions imposed in constructing the regularization, whose validity might be hard to check in real-world data sets.

There is a large literature on various modifications of CCA and as a possible starting point for explorations of numerous methods, one can consider gao2015minimax,gao2017sparse,uurtio2017tutorial,yang2019survey,tuzhilina2023canonical and references therein.

Cointegration

Beyond its primary objective of identifying correlated components across different data sets, CCA can also be applied to analyze various properties of time series. The idea is to manually construct two rectangular data sets out of a single data set and then implement CCA to infer some features of the data. In particular, in this section we explore how CCA can be used to test for the presence of non-exploding, or cointegrating, linear combinations within an exploding multidimensional time series.

In time series analysis data is typically allowed to be correlated across the $S$ dimension, which in this case represents time and is frequently denoted as $T$. Thus, we now depart from the i.i.d.\ column setting for $\mathbf U$ and $\mathbf V$ used in the previous three chapters. Going forward, the matrices $\mathbf U$ and $\mathbf V$ will exhibit strong dependencies, with the columns of one matrix also being interdependent.

Vector autoregression of order $\mathbf{1}$

The simplest yet nontrivial model for time evolution $X=(X_0,X_1,X_2,\dots,X_T)$ is obtained by considering the vector autoregression model (without an intercept) of order $1$, VAR($1$). The model is based on a sequence of i.i.d. mean zero vector errors $\{\varepsilon_t\}$ with non-degenerate covariance matrix $\Lambda$. The $K$--dimensional process $X_t$ is initialized at a fixed $X_0$ and solves an evolution equation:

where $\Delta X_t:=X_t-X_{t-1}$. Eq. (ref) is called error correction form of a vector autoregression. Eq. (ref) has two (unknown, in principle) deterministic parameters: $K\times K$ positive definite covariance matrix $\Lambda$ and $K\times K$ coefficient matrix $\Pi$. Throughout this chapter we assume that the errors $\varepsilon_t$ are Gaussian, however, many results are either proven or expected to extend to the non-Gaussian setting.

The properties of $X_t$ as a function of $t$ can be very different depending on the matrix $\Pi$. Let us illustrate this in the case $K=1$, so that the evolution equation (ref) becomes

where $\theta$ is a constant corresponding to $I+\Pi$ in (ref). Possible behaviors of (ref) are:

The autoregression model (ref) and its various generalizations are widely adopted as a workhorse for time series modeling, see, e.g., hamilton2020time. The reasons for that are simplicity of the linear model on one hand and diverse possible behaviors on the other hand. The $K$--dimensional extension (ref) has even richer behaviors. In particular, if $\Pi=0$, then $X_t$ is a multidimensional random walk, as in the $\theta=1$ case above. On the other hand, if $\Pi=-I_K$, then $X_t=\varepsilon_t$ is just a sequence of i.i.d.\ random variables, as in the $\theta=0$ case. As we will see below, many intermediate cases are important for applications and have non-trivial theoretical properties.

Cointegration

We start with an informal definition. A $K$--dimensional time-series $X_t$, $t=0,1,2,\dots,$ is said to be cointegrated, if its individual components have statistical properties of a unit root --- like the $\theta=1$ case in (ref) --- but there exists a linear combination of components which is stationary. The number of such linearly independent stationary combinations is called the cointegration rank.

The prefix “co-” in cointegration refers to a sense of “jointly”, as in co-movements. The intuition is that all coordinates follow a long-run equilibrium defined by their stationary linear combination. This is particularly intuitive in the $K=2$ case, where cointegration indicates that the first and second coordinates of $X_t$ tend to move together over time.

The importance of cointegration stems from the celebrated papers granger1981 and engle_granger1987. For example, they argue that monthly rates on $1$-month and $20$-year treasury bonds are cointegrated. A lot of other variables in macroeconomics and finance, such as price level, consumption, output, trade flows, and interest rates are non-stationary, and, thus, are potentially subject to cointegration. Figure (ref) illustrates this point with interest rates on U.S. Treasury bills of three-month and one-year maturities. The time series on the left show wide correlations and can be modeled as unit roots, while the difference of the two series, called spread, on the right exhibits much shorter correlations, consistent with a stationary process.

Going back to the VAR(1) model (ref), we can restate the presence of cointegration as a property of the matrix $\Pi$.

Note that for $r\ge 1$, there is no uniqueness in the choice of the cointegrating relationships, only the linear space spanned by them is uniquely determined: the matrix $\mathbf B$ in Theorem (ref) can be replaced with $Q \mathbf B$ for any non-degenerate $r\times r$ matrix $Q$. Similarly, in the decomposition $\Pi=AB$ we can replace $A\mapsto A Q^{-1}$, $B\mapsto Q B$.

Cointegration testing

For a non-stationary data set, one might ask whether it is cointegrated and, if so, what is the rank of cointegration. That is, following Theorem (ref), we can test if the rank of $\Pi$ is zero (i.e., $\Pi=0$), which would imply the absence of cointegration.

Similarly to Task (ref), we can approach this new question via Gaussian Maximum Likelihood. Namely, assuming that the $\varepsilon_t$ in (ref) are Gaussian and given observations $X=(X_0,X_1,\dots,X_T)$, we evaluate the likelihood function $L(X; \Lambda,\Pi)$, which computes the probability density at $X$. Next, we maximize this function over $\Lambda$ and $\Pi$ under the assumption $\mathrm{rank}(\Pi)\le r$, thus getting

The definition readily implies the monotonicity: $$ L_0(X)\le L_1(X)\le \dots\le L_K(X). $$ We expect that if $\Pi=0$ already fits the data well, maximizing over additional choices of $\Pi$ should not result in a substantial increase in $L_r(X)$ as a function of $r$. On the other hand, if the true rank of $\Pi$ is non-zero, a significant increase in $L_r(X)$ is anticipated. This leads us to the following statistical procedure, known as the likelihood ratio test.

The implementation of Procedure (ref) involves two steps. The first step is to compute (ref) as a function of $X$. This explicit computation goes back to Anderson and johansen1988, and we derive it in Theorem (ref) below. The second step requires understanding the typical magnitude of (ref) to be able to interpret when it becomes significantly large.

To understand why CCA might be relevant to identifying the rank of $\Pi$, let’s revisit (ref). Suppose that $\varepsilon_t$ is very small, while $\Pi$ has rank $1$. Then $\Delta X_t=\Pi X_{t-1}+\varepsilon_t$ implies that there is a linear combination of rows of $X_{t-1}$, which is nearly proportional to each row of $\Delta X_t$, resulting in a canonical correlation close to $1$. Similarly, if $\Pi$ has rank $r$ (still with $\varepsilon_t$ small), there would be $r$ canonical correlations close to $1$. This suggests that the number of large canonical correlations corresponds to the rank of $\Pi$. Theorem (ref) formalizes this observation and extends it to cases where $\varepsilon_t$ is not negligible.

For the second step in implementing Procedure (ref), we draw on the ideas that partly parallel the developments in Chapter 3. Specifically, we analyze the asymptotic behavior of the sum $\sum_{i=1}^r \ln(1-\lambda_i)$ in two regimes: for fixed $K$ and large $T$, and for $K$ and $T$ being simultaneously large.

Finite $K$, large $T$ asymptotics

We start from the classical asymptotic regime where the time series dimension $K$ is small and fixed, while the number of observations $T$ is large and is assumed to tend to infinity. Here, we aim to analyze $\sum_{i=1}^r \ln(1 - \lambda_i)$ under the null hypothesis $\Pi = 0$, which helps define “large” in Procedure (ref), and under alternatives with non-trivial $\Pi$, to assess the asymptotic power of the test. Let us start with the null hypothesis case.

Corollary (ref) clarifies the term “large” in Procedure (ref): when $K$ is small and $T$ is large, the likelihood ratio (ref) should be compared with quantiles of the random variable $-\frac{1}{2} \sum\limits_{i=1}^r \nu_i$.

The procedure’s effectiveness --- rejecting the null of no cointegration when $\Pi \neq 0$ --- relies on the fact that asymptotically $I(0)$ (stationary) and $I(1)$ (unit root) variables are uncorrelated: The covariance between them remains bounded, $I(0)$ variables have bounded variance, while $I(1)$ variables' variance grows to infinity over time. When $\Pi=0$, no stationary linear combination $\beta^{\mathsf T}X_{t-1}$ exists, and we are left with correlations between various $I(0)$ and $I(1)$ variables, which causes all sample canonical correlations to vanish asymptotically. In contrast, when $\Pi\neq0$ a stationary linear combination exists and the corresponding sample canonical correlation does not vanish.

Analyzing the procedure under $\Pi\neq0$ corresponds to assessing its power. Theorem (ref) and Corollary (ref) can be extended to allow for this scenario. Additionally, extensions accommodate higher-order autoregressions (i.e., including $X_{t-1},X_{t-2},\dots,X_{t-k}$ rather than only $X_{t-1}$ in (ref)) and inclusion of additional deterministic terms, see johansen1988,johansen1991.

Large $K$ regime

Soon after Theorem (ref) and its extensions were discovered, practitioners noticed that $K$ needs to be much smaller than $T$, in order for the result to be useful. The distance between finite sample and asymptotic limit is too large for intermediate values of $K$, e.g., for $K=10$, $T=100$, see, for example, ho_sorensen1996 and gonzalo_pitarakis. One possible way to address this distortion is by considering a different asymptotic regime. The development of alternative asymptotics based on $K$ and $T$ growing to infinity proportionally to each other is recent, see onatski_ecta, BG1, BG2. It turns out that the asymptotic behavior is governed by the same objects as in Theorems (ref) and (ref) of Chapter 3: Wachter distribution and Airy$_1$ random point process. However, the parameters change.

We illustrate Theorem (ref) using a real data set. Figure (ref) in Chapter 1 shows a histogram of the squared sample canonical correlations $\lambda_i$, calculated for weekly stock returns in S&P$100$ over ten years, so that $K=92$, $T=521$, see BG1 for the data description. The histogram shows a striking alignment with the density of Wachter distribution $\omega_{1+\tau,\,(1+\tau)/2}$ (orange curve), where $\tau=521/92\approx 5.66$.

Note that one should be careful with applying Theorem (ref) to $\sum_{i=1}^r \ln(1-\lambda_i)$ needed for Procedure (ref). The reason is that $\ln(1-x)$ explodes at $x=1$ and a priory, Theorem (ref) only implies a one-sided bound on $\lim_{K,T\to\infty} \frac{1}{K}\sum_{i=1}^{r}\ln(1-\lambda_i)$ in the regime of $r$ growing linearly with $K,T$, cf.\ onatski_ecta and BG2. In order to obtain more precise asymptotics, one needs better control on $\lambda_1$; this control was achieved in BG1 (see also onatski2019extreme) at the expense of slightly modifying the procedure for producing $\lambda_i$.

This modified procedure is discussed in detail in BG1. Conceptually, new detrending and demeaning steps remove potential non-zero means in the noise terms $\varepsilon_t$ making the approach quite natural (in practice there is no guarantee that the true mean is zero). Technically, these steps enable us to establish an asymptotic theorem:

Again, (ref) can be immediately applied to clarify the term “large” in Procedure (ref) (modified version of Procedure (ref) based on $\tilde\lambda_i$). One should use the quantiles of the random variables $\sum_{i=1}^r \mathfrak a_i$ to deduce what values are atypically large. The densities of these random variables for $r\le 3$ are plotted in Figure (ref) and tables of distribution functions for $r\le 10$ are contained in vignette_largevars.

Theorem (ref) was first proved in onatski_ecta by using delicate manipulations with resolvents to obtain a system of four equations for the Stieltjes transform of $\lim_{K,T\to\infty} \mu_{K,T}$ and then identifying the solution as the Wachter distribution. In this approach the Gaussianity of $\varepsilon_t$ is not important, and it allows slight modifications of the procedure (e.g., works both for $\lambda_i$ and $\tilde \lambda_i$). An alternative proof of Theorem (ref), as well as the proof of Theorem (ref) was obtained in BG1 exploring a different set of ideas, for which both Gaussianity of $\varepsilon_t$ and the precise details of Procedure (ref) are important\footnote{We expect Theorem (ref) to extend to much wider generality, but this was not proven as of 2024.}. The latter approach builds on a connection to Jacobi ensemble:

Let us emphasize the differences between Theorem (ref) and the previous appearance of the Jacobi ensemble in CCA in Theorem (ref). In the latter theorem, the rectangular matrices $\mathbf U$ and $\mathbf V$ are independent; in contrast, in the former theorem one matrix is a function of another: up to minor corrections, $\tilde \mathbf U$ of Procedure (ref) is the time-increment of $\tilde \mathbf V$. Moreover, Theorem (ref) is an exact match with the Jacobi ensemble, while Theorem (ref) is only an approximation (and no exact match is known or expected). The idea of the proof of Theorem (ref) is to replace the time-shift operator implicitly entering Procedure (ref) by a uniformly random orthogonal operator and check that this replacement does not change the canonical correlations much, see BG1 for details.

Once Theorem (ref) is established, Theorems (ref) and (ref) become applications of the asymptotic results for the Jacobi ensemble which we discussed in Chapter 3.

Outlook: signal plus noise for cointegration

So far we only explored the large $K$, $T$ asymptotics of the objects related to cointegration tests under the null hypothesis of no cointegration, $\Pi=0$. Drawing from the ideas in Chapter 4, one might hope to develop a signal-plus-noise framework for cases where $\Pi\ne 0$. As of 2024, such a theory has yet to be fully established. Nevertheless, Monte Carlo simulations (see Figure (ref)) are encouraging and suggest that this theory may soon emerge. Looking at the real data, the presence of three correlations extending to the right of the support of the Wachter distribution for cryptocurrencies in Figure (ref) of Chapter 1 likely signals the presence of three cointegrating relationships.

Exercises

Problem set 0 --- before Chapter 1

A refresher on the linear algebra

{\bf Problem 1.} Consider two pairs of lines in $\mathbb R^n$: $(\mathbf u,\mathbf v)$ and $(\mathbf u',\mathbf v')$. Suppose that the angle between $\mathbf u$ and $\mathbf v$ is $0\le \phi\le \frac{\pi}{2}$ and the angle between $\mathbf u'$ and $\mathbf v'$ is $0\le \phi'\le \frac{\pi}{2}$. Prove that $\phi=\phi'$ if and only if there exists an orthogonal transformation of $\mathbb R^n$ (i.e., $n\times n$ orthogonal matrix $O$), which maps $\mathbf u$ to $\mathbf u'$ and $\mathbf v$ to $\mathbf v'$.

{\bf Problem 2.} Let $A$ be a $k\times n$ real matrix. Treat the rows of this matrix as vectors in $\mathbb R^{n}$, let $\mathbf U\subset R^n$ be the linear space spanned by the rows of $A$, and let $P_{\mathbf U}$ be the orthogonal projection on $\mathbf U$. Viewing $P_\mathbf U$ as a $n\times n$ matrix, find an expression for it in terms of $A$.

{\bf Remark.} For simplicity, feel free to assume $k\le n$ and that $A$ has rank $k$. If you want to challenge yourself, you can also try to develop a formula for the case when rank of $A$ is smaller than $k$.

{\bf Problem 3.} Consider two lines in $\mathbb R^n$, $\mathbf u$ and $\mathbf v$ with angle $\phi$ between them. Let $P_\mathbf u$ be orthogonal projection on $\mathbf u$ and let $P_\mathbf v$ be orthogonal projection on $\mathbf v$. Prove that $P_\mathbf v P_\mathbf u$ has a unique non-zero eigenvalue and this eigenvalue equals $\cos^2\phi$.

{\bf Problem 4.} Let $\xi_1,\xi_2,\dots,\xi_n$ be a random vector. Assume that $\mathbb E \xi_i=0$ and $\mathbb E \xi_i^2<\infty$ for all $1\le i \le n$. Prove that there exists a deterministic $n\times n$ matrix $A$ and a random vector $\eta_1,\eta_2,\dots,\eta_n$, such that:

Problem set 1 --- between Chapters 1 and 2

Transformations of Gaussian vectors

{\bf Problem 1.} Compute $\displaystyle \int_{-\infty}^{\infty} e^{-x^2}\, \mathrm d x$.

{\bf Hint.} One possible approach is to square the expression and treat it as a 2-dim integral.

{\bf Problem 2.} Given a $n\times n$ positive-definite matrix $\Phi$, we define an $n$--dimensional mean $0$ Gaussian vector $(\xi_1,\xi_2,\dots,\xi_n)$ through its probability density with respect to the Lebesgue measure: $$ \rho(x_1,\dots,x_n)= \frac{1}{Z(\Phi)} \exp\left( - \frac{1}{2}

\Phi

\right) $$

{\bf Notation.} Such a Gaussian vector is denoted $\mathcal N(0,\Lambda)$, where $0$ indicates that the mean is $0$ (and alternatively, if we add a deterministic vector $\mu$, then the mean would be $\mu$), and $\Lambda$ is the covariance matrix. In particular, in the one-dimensional case $\mathcal N(0,1)$ is the standard Gaussian random variable of mean $0$ and variance $1$.

{\bf Problem 3.} Let $(\xi_1,\dots,\xi_n)$ be a random vector with i.i.d.\ $\mathcal N(0,1)$ components. Find the probability density function of $\xi_1^2+\xi_2^2+\dots+\xi_n^2$.

{\bf Hint.} The answer for $n=2$ is particularly simple.

{\bf Problem 4.} Let $(\xi_1,\dots,\xi_n)$ be a vector with i.i.d.\ $\mathcal N(0,\sigma^2_{\xi})$ components and let $(\eta_1,\dots,\eta_n)$ be another independent vector with i.i.d.\ $\mathcal N(0,\sigma^2_{\eta})$ components. Find the probability distribution (its density function) of the squared correlation coefficient between $\xi$ and $\eta$, which is $$ \frac{(\xi_1 \eta_1 + \xi_2 \eta_2+\dots +\xi_n \eta_n)^2}{(\xi_1^2+\xi_2^2+\dots+\xi_n^2) \cdot (\eta_1^2+\eta_2^2+\dots+\eta_n^2)}. $$

{\bf Hint.} You might try to rotate $\mathbb R^{n}$ and show first that the distribution of this squared correlation coefficient is the same as the distribution of the squared correlation coefficient between vector $\xi$ and deterministic vector $(1,0,0,\dots,0)$.

Problem set 2 --- between Chapters 2 and 3

MLE and matrix algebra

{\bf Problem 1.} Let $\xi_1,\xi_2,\dots,\xi_S$ be i.i.d.\ $\mathcal N(0,\sigma^2)$, where $\sigma^2$ is unknown parameter, which we would like to estimate. Write down the joint density of $\xi_1,\dots,\xi_S$ as a function $L(x_1,\dots,x_S;\sigma^2)$. Find the Gaussian MLE $\hat \sigma^2$ by maximizing $$ \hat \sigma^2=\mathrm{argmax}_{\sigma^2}\left[ L(\xi_1,\dots,\xi_S;\sigma^2)\right]. $$ Use the Law of Large Numbers to prove that $\lim_{S\to\infty} \hat \sigma^2 = \sigma^2$.

{\bf Problem 2.} Suppose that $D$ is a positive-definite symmetric $N\times N$ matrix and let $\gamma>0$. Consider a function $$ f(G)=-\gamma\ln\det(G)-\mathrm{Trace}(G^{-1} D). $$ Prove that the maximum of $f(G)$ over all $N\times N$ positive-definite symmetric $N\times N$ matrices occurs at $G=\frac{1}{\gamma}D$.

{\bf Hint.} First, consider $D=I_N$ and represent trace and determinant through eigenvalues. Then reduce general $D$ to $D=I_N$ by a change of variables.

{\bf Problem 3.} Let $\xi_1,\xi_2,\dots,\xi_S$ be i.i.d.\ $N$-dimensional vectors $\mathcal N(0,\Lambda)$, where $\Lambda$ is an unknown covariance matrix, which we would like to estimate. Write down the joint density of $\xi_1,\dots,\xi_S$ as a function $L(x_1,\dots,x_S; \Lambda)$. Find the Gaussian MLE $\hat \Lambda$ by maximizing $$ \hat \Lambda=\mathrm{argmax}_{\Lambda}\left[ L(\xi_1,\dots,\xi_S; \Lambda)\right]. $$ Use the Law of Large Numbers to prove that $\lim_{S\to\infty} \hat \Lambda = \Lambda$.

{\bf Hint.} Use Problem 2 to find $\mathrm{argmax}$.

{\bf Problem 4.} Let $R$ be a $N\times N$ matrix. Define $$ O=\exp(R)=I_N + R + \frac{R^2}{2!} + \frac{R^3}{3!} + \dots. $$

Problem set 3 --- between Chapters 3 and 4

Limits under independence

{\bf Problem 1.} Let $\xi$ and $\eta$ be two independent mean $0$ random variables with finite fourth moments. Let $\xi_1,\xi_2,\dots$ and $\eta_1,\eta_2,\dots$ be i.i.d.\ samples of $\xi$ and $\eta$. Let $\rho^2(S)$ denote the squared sample correlation coefficient between $(\xi_1,\dots,\xi_S)$ and $(\eta_1,\dots,\eta_S)$. Using LLN and CLT investigate the behavior of $\rho^2(S)$ as $S\to\infty$.

{\bf Hint.} You should get the same answer as in the Gaussian case discussed in Chapter 3.

{\bf Problem 2.} Recall that we showed in Chapter 2 that in the Gaussian case the single squared sample canonical correlation $\hat{\mathfrak c}^2$ between $1\times S$ matrix $\mathbf U$ and $M\times S$ matrix $\mathbf V$ is distributed as $\mathcal J(1; \frac{M}{2},\frac{S-M}{2})$, i.e., with Beta-density proportional to

Suppose that $M,S\to\infty$ in such a way that $S/M\to \tau_M>1$. Prove that there exist constants $a$ and $b$ such that $$ \lim_{M,S\to\infty} \left[\frac{\hat{\mathfrak c}^2 - a}{b}\sqrt{S}\right]\stackrel{d}{=} \mathcal N(0,1). $$ Find formulas for $a$ and $b$ in terms of $\tau_M$.

{\bf Hint.} Taylor-expand the density near its maximum.

{\bf Problem 3.} Extend the asymptotics of Problem 2 to the case when $\mathbf U$ is $K\times S$ and $K$ stays fixed as $S\to\infty$.

{\bf Problem 4.} Let $\xi$ be a random variable $\mathcal N(m,\sigma^2)$. Let $\mathbf {Poly}$ denote the class of all continuously differentiable functions $f(x)$ on $\mathbb R$, which grow at most polynomially as $x\to\infty$.

Problem set 4 --- between Chapters 4 and 5

Spiked GOE (simpler version of Theorem (ref))

{\bf Problem 1.} Take $N\times N$ real symmetric matrix $\mathbf B$ with distinct eigenvalues $\lambda_1>\lambda_2>\dots>\lambda_N$ and normalized eigenvectors $\mathbf u_i$, $1\le i \le N$, satisfying $\langle \mathbf u_i, \mathbf u_j \rangle =\delta_{i=j}$. In addition, take a column-vector $\mathbf u^*$. For a real constant $\theta$ define $$ \mathbf A=\theta\cdot \mathbf u^* (\mathbf u^*)^\mathsf T + \mathbf B. $$ Suppose that $\langle \mathbf u_i, \mathbf u^*\rangle\ne 0$ for all $i$. Prove that each eigenvalue $a$ of $\mathbf A$ satisfies an equation $$ \frac{1}{\theta }=\sum_{i=1}^N \frac{ \langle \mathbf u_i, \mathbf u^*\rangle^2}{a-\lambda_i}. $$

{\bf Hint.} Search for an eigenvector of the form $\sum_{i=1}^N \alpha_i \mathbf u_i$.

{\bf Problem 2.} As continuation of Problem 1, suppose that $\mathbf u^*$ has i.i.d.\ $\mathcal N(0,\frac{1}{N})$ components. Prove that $ \langle \mathbf u_i, \mathbf u^*\rangle$, $i=1,2,\dots,N$ are also i.i.d.\ $\mathcal N(0,\frac{1}{N})$.

{\bf Problem 3.} Suppose that $\mathbf B$ is the rescaled GOE, which means that $\mathbf B=\frac{1}{2\sqrt{N}}(X+X^*)$, where $X$ is a matrix of i.i.d.\ $\mathcal N(0,1)$. You can use without proof the semicircle law in the form: if $\lambda_1\ge \dots\ge \lambda_N$ are eigenvalues of $\mathbf B$ and $G_N(z)=\frac{1}{N}\sum_{i=1}^N \frac{1}{z-\lambda_i}$, then: $$\lim_{N\to\infty} \lambda_1=2, \qquad \lim_{N\to\infty}\lambda_N = -2, \qquad \lim_{N\to\infty} G_N(z)=\frac{1}{2}\bigl(z-\sqrt{z^2-4}\bigr).$$ Suppose $\theta>1$. Prove that (for $\mathbf u^*$ as in Problem 2) as $N\to\infty$, the relation between $\theta$ and $a>2$ of Problem $1$ becomes $a=\theta+\frac{1}{\theta}$.

{\bf Problem 4.} In the setting of Problem 3 compute the limit as $N\to\infty$ of the angle between $\mathbf u^*$ and the eigenvector of $\mathbf A$ with the largest eigenvalue $a\approx \theta + \frac{1}{\theta}$. How does it behave as $\theta\to 1$? As $\theta\to+\infty$?

Problem set 5 --- after Chapter 5

Examples and properties of autoregressions

{\bf Problem 1.} Let $\varepsilon_t$ be i.i.d.\ $\mathcal N(0,1)$. Fix a deterministic parameter $\theta$ and let $x_t$, $t=0,1,2,\dots$ be a solution to $$ x_t=\theta x_{t-1}+\varepsilon_t,\qquad t=1,2,\dots,\qquad x_0=0. $$ For all values of $\theta$, compute explicitly the correlation coefficient between $x_t$ and $x_{t+s}$ as a function of $t$ and $s$. Prove that the coefficient:

{\bf Problem 2.} In the setting of previous problem with $|\theta|<1$, find out what should be the distribution of (random) $x_0=\xi$, so that $x_t$ would become a stationary process in $t$. Prove that for any other initial condition $x_0$, one has $$ \lim_{t\to\infty} x_t\stackrel{d}{=} \xi. $$

{\bf Problem 3.} Suppose that a two-dimensional vector $X_t$ solves an evolution equation: $$ \Delta X_t=

X_{t-1} + \varepsilon_t, \qquad \Delta X_t = X_t- X_{t-1}, \qquad t=1,2,\dots, $$ where $\varepsilon_t$ are i.i.d.\ with independent $\mathcal N(0,1)$ components. Show that $X_t$ is cointegrated, i.e.,

{\bf Problem 4.} Suppose that a two-dimensional vector $X_t$ solves an evolution equation: $$ \Delta X_t=

X_{t-1} + \varepsilon_t, \qquad \Delta X_t = X_t- X_{t-1}, \qquad t=1,2,\dots, $$ where $\varepsilon_t$ are i.i.d.\ with independent $\mathcal N(0,1)$ components. Compute the speed of growth of the variances of each of the two components of $X_t$ as $t\to\infty$.

{\bf Remark.} We say that one of the components of $X_t$ is $I(1)$ process, while another one is $I(2)$ process; the same terminology refers to the stationary processes as $I(0)$. The value $d$ in $I(d)$ refers to how many times one needs to apply $\Delta$ to $X_t$ to transform it to a stationary process.