Unsupervised Learning II — PCA, SVD & Factor Analysis

Compressing the cross-section of returns into its underlying factors

Clustering grouped the observations; dimensionality reduction compresses the features. Principal component analysis (PCA) finds the orthogonal directions of maximum variance in the data — a new, ranked coordinate system in which the first few axes capture most of the variation, so a high-dimensional dataset can be summarised by a handful of numbers. The clean way to compute it is the singular value decomposition of the centered data, $X = U S V^\top$: the columns of $V$ are the principal directions, $US$ are the coordinates, and $S_k^2$ is the variance along axis $k$.

The finance application is one of the most useful facts in quantitative finance: the returns of many assets are driven by a few common factors. Run PCA on a cross-section of stock returns and the first component is essentially the market; the next few are style/sector contrasts; the rest is idiosyncratic noise. That low-dimensional factor structure underlies index models, risk decomposition, and the covariance-shrinkage of the asset-risk arc.

We build PCA from scratch via SVD (validated against scikit-learn), apply it to 48 stocks to recover the market factor and principal portfolios, use it to reveal the low-rank structure of the covariance matrix, and contrast it with factor analysis — its probabilistic cousin that separates common from idiosyncratic variance, the frequentist sibling of your Bayesian factor models. Python-lead.

1. The data — a cross-section of returns

Weekly returns of 48 large-cap US stocks (the asset-risk arc's stocks_weekly.csv, 2019–2024). Stacked as a 312×48 matrix, the striking feature is co-movement: the correlation heatmap is overwhelmingly positive — when the market moves, most stocks move together. That shared variation is what PCA will distil into a few factors. We standardise each stock (unit variance) so the analysis is about correlation structure, not which stock happens to be most volatile.

import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings
warnings.filterwarnings("ignore")
from pca import pca, reconstruct
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
d=pd.read_csv("stocks_weekly.csv",index_col=0); R=d.values/100.0; tickers=list(d.columns); N=R.shape[1]
dates=pd.to_datetime(d.index)
print(f"{R.shape[0]} weeks x {N} stocks;  mean pairwise correlation = {np.corrcoef(R.T)[np.triu_indices(N,1)].mean():.2f}")
fig,ax=plt.subplots(figsize=(6.2,5.2)); im=ax.imshow(np.corrcoef(R.T),cmap="RdBu_r",vmin=-1,vmax=1)
ax.set_title("Stock return correlations — overwhelmingly positive co-movement"); ax.set_xticks([]); ax.set_yticks([])
plt.colorbar(im,ax=ax,fraction=0.046,label="correlation"); plt.tight_layout(); plt.show()
print("Almost every pair is positively correlated -- a shared driver (the market) plus finer structure. PCA extracts them.")

312 weeks x 48 stocks;  mean pairwise correlation = 0.38

Almost every pair is positively correlated -- a shared driver (the market) plus finer structure. PCA extracts them.

2. PCA from scratch via SVD

pca.py centers (and here standardises) the data and takes its SVD; the variance explained by each component is $S_k^2/(n-1)$. We confirm the explained-variance ratios match scikit-learn, then plot the scree: the first component alone captures a large share of all the co-movement, and the curve falls off fast — the hallmark of low-dimensional factor structure.

from sklearn.decomposition import PCA
Rz=(R-R.mean(0))/R.std(0,ddof=1)
res=pca(R,standardize=True); sk=PCA().fit(Rz)
print("explained-variance ratio (first 5):")
print("  from-scratch:", np.round(res['explained_variance_ratio'][:5],4))
print("  scikit-learn:", np.round(sk.explained_variance_ratio_[:5],4))
evr=res['explained_variance_ratio']; cum=np.cumsum(evr)
print(f"\nPC1 explains {evr[0]:.1%};  first 3 PCs {cum[2]:.1%};  first 10 {cum[9]:.1%}  of the {N}-stock covariation")
fig,ax=plt.subplots(1,2,figsize=(13,4))
ax[0].bar(range(1,16),evr[:15]*100,color=BLUE); ax[0].set_xlabel("principal component"); ax[0].set_ylabel("variance explained (%)"); ax[0].set_title("Scree — PC1 dominates")
ax[1].plot(range(1,N+1),cum*100,"o-",color=GREEN,lw=2,ms=3); ax[1].axhline(90,color=RED,ls=":"); ax[1].set_xlabel("number of components"); ax[1].set_ylabel("cumulative variance (%)"); ax[1].set_title(f"{np.argmax(cum>=0.9)+1} PCs explain 90% of {N} stocks")
plt.tight_layout(); plt.show()
print("From-scratch SVD matches scikit-learn exactly. One component captures ~40% of all co-movement across 48 stocks --")
print("the market. The effective dimension of the cross-section is a fraction of its 48 nominal dimensions.")

explained-variance ratio (first 5):
  from-scratch: [0.4002 0.0925 0.0694 0.0329 0.03  ]
  scikit-learn: [0.4002 0.0925 0.0694 0.0329 0.03  ]

PC1 explains 40.0%;  first 3 PCs 56.2%;  first 10 73.1%  of the 48-stock covariation

From-scratch SVD matches scikit-learn exactly. One component captures ~40% of all co-movement across 48 stocks --
the market. The effective dimension of the cross-section is a fraction of its 48 nominal dimensions.

3. The market factor and principal portfolios

Each principal direction is a portfolio — a set of weights on the 48 stocks — and its score is that portfolio's return. Reading them off:

• PC1 is the market. Its loadings are (up to an arbitrary overall sign) all the same sign — a long-everything portfolio — and its return correlates ~1.0 with the equal-weight market index. PCA rediscovers the market factor from the data alone, with no index supplied.
• Higher components are long/short style factors. PC2, PC3… have mixed-sign loadings — they go long one group of stocks and short another, capturing sector or style contrasts orthogonal to the market.

This is the statistical foundation of factor investing: the CAPM/Fama-French factors are hypothesised portfolios; PCA discovers the dominant ones empirically.

l1=res['loadings'][0].copy(); s1=res['scores'][:,0].copy()
if l1.mean()<0: l1,s1=-l1,-s1                                    # fix arbitrary sign so the market factor is positive
mkt=Rz.mean(1)
print(f"PC1 loadings all positive: {np.mean(l1>0)*100:.0f}%   corr(PC1 return, equal-weight market) = {abs(np.corrcoef(s1,mkt)[0,1]):.3f}")
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
o=np.argsort(l1); ax[0].bar(range(N),l1[o],color=BLUE); ax[0].set_title("PC1 loadings — all one sign = the MARKET factor"); ax[0].set_xlabel("stock (sorted)"); ax[0].set_ylabel("loading"); ax[0].axhline(0,color="k",lw=.5)
l2=res['loadings'][1]; o2=np.argsort(l2)
ax[1].bar(range(N),l2[o2],color=[GREEN if v>0 else RED for v in l2[o2]]); ax[1].set_title("PC2 loadings — mixed sign = a long/short style factor"); ax[1].set_xlabel("stock (sorted)"); ax[1].set_ylabel("loading"); ax[1].axhline(0,color="k",lw=.5)
plt.tight_layout(); plt.show()
top=np.argsort(l2)[::-1][:3]; bot=np.argsort(l2)[:3]
print(f"PC2 goes long {', '.join(tickers[i] for i in top)} and short {', '.join(tickers[i] for i in bot)} -- a style/sector contrast.")
print("PC1 is a near-perfect proxy for the market; the higher PCs are orthogonal long/short bets -- principal portfolios.")

PC1 loadings all positive: 100%   corr(PC1 return, equal-weight market) = 0.998

PC2 goes long CVX, XOM, COP and short AMZN, ADBE, NVDA -- a style/sector contrast.
PC1 is a near-perfect proxy for the market; the higher PCs are orthogonal long/short bets -- principal portfolios.

4. Low-rank structure — reconstructing the covariance

Because a few components carry most of the variance, a low-rank approximation of the correlation matrix recovers most of its structure. The right construction is a spectral truncation — keep the top $k$ eigenvalue/eigenvector pairs of $C$ itself, $C_k=\sum_{j\le k}\lambda_j v_jv_j^\top$, which is its best rank-$k$ approximation in Frobenius norm. Below, the full 48×48 matrix and its rank-3 truncation next to it, with the error curve reported both for the whole matrix and for the off-diagonal entries alone (the diagonal is 1 by construction and flatters any reconstruction).

The cell also shows why the tempting shortcut — rebuild the data at rank $k$ and correlate it — answers a different question and cannot be read as reconstruction quality. This is the same fact the asset-risk arc exploited from the other side: the sample covariance is noisy and near-singular, and Ledoit–Wolf shrinkage pulls it toward exactly this low-rank/structured target. It is also the factor-zoo lesson again — many assets, few independent dimensions.

full=np.corrcoef(Rz.T)
# The object to approximate is the CORRELATION MATRIX itself, so truncate its spectrum:
# C_k = sum_{j<=k} lambda_j v_j v_j', the best rank-k approximation of C in Frobenius norm.
_w,_V=np.linalg.eigh(full); _o=np.argsort(_w)[::-1]; _w,_V=_w[_o],_V[:,_o]
def corr_from_k(k): return (_V[:,:k]*_w[:k])@_V[:,:k].T
# It is worth seeing why the tempting shortcut does not work: rebuild the DATA at rank k and correlate it.
def corr_of_lowrank_data(k): return np.corrcoef(reconstruct(res,k).T)

iu=np.triu_indices(N,1)
r3=corr_from_k(3)
fig,ax=plt.subplots(1,3,figsize=(15,4.2))
for a,(M,t) in zip(ax[:2],[(full,"full correlation matrix"),(r3,"rank-3 spectral reconstruction")]):
    im=a.imshow(M,cmap="RdBu_r",vmin=-1,vmax=1); a.set_title(t); a.set_xticks([]); a.set_yticks([]); plt.colorbar(im,ax=a,fraction=0.046)
err=[np.linalg.norm(full-corr_from_k(k),"fro")/np.linalg.norm(full,"fro") for k in range(1,21)]
erro=[np.linalg.norm(full[iu]-corr_from_k(k)[iu])/np.linalg.norm(full[iu]) for k in range(1,21)]
ax[2].plot(range(1,21),err,"o-",color=PURP,lw=2,ms=3,label="whole matrix")
ax[2].plot(range(1,21),erro,"o-",color=ORANGE,lw=2,ms=3,label="off-diagonal only")
ax[2].set_xlabel("components used"); ax[2].set_ylabel("relative reconstruction error")
ax[2].set_title("A few factors rebuild the correlation matrix"); ax[2].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"Rank-3 spectral reconstruction: {err[2]:.1%} relative error on the whole matrix, {erro[2]:.1%} on the off-diagonal")
print(f"entries alone (the diagonal is 1 by construction and flatters any reconstruction). Rank-10 brings those to")
print(f"{err[9]:.1%} and {erro[9]:.1%}. Real low-rank structure, then -- three numbers per stock carry most of a 48x48")
print("correlation matrix -- but 'almost perfectly' would be overselling a 12% off-diagonal error.")
print("\nOne trap worth naming, because the obvious alternative is wrong. It is tempting to rebuild the DATA at rank k")
print("and take its correlation matrix. That measures something else entirely:")
for _k in (1,2,3):
    _M=corr_of_lowrank_data(_k)
    print(f"   rank {_k}: mean |correlation| of the rebuilt data {np.abs(_M[iu]).mean():.3f}   vs {np.abs(full[iu]).mean():.3f} in the real data")
print("A rank-k matrix has only k independent directions, and correlating it rescales every column by its OWN")
print("reconstructed standard deviation -- so the correlations are driven toward +/-1 no matter how good the fit is.")
print("At rank 1 every pair reads exactly 1.000. The spectral truncation above is the quantity actually being claimed.")
print("\nThis low-rank structure is what Ledoit-Wolf shrinkage targets in the high-dimensional portfolio notebook -- PCA")
print("shows it directly; shrinkage uses it to stabilize the noisy sample covariance.")

Rank-3 spectral reconstruction: 19.5% relative error on the whole matrix, 12.2% on the off-diagonal
entries alone (the diagonal is 1 by construction and flatters any reconstruction). Rank-10 brings those to
11.9% and 7.4%. Real low-rank structure, then -- three numbers per stock carry most of a 48x48
correlation matrix -- but 'almost perfectly' would be overselling a 12% off-diagonal error.

One trap worth naming, because the obvious alternative is wrong. It is tempting to rebuild the DATA at rank k
and take its correlation matrix. That measures something else entirely:
   rank 1: mean |correlation| of the rebuilt data 1.000   vs 0.375 in the real data
   rank 2: mean |correlation| of the rebuilt data 0.807   vs 0.375 in the real data
   rank 3: mean |correlation| of the rebuilt data 0.699   vs 0.375 in the real data
A rank-k matrix has only k independent directions, and correlating it rescales every column by its OWN
reconstructed standard deviation -- so the correlations are driven toward +/-1 no matter how good the fit is.
At rank 1 every pair reads exactly 1.000. The spectral truncation above is the quantity actually being claimed.

This low-rank structure is what Ledoit-Wolf shrinkage targets in the high-dimensional portfolio notebook -- PCA
shows it directly; shrinkage uses it to stabilize the noisy sample covariance.

5. PCA vs factor analysis — variance vs latent factors

PCA and factor analysis (FA) are often confused but answer different questions. PCA finds directions of maximum total variance, mixing common co-movement with each stock's idiosyncratic noise. FA is an explicit probabilistic model, $X = Lf + \varepsilon$: a few common latent factors $f$ plus per-stock unique noise $\varepsilon$, and it separates the two — estimating how much of each stock's variance is shared (communality) vs idiosyncratic (uniqueness). For returns that distinction is meaningful: the shared part is systematic (undiversifiable) risk, the unique part is diversifiable. FA is the frequentist face of the Bayesian factor models in your IRT and factor-VAR arcs; PCA is the quick, assumption-light summary.

from sklearn.decomposition import FactorAnalysis
fa=FactorAnalysis(n_components=3,random_state=0).fit(Rz)
uniq=fa.noise_variance_                                          # idiosyncratic (diversifiable) variance per stock
commun=1-uniq                                                    # communality (systematic share), std'd data -> total var 1
fig,ax=plt.subplots(1,2,figsize=(13,4.2))
o=np.argsort(commun)
ax[0].bar(range(N),commun[o],color=GREEN,label="communality (systematic)"); ax[0].bar(range(N),uniq[o],bottom=commun[o],color=GREY,label="uniqueness (idiosyncratic)")
ax[0].set_title("Factor analysis splits each stock's variance"); ax[0].set_xlabel("stock (sorted by communality)"); ax[0].set_ylabel("variance share"); ax[0].legend(fontsize=8)
# compare FA loading-1 to PCA loading-1
fl=fa.components_[0].copy();
if np.corrcoef(fl,l1)[0,1]<0: fl=-fl
ax[1].scatter(l1,fl,s=18,color=BLUE); ax[1].set_xlabel("PCA PC1 loading"); ax[1].set_ylabel("FA factor-1 loading"); ax[1].set_title(f"PCA vs FA first factor (corr {abs(np.corrcoef(l1,fl)[0,1]):.2f})")
plt.tight_layout(); plt.show()
print(f"Average systematic (undiversifiable) share across stocks: {commun.mean():.0%}; the rest is idiosyncratic.")
print("PCA's PC1 and FA's first factor agree closely on the market direction, but FA additionally quantifies how much of")
print("each stock is systematic vs diversifiable -- the split that matters for risk, and the object Bayesian factor models put a posterior on.")

Average systematic (undiversifiable) share across stocks: 53%; the rest is idiosyncratic.
PCA's PC1 and FA's first factor agree closely on the market direction, but FA additionally quantifies how much of
each stock is systematic vs diversifiable -- the split that matters for risk, and the object Bayesian factor models put a posterior on.

6. Summary

PCA compresses correlated features into a few variance-maximising directions; on a cross-section of returns those directions are the market and style factors. Built from scratch via SVD (matched to scikit-learn), it showed the defining fact of quantitative finance: 48 stocks are driven by a handful of components — PC1 is the market (≈40% of variance, all-positive loadings, ~1.0 correlation with the index), higher PCs are long/short style portfolios, and a rank-3 spectral truncation reproduces the correlation matrix to about 19% relative error, 12% off the diagonal. Factor analysis then split each stock's variance into a systematic (undiversifiable) and an idiosyncratic (diversifiable) part — the distinction PCA blurs.

Cross-links run throughout the portfolio: the low-rank structure PCA reveals is exactly what Ledoit–Wolf shrinkage targets in the high-dimensional-portfolio notebook; the "few independent dimensions" echo the factor-zoo result; and FA is the frequentist sibling of the Bayesian factor models in the factor-VAR and IRT arcs.

Next: autoencoders — a neural network that learns a nonlinear compression, generalising PCA (a linear autoencoder with squared-error loss recovers PCA exactly), with an anomaly-detection application.