Unsupervised Learning III — Autoencoders

Neural compression: linear autoencoders recover PCA, nonlinear ones surpass it

An autoencoder is a neural network trained to copy its input to its output through a bottleneck. An encoder compresses the input to a low-dimensional latent code; a decoder reconstructs the input from that code; and the network minimises the reconstruction error. Because information has to squeeze through the narrow bottleneck, the network is forced to learn an efficient, compressed representation — dimensionality reduction, learned by backprop.

The link to the previous notebook is exact: a linear autoencoder with squared-error loss learns the same subspace as PCA. The autoencoder's power comes from making the encoder and decoder nonlinear — then it can capture curved structure that PCA's straight axes cannot, compressing the same information into fewer dimensions. We show that ladder — linear AE ≈ PCA, then nonlinear AE > PCA on images — visualise the learned latent space, and put the reconstruction error to work as a finance anomaly detector. Python-only (deep learning).

1. Linear autoencoder = PCA

Start with the tightest possible link. A linear autoencoder — encoder $x\mapsto W_e x$, decoder $z\mapsto W_d z$, no activations, squared-error loss — is provably equivalent to PCA: at the optimum it spans the same top-$k$ principal subspace, so its reconstruction error equals PCA's. We train a 2-unit-bottleneck linear autoencoder on the 48-stock returns from the PCA notebook and confirm its reconstruction error lands on the PCA rank-2 value.

It is worth being precise about what the theorem does and does not say, because the two are easy to run together. The subspace is the same; the axes inside it are not. PCA returns ordered, orthogonal components; the autoencoder returns whatever basis of that subspace gradient descent happened to land on, and nothing in the squared-error loss prefers one basis over another. Both facts are checked below. This is the anchor: whatever an autoencoder adds beyond PCA comes entirely from nonlinearity.

import os
os.environ["KMP_DUPLICATE_LIB_OK"]="TRUE"
import numpy as np, pandas as pd, matplotlib.pyplot as plt, time, warnings
warnings.filterwarnings("ignore")
import torch, torch.nn as nn; torch.set_num_threads(2)
from scipy.stats import spearmanr
from sklearn.decomposition import PCA
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
d=pd.read_csv("stocks_weekly.csv",index_col=0); dates=pd.to_datetime(d.index)
R=d.values/100.0; R=(R-R.mean(0))/R.std(0); N=R.shape[1]
pca2=PCA(2).fit(R); pca_err=np.mean((R-pca2.inverse_transform(pca2.transform(R)))**2)
Xt=torch.tensor(R,dtype=torch.float32); torch.manual_seed(0)
lae=nn.Sequential(nn.Linear(N,2),nn.Linear(2,N)); opt=torch.optim.Adam(lae.parameters(),1e-2)
for e in range(2000): opt.zero_grad(); ((lae(Xt)-Xt)**2).mean().backward(); opt.step()
with torch.no_grad(): lae_err=float(((lae(Xt)-Xt)**2).mean())
print(f"48-stock returns, 2-D bottleneck:")
print(f"  PCA rank-2 reconstruction MSE : {pca_err:.6f}")
print(f"  linear autoencoder MSE        : {lae_err:.6f}   -> agree to {abs(pca_err-lae_err)/pca_err:.1e} relative")

# What is identical is the SUBSPACE. The axes inside it are not PCA's, and there is no reason they should be.
from scipy.linalg import subspace_angles
We=lae[0].weight.detach().numpy(); Vp=pca2.components_
_ang=np.degrees(subspace_angles(We.T,Vp.T))
_u=lambda A: A/np.linalg.norm(A,axis=1,keepdims=True)
_c=np.abs(_u(We)@_u(Vp).T); _g=np.abs(_u(We)@_u(We).T)
print(f"\nPrincipal angles between the encoder's subspace and PCA's top-2 subspace: {_ang[0]:.3f} and {_ang[1]:.3f} degrees.")
print("The theorem is about that subspace, not about the axes -- and the axes are visibly not PCA's:")
print(f"  |cos| of encoder row 1 with (PC1, PC2): {_c[0,0]:.2f}, {_c[0,1]:.2f}")
print(f"  |cos| of encoder row 2 with (PC1, PC2): {_c[1,0]:.2f}, {_c[1,1]:.2f}")
print(f"  the two encoder rows are not even orthogonal to each other (|cos| = {_g[0,1]:.2f}).")
print("A linear autoencoder recovers PCA's subspace, so it reconstructs identically; it does not recover PCA's ordered,")
print("orthogonal components, and nothing in the loss asks it to. Any basis of the right subspace is an optimum.")
print("\nSo an autoencoder only earns its keep once encoder/decoder are NONLINEAR -- which is the rest of this notebook.")

48-stock returns, 2-D bottleneck:
  PCA rank-2 reconstruction MSE : 0.507306
  linear autoencoder MSE        : 0.507306   -> agree to 5.9e-08 relative

Principal angles between the encoder's subspace and PCA's top-2 subspace: 0.000 and 0.000 degrees.
The theorem is about that subspace, not about the axes -- and the axes are visibly not PCA's:
  |cos| of encoder row 1 with (PC1, PC2): 0.66, 0.75
  |cos| of encoder row 2 with (PC1, PC2): 0.73, 0.69
  the two encoder rows are not even orthogonal to each other (|cos| = 0.04).
A linear autoencoder recovers PCA's subspace, so it reconstructs identically; it does not recover PCA's ordered,
orthogonal components, and nothing in the loss asks it to. Any basis of the right subspace is an optimum.

So an autoencoder only earns its keep once encoder/decoder are NONLINEAR -- which is the rest of this notebook.

2. Nonlinear autoencoder on images

Now the nonlinear version, on Fashion-MNIST (28×28 = 784 pixels). The encoder 784→128→32 (ReLU) squeezes each image to a 32-number code; the decoder mirrors it back to 784 pixels. Trained to minimise pixel reconstruction error, it learns a nonlinear 32-dimensional manifold of clothing images. At the same bottleneck size, the nonlinear autoencoder reconstructs more faithfully than PCA-32, because it can bend around the image manifold rather than fit a flat subspace. The grid compares originals, PCA-32, and autoencoder reconstructions.

import torchvision, torchvision.transforms as T
tr=torchvision.datasets.FashionMNIST("data",train=True,download=True,transform=T.ToTensor())
te=torchvision.datasets.FashionMNIST("data",train=False,download=True,transform=T.ToTensor())
Xtr=tr.data.view(-1,784).float().numpy()/255; Xte=te.data.view(-1,784).float().numpy()/255; yte=te.targets.numpy(); classes=tr.classes
p32=PCA(32).fit(Xtr); pca_mse=np.mean((Xte-p32.inverse_transform(p32.transform(Xte)))**2)
class AE(nn.Module):
    def __init__(s,b=32):
        super().__init__(); s.enc=nn.Sequential(nn.Linear(784,128),nn.ReLU(),nn.Linear(128,b))
        s.dec=nn.Sequential(nn.Linear(b,128),nn.ReLU(),nn.Linear(128,784),nn.Sigmoid())
    def forward(s,x): return s.dec(s.enc(x))
torch.manual_seed(0); ae=AE(32); opt=torch.optim.Adam(ae.parameters(),1e-3); lf=nn.MSELoss()
Xtrt=torch.tensor(Xtr,dtype=torch.float32); t=time.time()
for ep in range(10):
    pm=torch.randperm(len(Xtrt))
    for b in range(0,len(Xtrt),256): bi=pm[b:b+256]; opt.zero_grad(); lf(ae(Xtrt[bi]),Xtrt[bi]).backward(); opt.step()
Xtet=torch.tensor(Xte,dtype=torch.float32)
with torch.no_grad(): ae_rec=ae(Xtet).numpy(); ae_mse=np.mean((Xte-ae_rec)**2)
pca_rec=p32.inverse_transform(p32.transform(Xte))
print(f"FashionMNIST, bottleneck=32:  PCA MSE {pca_mse:.4f}  |  nonlinear AE MSE {ae_mse:.4f}   ({time.time()-t:.0f}s)")
idx=[np.where(yte==c)[0][0] for c in range(8)]
fig,ax=plt.subplots(3,8,figsize=(13,4.6))
for j,i in enumerate(idx):
    ax[0,j].imshow(Xte[i].reshape(28,28),cmap="gray"); ax[1,j].imshow(pca_rec[i].reshape(28,28),cmap="gray"); ax[2,j].imshow(ae_rec[i].reshape(28,28),cmap="gray")
    for r in range(3): ax[r,j].axis("off")
    ax[0,j].set_title(classes[yte[i]],fontsize=7)
for r,lab in enumerate(["original","PCA-32","autoencoder-32"]): ax[r,0].set_ylabel(lab,fontsize=9,rotation=90); ax[r,0].axis("on"); ax[r,0].set_xticks([]); ax[r,0].set_yticks([])
plt.tight_layout(); plt.show()
print(f"At the same 32-D bottleneck the nonlinear AE reconstructs ~{100*(pca_mse-ae_mse)/pca_mse:.0f}% more faithfully than PCA --")
print("sharper edges and cleaner shapes -- because it bends to the nonlinear manifold of clothing images.")

FashionMNIST, bottleneck=32:  PCA MSE 0.0151  |  nonlinear AE MSE 0.0130   (6s)

At the same 32-D bottleneck the nonlinear AE reconstructs ~14% more faithfully than PCA --
sharper edges and cleaner shapes -- because it bends to the nonlinear manifold of clothing images.

3. The learned latent space

Squeeze the bottleneck to 2 dimensions and the autoencoder learns a map we can plot. Encoding the test set and colouring by (held-out) class label shows the classes separating into regions — footwear here, bags there, upper-body garments clustered together — discovered with no labels used in training. Against PCA's 2 linear components, the nonlinear autoencoder pulls the classes apart more cleanly. (The next notebook, t-SNE/UMAP, specialises in exactly this 2-D visualisation.)

torch.manual_seed(0); ae2=AE(2); o2=torch.optim.Adam(ae2.parameters(),1e-3)
for ep in range(10):
    pm=torch.randperm(len(Xtrt))
    for b in range(0,len(Xtrt),256): bi=pm[b:b+256]; o2.zero_grad(); lf(ae2(Xtrt[bi]),Xtrt[bi]).backward(); o2.step()
with torch.no_grad(): Zae=ae2.enc(Xtet).numpy()
Zpca=PCA(2).fit(Xtr).transform(Xte)
s=np.random.default_rng(0).choice(len(yte),3000,replace=False)
fig,ax=plt.subplots(1,2,figsize=(13,5))
for a,(Zp,t_) in zip(ax,[(Zpca,"PCA (2 linear components)"),(Zae,"Autoencoder (2-D nonlinear latent)")]):
    sc=a.scatter(Zp[s,0],Zp[s,1],c=yte[s],cmap="tab10",s=6,alpha=.6); a.set_title(t_); a.set_xticks([]); a.set_yticks([])
cb=fig.colorbar(sc,ax=ax,fraction=0.025,ticks=range(10)); cb.ax.set_yticklabels(classes,fontsize=7)
plt.show()
# "separates more cleanly" is an eyeball claim unless it is scored, so score it: how well can the held-out
# class label be read off each 2-D embedding, and how well separated are the classes within it?
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import cross_val_score
from sklearn.metrics import silhouette_score
from sklearn.preprocessing import StandardScaler
print("Both embed 784-D images into 2-D. Scoring the separation rather than eyeballing it:")
for _nm,_Z in [("PCA-2",Zpca),("AE-2",Zae)]:
    _Zs=StandardScaler().fit_transform(_Z)
    _a=cross_val_score(KNeighborsClassifier(15),_Zs[s],yte[s],cv=5).mean(); _si=silhouette_score(_Zs[s],yte[s])
    print(f"   {_nm:6s}: 15-NN accuracy on the class label {_a:.3f}   silhouette by true class {_si:+.3f}")
print("   (10 classes, so 0.100 is chance. Neither map ever saw a label -- the labels only score the result.)")
print("The nonlinear embedding is the better one, and by a clear margin rather than a matter of taste.")

Both embed 784-D images into 2-D. Scoring the separation rather than eyeballing it:
   PCA-2 : 15-NN accuracy on the class label 0.525   silhouette by true class -0.026

   AE-2  : 15-NN accuracy on the class label 0.673   silhouette by true class +0.080
   (10 classes, so 0.100 is chance. Neither map ever saw a label -- the labels only score the result.)
The nonlinear embedding is the better one, and by a clear margin rather than a matter of taste.

4. Anomaly detection — reconstruction error as a novelty score

An autoencoder trained on "normal" data reconstructs normal inputs well and unusual inputs badly — so the reconstruction error is an unsupervised anomaly score. We train one on the cross-section of 48 stock returns (each week is a 48-dimensional vector) and score each week by how poorly it reconstructs.

What the two plots below show, since a dated line chart on a finance page is usually one asset or an index and this is neither. Every week is a 48-number vector — that week's return for each stock. The network squeezes those 48 numbers through a 4-dimensional bottleneck and rebuilds them, and the plotted value is the squared rebuilding error averaged across all 48 stocks for that week (((net(X)-X)**2).mean(1), where the .mean(1) is the average over stocks). So a single point on the line answers: how unlike the usual pattern was the shape of this week's cross-section?

That is a different question from “how much did the market move”, and the difference is worth holding on to when reading the peaks. A week in which all 48 stocks fall together can score low: a common market move is one direction, and one direction is the easiest thing for a 4-dimensional bottleneck to capture — it is essentially PC1, which the PCA notebook shows carries 40% of the covariation on its own. What scores high is a week whose internal arrangement is unusual: dispersion far outside the normal range, or sectors moving against each other in a combination the network has not learned to express in four numbers. The score is about the shape of the cross-section, not its average.

The method has a trap in it, and it is worth walking into deliberately because it is easy to miss and it reverses the answer. The obvious implementation trains on every week and then scores those same weeks. But this network has 1,748 parameters and 312 weeks to fit, and a network with capacity to spare does not learn "normal" — it learns these observations, the extreme ones included. An outlier that the network has memorised reconstructs beautifully, so it scores as maximally normal. The anomaly ranking inverts exactly where it is supposed to work.

The fix is the same one supervised learning uses without thinking about it: score data the model was not fitted on. We compute both — the in-sample score and a 6-fold out-of-sample score, where every week is reconstructed by a network trained without it — and compare them against a linear baseline the method should be made to beat, PCA at the same rank.

def anom_net(seed=0):
    torch.manual_seed(seed)
    return nn.Sequential(nn.Linear(N,16),nn.ReLU(),nn.Linear(16,4),nn.ReLU(),nn.Linear(4,16),nn.ReLU(),nn.Linear(16,N))
def fit(net,Xtrain,steps=1500):
    o=torch.optim.Adam(net.parameters(),5e-3)
    for _ in range(steps): o.zero_grad(); ((net(Xtrain)-Xtrain)**2).mean().backward(); o.step()
    return net
npar=sum(p.numel() for p in anom_net().parameters())

net=fit(anom_net(),Xt)                                          # the in-sample version
with torch.no_grad(): rec_in=((net(Xt)-Xt)**2).mean(1).numpy()

folds=np.random.default_rng(0).permutation(len(Xt))%6            # every week scored by a net trained without it
rec_oos=np.zeros(len(Xt))
for f in range(6):
    tr_i=np.where(folds!=f)[0]; te_i=np.where(folds==f)[0]
    nf=fit(anom_net(),Xt[tr_i])
    with torch.no_grad(): rec_oos[te_i]=((nf(Xt[te_i])-Xt[te_i])**2).mean(1).numpy()

p4=PCA(4).fit(R); rec_pca=((R-p4.inverse_transform(p4.transform(R)))**2).mean(1)   # linear baseline
covid=int(np.argmin(np.abs((dates-pd.Timestamp("2020-03-20")).days)))
pct=lambda v,i:(v<v[i]).mean()

fig,ax=plt.subplots(1,2,figsize=(14,4.2),sharey=False)
for a,(v,t_) in zip(ax,[(rec_in,"Scored in-sample (the trap)"),(rec_oos,"Scored out-of-sample (6-fold)")]):
    a.plot(dates,v,color=BLUE,lw=.9); a.axhline(np.quantile(v,.95),color=RED,ls="--",lw=1,label="95th percentile")
    _t=np.argsort(v)[::-1][:5]; a.scatter(dates[_t],v[_t],color=RED,zorder=5,s=36)
    a.scatter(dates[covid],v[covid],color=GREEN,s=90,zorder=6,edgecolor="k",lw=.6,
              label=f"COVID crash: {pct(v,covid):.1%} percentile")
    a.set_ylabel("reconstruction error\n(mean over the 48 stocks)"); a.set_title(t_); a.legend(fontsize=8)
plt.tight_layout(); plt.show()

print(f"The network has {npar} parameters and {len(Xt)} weeks to fit, and the consequence is not subtle:")
print(f"  COVID crash week ({dates[covid].date()}), scored in-sample     : {pct(rec_in,covid):>6.1%} percentile")
print(f"  the same week, scored out-of-sample                  : {pct(rec_oos,covid):>6.1%} percentile")
print(f"  the same week, under a plain rank-4 PCA baseline     : {pct(rec_pca,covid):>6.1%} percentile")
print("Trained on it, the network reproduces that week almost exactly and reports it as one of the most ordinary in")
print("the sample. Held out, it is the most extraordinary. Nothing about the week changed -- only whether the model")
print("had already seen it.")
_ti=np.argsort(rec_in)[::-1][:5]; _to=np.argsort(rec_oos)[::-1][:5]
print(f"\n  top 5, in-sample     : {[str(dates[i].date()) for i in _ti]}")
print(f"  top 5, out-of-sample : {[str(dates[i].date()) for i in _to]}")
print(f"  weeks in common: {len(set(_ti)&set(_to))} of 5 -- the two rankings agree overall (Spearman "
      f"{spearmanr(rec_in,rec_oos).statistic:+.2f}) and disagree completely about the tail, which is the only part anyone uses.")
_mag=np.abs(R).mean(1); _ext=np.argsort(_mag)[::-1][:10]
print(f"\nThe signature of the problem: the ten largest-magnitude weeks average the {np.mean([pct(rec_in,i) for i in _ext]):.0%} percentile of")
print(f"in-sample error but the {np.mean([pct(rec_oos,i) for i in _ext]):.0%} percentile out-of-sample. Memorisation is strongest exactly where the")
print("data is most extreme, so the score is pushed down hardest on the points the detector exists to find.")
print(f"\nAnd it depends on the draw. Refitting the in-sample version under five different seeds puts that same COVID")
_pcts=[]
for _sd in range(5):
    _n=fit(anom_net(_sd),Xt)
    with torch.no_grad(): _r=((_n(Xt)-Xt)**2).mean(1).numpy()
    _pcts.append(pct(_r,covid))
print(f"week anywhere from the {min(_pcts):.0%} to the {max(_pcts):.0%} percentile ({', '.join(f'{v:.0%}' for v in _pcts)}) -- a conclusion")
print("that moves that much with the random seed was never carrying information about the market.")
print(f"\nDone properly, what does the detector find? Out-of-sample it agrees closely with the linear baseline")
print(f"(Spearman {spearmanr(rec_oos,rec_pca).statistic:+.2f} with rank-4 PCA) and correlates {spearmanr(rec_oos,_mag).statistic:+.2f} with sheer weekly magnitude, and the")
print(f"weeks it flags are the March 2020 crash. The nonlinearity is not buying a distinct signal on this cross-section:")
print("reconstruction error is a real anomaly score, but here it is largely measuring how big the week was, and a rank-4")
print("linear projection measures the same thing for none of the cost. That is a useful thing to have established")
print("rather than assumed -- and it is only visible once the score is computed on data the model was not fitted to.")

The network has 1748 parameters and 312 weeks to fit, and the consequence is not subtle:
  COVID crash week (2020-03-17), scored in-sample     :   5.8% percentile
  the same week, scored out-of-sample                  :  99.7% percentile
  the same week, under a plain rank-4 PCA baseline     :  99.4% percentile
Trained on it, the network reproduces that week almost exactly and reports it as one of the most ordinary in
the sample. Held out, it is the most extraordinary. Nothing about the week changed -- only whether the model
had already seen it.

  top 5, in-sample     : ['2022-10-25', '2022-04-26', '2022-02-01', '2024-09-10', '2024-11-05']
  top 5, out-of-sample : ['2020-03-17', '2020-03-24', '2022-05-17', '2020-03-10', '2020-03-03']
  weeks in common: 0 of 5 -- the two rankings agree overall (Spearman +0.83) and disagree completely about the tail, which is the only part anyone uses.

The signature of the problem: the ten largest-magnitude weeks average the 44% percentile of
in-sample error but the 91% percentile out-of-sample. Memorisation is strongest exactly where the
data is most extreme, so the score is pushed down hardest on the points the detector exists to find.

And it depends on the draw. Refitting the in-sample version under five different seeds puts that same COVID

week anywhere from the 0% to the 36% percentile (6%, 0%, 36%, 0%, 6%) -- a conclusion
that moves that much with the random seed was never carrying information about the market.

Done properly, what does the detector find? Out-of-sample it agrees closely with the linear baseline
(Spearman +0.91 with rank-4 PCA) and correlates +0.71 with sheer weekly magnitude, and the
weeks it flags are the March 2020 crash. The nonlinearity is not buying a distinct signal on this cross-section:
reconstruction error is a real anomaly score, but here it is largely measuring how big the week was, and a rank-4
linear projection measures the same thing for none of the cost. That is a useful thing to have established
rather than assumed -- and it is only visible once the score is computed on data the model was not fitted to.

5. Summary

An autoencoder learns compression by reconstructing its input through a bottleneck. We saw the full ladder:

• a linear autoencoder reproduces PCA's subspace exactly — identical reconstruction error to seven decimal places, principal angles of 0.000° — while its axes are an arbitrary basis of that subspace rather than PCA's ordered orthogonal components, because nothing in the loss asks for them;
• a nonlinear autoencoder beats PCA at the same bottleneck on Fashion-MNIST (14% lower test error, visibly sharper reconstructions) by bending to the data manifold;
• its 2-D latent space separates the image classes more cleanly than PCA, scored rather than eyeballed — 0.673 against 0.525 nearest-neighbour accuracy on labels neither map ever saw;
• its reconstruction error is an unsupervised anomaly score — but only when computed on data the network was not fitted to. Scored in-sample, a 1,748-parameter network on 312 weeks memorises its outliers and rates the COVID crash at the 6th percentile, the most normal end of the scale; held out, the same week sits at the 99.7th. The two rankings share none of their top five. Done properly the score tracks a rank-4 PCA baseline at Spearman 0.91, so on this cross-section the nonlinearity is not buying a distinct signal — a conclusion worth having tested rather than assumed.

The general lesson is the one this whole subsection turns on: unsupervised methods have no held-out label to be wrong against, so the discipline has to be supplied deliberately — a baseline the method must beat, and a score computed on data the model has not already fitted. Without both, a flexible model will report structure it has memorised.

Cross-links: autoencoders generalise the PCA/factor notebook (nonlinear vs linear compression); they reuse the neural-network machinery (encoder/decoder, backprop, Adam) from the deep-learning subsection; and the anomaly application sits alongside the volatility arc, which measures the magnitude of market moves directly. A probabilistic autoencoder (the variational autoencoder) would add a prior on the latent code and a full generative model — the Bayesian extension, and the bridge to generative modelling.

Next: t-SNE and UMAP — manifold-learning methods built specifically to visualise high-dimensional data in two dimensions, which we compare head-to-head with the PCA and autoencoder embeddings seen here.