Computing Canonical Roots From Sample Covariances

From Paired Observations to Canonical Roots

Canonical correlation analysis starts with paired observations on two variable sets. Computing its sample roots requires more than entering a formula into an eigensolver: the covariance blocks must be formed using a consistent centering and divisor, the eigenvalues must be ordered and interpreted correctly, and the resulting directions must satisfy the intended normalization. These checks are especially important when computations use rounded data or software output.

Let $x_i\in\mathbb{R}^p$ and $y_i\in\mathbb{R}^q$ denote the two sets observed on the same $n$ units. Write $\bar{x}$ and $\bar{y}$ for their sample means. Use the usual sample covariance divisor $n-1$, and assume the within-set sample covariance matrices are positive definite when solving the generalized eigenproblem below.

Definition: The sample within-set covariance matrices and cross-set covariance matrix are $$ S_{XX}=\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})(x_i-\bar{x})^\top,\qquad S_{YY}=\frac{1}{n-1}\sum_{i=1}^{n}(y_i-\bar{y})(y_i-\bar{y})^\top, $$ $$ S_{XY}=\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y})^\top,\qquad S_{YX}=S_{XY}^\top. $$ The matrix $S_{XX}$ is $p\times p$, $S_{YY}$ is $q\times q$, and $S_{XY}$ is $p\times q$. The joint sample covariance matrix has blocks $S_{XX},S_{XY},S_{YX},S_{YY}$.

Each row must pair the $X$- and $Y$-measurements from the same unit. The same $n$ observations and the same divisor must be used for all four blocks. Once the blocks are assembled, the squared canonical roots are the generalized eigenvalues $\lambda$ in $S_{XY}S_{YY}^{-1}S_{YX}a=\lambda S_{XX}a$. The canonical correlations are the nonnegative square roots $\rho=\sqrt{\lambda}$.

Why the Generalized Eigenproblem Gives the Roots

The two canonical directions $a$ and $b$ are scaled so that their sample scores have unit variance: $a^\top S_{XX}a=1$ and $b^\top S_{YY}b=1$. At a canonical pair with correlation $\rho$, their sample covariance is $a^\top S_{XY}b=\rho$. The corresponding stationarity equations, established in the earlier material on canonical correlation, are

S_{XY}b=\rho S_{XX}a,\qquad S_{YX}a=\rho S_{YY}b.

The first equation gives $a$ for a proposed $b$, and the second gives $b$ for a proposed $a$. Eliminating $b$ yields the eigenproblem used in computation.

Theorem: Suppose $S_{XX}$ and $S_{YY}$ are positive definite and arise as the diagonal blocks of a positive semidefinite joint sample covariance matrix. The squared canonical roots are the generalized eigenvalues of $$ S_{XY}S_{YY}^{-1}S_{YX}a=\lambda S_{XX}a. $$ They are real and nonnegative. They can be ordered as $\lambda_1\geq\lambda_2\geq\cdots\geq 0$; the canonical correlations are $\rho_j=\sqrt{\lambda_j}$.

Proof. The second stationarity equation gives $b=\rho^{-1}S_{YY}^{-1}S_{YX}a$ for a positive root. Substituting into the first gives $S_{XY}S_{YY}^{-1}S_{YX}a=\rho^2S_{XX}a$. Thus $\lambda=\rho^2$. To verify reality and nonnegativity, define $C=S_{XX}^{-1/2}S_{XY}S_{YY}^{-1/2}$, where the inverse square roots exist by positive definiteness. The matrix in the ordinary eigenproblem corresponding to the generalized one is $S_{XX}^{-1/2}S_{XY}S_{YY}^{-1}S_{YX}S_{XX}^{-1/2}=CC^\top$. This matrix is symmetric positive semidefinite because, for every vector $v$, $v^\top CC^\top v=\|C^\top v\|^2\geq 0$. Its eigenvalues are therefore real and nonnegative, and these are precisely the generalized eigenvalues. Sorting them in descending order sorts the squared canonical correlations in descending order as well. $\square$

For a computed generalized eigenvector $a$ with $\lambda>0$, one can recover the paired direction from $b=\lambda^{-1/2}S_{YY}^{-1}S_{YX}a$. Rescale the directions so that both score variances equal one, and choose their signs together so the score covariance is nonnegative. A numerical solver may return eigenvectors in arbitrary scale and sign, so neither normalization nor ordering should be assumed from the raw output.

Worked Computations

Worked Example: Building Covariance Blocks From Paired Data

Consider five paired observations described by four centered columns $U,V,W,Z$. Their entries are specified by $U=\sqrt{2}(1,-1,0,0,0)^\top$, $V=\frac{2}{\sqrt{6}}(1,1,-2,0,0)^\top$, $W=\frac{1}{\sqrt{3}}(1,1,1,-3,0)^\top$, and $Z=\frac{1}{\sqrt{5}}(1,1,1,1,-4)^\top$. Take $X_1=U,\ X_2=V,\ Y_1=0.8U+0.6W,\ Y_2=0.3V+\sqrt{0.91}Z$. These formulas specify the paired values for each of the five units.

Each column sums to zero. Their pairwise inner products are zero, and each squared length is four: for example, $U^\top U=2(1^2+(-1)^2)=4$, while $V^\top V=(4/6)(1+1+4)=4$. The same calculation gives $W^\top W=(1/3)(1+1+1+9)=4$ and $Z^\top Z=(1/5)(1+1+1+1+16)=4$. Since $n-1=4$, each column has sample variance one, and distinct columns have sample covariance zero.

Consequently $S_{XX}=I_2$. The $Y$-columns also have sample variance one: $0.8^2+0.6^2=1$ and $0.3^2+0.91=1$. They are uncorrelated because they use disjoint orthogonal columns. The cross-covariances give

S_{YY}=I_2,\qquad S_{XY}= \begin{pmatrix} 0.8&0\\ 0&0.3 \end{pmatrix}.

Thus $S_{XY}S_{YY}^{-1}S_{YX}=\operatorname{diag}(0.64,0.09)$, and $S_{XX}=I_2$. The generalized eigenvalues are $0.64$ and $0.09$, so the ordered canonical correlations are $0.8$ and $0.3$. The corresponding directions are the two coordinate vectors in each set, and their score variances are one. This calculation also checks that the larger root is listed first.

Worked Example: Solving a Two-Variable-Set Generalized Eigenproblem

Suppose the sample covariance blocks are $S_{XX}=\begin{pmatrix}4&1\\1&2\end{pmatrix}$, $S_{YY}=(9)$, and $S_{XY}=\begin{pmatrix}2\\3\end{pmatrix}$. Both within-set covariance matrices are positive definite: $S_{XX}$ has positive leading entry and determinant $8-1=7>0$, and $S_{YY}=9>0$.

The inverse of $S_{XX}$ is $\frac{1}{7}\begin{pmatrix}2&-1\\-1&4\end{pmatrix}$. The one-dimensional $Y$-set makes the root calculation especially direct:

S_{XY}S_{YY}^{-1}S_{YX} =\frac{1}{9}\begin{pmatrix}2\\3\end{pmatrix}(2\ \ 3) =\frac{1}{9}\begin{pmatrix}4&6\\6&9\end{pmatrix}.

The nonzero generalized eigenvalue is $\lambda=\frac{1}{9}S_{YX}S_{XX}^{-1}S_{XY}$. In detail, $S_{XX}^{-1}S_{XY}=\frac{1}{7}(1,10)^\top$, so $S_{YX}S_{XX}^{-1}S_{XY}=(2,3)\frac{1}{7}(1,10)^\top=32/7$? Computing the dot product carefully gives $(2+30)/7=32/7$. Hence $\lambda=32/63$, and $\rho=\sqrt{32/63}\approx0.7127$.

An eigenvector direction is proportional to $S_{XX}^{-1}S_{XY}=(1,10)^\top$. Its squared score scale is $(1,10)S_{XX}(1,10)^\top=4+20+200=224$, so set $a=(1,10)^\top/\sqrt{224}$. The paired $Y$-direction can be chosen as $b=1/3$, which has variance $b^2S_{YY}=1$. The score covariance is $a^\top S_{XY}b=((2+30)/\sqrt{224})(1/3)=32/(3\sqrt{224})\approx0.7127$, agreeing with the root. The joint covariance is admissible: $S_{YY}-S_{YX}S_{XX}^{-1}S_{XY}=9-32/7=31/7>0$.

Worked Example: Checking Score Normalization and a Cross-Covariance

Let $S_{XX}=S_{YY}=I_2$ and $S_{XY}=\begin{pmatrix}0.5&0\\0&0.2\end{pmatrix}$. The generalized eigenproblem has roots $0.25$ and $0.04$, so the leading canonical correlation is $0.5$. Consider the directions $a=(1/\sqrt{2},0)^\top$ and $b=(1,0)^\top$. Their score variances are $a^\top S_{XX}a=1/2$ and $b^\top S_{YY}b=1$, so $a$ is not yet normalized.

The sample covariance of these unnormalized scores is $a^\top S_{XY}b=(1/\sqrt{2})(1/2)=1/(2\sqrt{2})\approx0.3536$, not $0.7071$. To normalize the first score, multiply $a$ by $\sqrt{2}$, giving $a_*=(1,0)^\top$. Then $a_*^\top S_{XX}a_*=1$, and $a_*^\top S_{XY}b=0.5$, the leading canonical correlation. This distinguishes a raw score covariance from the correlation obtained after both score variances are set to one.

For an exact residual check, the leading direction satisfies $S_{XY}S_{YY}^{-1}S_{YX}a_*=0.25a_*$, because the left side is $\operatorname{diag}(0.25,0.04)(1,0)^\top=(0.25,0)^\top$. Its generalized-eigenproblem residual is therefore the zero vector. The second root is smaller, as required by descending root order.

Numerical Checks and Common Pitfalls

In an actual calculation, form the centered data matrices and compute all covariance blocks from them, or verify supplied blocks against that construction. If $X_c$ and $Y_c$ have centered observations in rows, then $S_{XX}=X_c^\top X_c/(n-1)$, $S_{YY}=Y_c^\top Y_c/(n-1)$, and $S_{XY}=X_c^\top Y_c/(n-1)$. Check matrix dimensions before multiplication: $S_{XY}S_{YY}^{-1}S_{YX}$ must be $p\times p$.

For a candidate eigenpair $(\lambda,a)$, compute the generalized residual $r=S_{XY}S_{YY}^{-1}S_{YX}a-\lambda S_{XX}a$. A small residual relative to the sizes of the two terms supports the numerical solution; for example, one can report $\|r\|/(\|S_{XY}S_{YY}^{-1}S_{YX}a\|+|\lambda|\|S_{XX}a\|)$ when the denominator is nonzero. Then verify $a^\top S_{XX}a=1$, recover $b$, and verify $b^\top S_{YY}b=1$ and $a^\top S_{XY}b=\sqrt{\lambda}$, allowing for numerical tolerance. There is no universal residual cutoff independent of scale and software precision.

A frequent error is to report eigenvalues as canonical correlations. The eigenvalues are squared correlations, so a root of $0.25$ means a canonical correlation of $0.5$, not $0.25$. Another is to sort by the magnitude of an unsquared quantity or leave the solver's output order unchecked. Report the nonnegative roots in descending order, retaining the associated directions with each root.

Key takeaway: Build every covariance block from the same centered paired observations and divisor; solve for generalized eigenvalues, take their nonnegative square roots, sort them in descending order, and verify both eigenproblem residuals and unit score variances.

Check Your Understanding

Use the definitions and checks above to answer these questions.

What are the dimensions of $S_{XY}$ and $S_{XY}S_{YY}^{-1}S_{YX}$ when the two sets have $p$ and $q$ variables?
If a generalized eigenvalue is $0.36$, what is the corresponding canonical correlation?
Why must the $X$- and $Y$-observations be paired in the same order when calculating $S_{XY}$?
What two score-variance checks should be made after recovering $a$ and $b$?
Explain why $a^\top S_{XY}b$ need not equal the canonical correlation before both directions are normalized.

Computing Canonical Roots From Sample Covariances

What You'll Learn