Principal components by Lagrange multipliers
Section of Principal component analysis, part of Statistical learning
Principal components by Lagrange multipliers
| Type | Section |
| Taught in | Principal component analysis |
| Part | Statistical learning |
Let \mathbf{x}_1,\ldots,\mathbf{x}_n\in\mathbb{R}^p have sample mean \bar{\mathbf{x}} and sample covariance S. The first direction is the coefficient vector of a linear combination whose sample variance is as large as possible. The combination is
Y_1=a_{11}X_1+\cdots+a_{p1}X_p=\mathbf{a}_1^\top\mathbf{X}.
The sample version of \operatorname{Var}(\mathbf{a}^\top\mathbf{X})=\mathbf{a}^\top\Sigma\mathbf{a}, with S in place of \Sigma, writes that variance as a quadratic form,
s^2_{Y_1}=\sum_{j=1}^p\sum_{k=1}^p a_{j1}a_{k1}s_{jk}=\mathbf{a}_1^\top S\mathbf{a}_1.
If there is a direction with positive variance, the unconstrained optimisation has no finite answer. Replace \mathbf{a}_1 by c\mathbf{a}_1 and the variance becomes c^2 times as large. That is a change of ruler, not a better summary. Fix \mathbf{a}_1^\top\mathbf{a}_1=1 so that orientations can be compared. A Lagrange multiplier enforces the constraint:
u_1=\mathbf{a}_1^\top S\mathbf{a}_1-\lambda(\mathbf{a}_1^\top\mathbf{a}_1-1).
Differentiate in the direction. For symmetric S the gradient of \mathbf{a}^\top S\mathbf{a} is 2S\mathbf{a}, so
\frac{\partial u_1}{\partial\mathbf{a}_1}=2S\mathbf{a}_1-2\lambda\mathbf{a}_1.
Set the gradient to zero:
S\mathbf{a}_1=\lambda\mathbf{a}_1.
The stationary points are the eigenvectors of S, and the multipliers are the eigenvalues. Substitute a stationary point back into the objective:
\mathbf{a}_1^\top S\mathbf{a}_1=\lambda\,\mathbf{a}_1^\top\mathbf{a}_1=\lambda.
At a unit eigenvector, the objective equals its eigenvalue. The objective is continuous and the unit sphere is compact, so a maximum exists. The largest eigenvalue gives it. If that eigenvalue is repeated, several directions are equally good. Even with a unique largest eigenvalue, \mathbf{a}_1 and -\mathbf{a}_1 describe the same axis.
A second direction that is not constrained will select the same eigenvector again. Later directions are required to be orthogonal to those already chosen. One construction removes the first component from the centred data and repeats the maximisation on what remains. The theorem reaches the same directions in one induction. Any unit vector is written in the eigenbasis of S, so its sample variance becomes a weighted average of the eigenvalues. Orthogonality to the first k-1 eigenvectors sets the first k-1 weights to zero, and the largest remaining eigenvalue is then the maximum.
Theorem 1 (Principal component directions) Let S=\Gamma\Lambda\Gamma^\top with \lambda_1\ge\cdots\ge\lambda_p\ge 0 and \Gamma=(\boldsymbol{\gamma}_1,\ldots,\boldsymbol{\gamma}_p) orthogonal. The jth principal component direction is \mathbf{a}_j=\boldsymbol{\gamma}_j. The jth principal component Y_j=\boldsymbol{\gamma}_j^\top\mathbf{X} has sample variance \lambda_j. Distinct components are uncorrelated, and the total variance is preserved:
\sum_{j=1}^p\operatorname{Var}(Y_j)=\sum_{j=1}^p\lambda_j=\operatorname{tr}(S)=\sum_{j=1}^p s_{jj}.
Any unit vector can be written \mathbf{a}=\Gamma\mathbf{b} with \mathbf{b}=\Gamma^\top\mathbf{a} and \mathbf{b}^\top\mathbf{b}=1. Then
\mathbf{a}^\top S\mathbf{a}=\mathbf{b}^\top\Gamma^\top\Gamma\Lambda\Gamma^\top\Gamma\mathbf{b}=\sum_{j=1}^p\lambda_j b_j^2.
For j=1 the right-hand side is a weighted average of the \lambda_j with weights b_j^2 that sum to one, so it is at most \lambda_1, with equality for \mathbf{b}=\mathbf{e}_1, that is \mathbf{a}=\boldsymbol{\gamma}_1. Suppose \mathbf{a}_1=\boldsymbol{\gamma}_1,\ldots,\mathbf{a}_{k-1}=\boldsymbol{\gamma}_{k-1}. Orthogonality of \mathbf{a} to these directions means b_1=\cdots=b_{k-1}=0, so
\mathbf{a}^\top S\mathbf{a}=\sum_{j\ge k}\lambda_j b_j^2\le\lambda_k,
again attained at \mathbf{b}=\mathbf{e}_k. By induction \mathbf{a}_j=\boldsymbol{\gamma}_j for every j. The variance of Y_j is \boldsymbol{\gamma}_j^\top S\boldsymbol{\gamma}_j=\lambda_j. For j\neq k,
\operatorname{Cov}(Y_j,Y_k)=\boldsymbol{\gamma}_j^\top S\boldsymbol{\gamma}_k=\lambda_k\boldsymbol{\gamma}_j^\top\boldsymbol{\gamma}_k=0.
Finally \sum_j\lambda_j=\operatorname{tr}(\Lambda)=\operatorname{tr}(\Gamma^\top S\Gamma)=\operatorname{tr}(S\Gamma\Gamma^\top)=\operatorname{tr}(S), where the trace cycles.
The calculation has diagonalised the sample covariance. It guarantees zero sample covariance between different score columns. It does not give independence, and it does not give zero covariance for arbitrary future observations. In the theorem, \operatorname{Var} and \operatorname{Cov} are the empirical score variances and covariances computed with divisor n-1. For a population analysis, repeat the argument with \Sigma in place of S. Under an additional Gaussian population model, the population principal components are independent, because they are uncorrelated linear images of a normal vector. Without normality, uncorrelated is the conclusion the argument has earned.
The proportion of variance explained by the first k components is
\frac{\sum_{j\le k}\lambda_j}{\sum_{j\le p}\lambda_j}.
Reconstruction
Main article: Reconstruction from principal components
The sum of all component variances equals \operatorname{tr}(S), which does not depend on the basis. Maximising the variance kept by the first k components is therefore the same problem as minimising the variance left out. Let \Gamma_k=(\boldsymbol{\gamma}_1,\ldots,\boldsymbol{\gamma}_k). Replacing each centred observation by its projection \Gamma_k\Gamma_k^\top(\mathbf{x}_i-\bar{\mathbf{x}}) leaves residuals whose total squared length, divided by n-1, equals
\sum_{j=k+1}^p\lambda_j.
That value is the smallest residual attainable by any rank-k linear approximation. The leading components are at once the highest-variance coordinates and the best linear summary of the cloud in squared error. A later autoencoder network generalises the second description: a hidden layer smaller than the input tries to reproduce the input after passing through a lower-dimensional representation. Factor analysis and non-negative matrix factorisation are further representations, with different constraints.
Used in
See also
Principal component directions and scores · Principal component directions · Scaling · Reconstruction from principal components · Choosing the number of principal components · How many components, on a correlation with no response · Singular value decomposition · Iris, with and without scaling