---
title: Linear Factor Models
module: Generative Models
moduleNumber: 7
lessonNumber: 1
order: 701
summary: >
  The simplest generative models share one template: a latent variable drawn
  from a fixed prior, run through a linear decoder, plus noise. Probabilistic
  PCA, factor analysis, independent component analysis, and sparse coding are all
  this template with a different prior on the latents and a different noise
  model. We derive each marginal, see why ICA needs non-Gaussianity to
  identify its sources, and show how sparse coding learns Gabor-like dictionary atoms.
topics: [Generative Models]
sources:
  - book: Goodfellow
    ref: "Ch. 13 — Linear Factor Models"
  - book: Goodfellow
    ref: "§13.1 Probabilistic PCA and Factor Analysis; §13.2 ICA; §13.4 Sparse Coding"
---

A generative model proposes a process for how the data was made and then fits that
process to the data. The simplest such process is **linear**: a small set of
hidden causes $h$, drawn from a fixed prior, with each observation $x$ a
linear mix of those causes plus noise. Every model in this lesson (probabilistic
PCA, factor analysis, independent component analysis, sparse coding) instantiates
this process; they differ only in the prior $p(h)$ and the noise.[^gf-lfm]

## The shared template

A **linear factor model** posits a latent vector $h \in \mathbb{R}^k$ and an
observed vector $x \in \mathbb{R}^d$ related by a linear decoder:

$$
h \sim p(h),
\qquad
x = W h + b + \text{noise},
$$

where $W \in \mathbb{R}^{d \times k}$ is the **factor-loading matrix**, $b$ the
mean, and the noise is independent of $h$. The columns of $W$ are the **factors**:
$x = \sum_{j=1}^{k} h_j W_{:,j} + b + \text{noise}$ reads each observation as a
weighted sum of a fixed set of directions, the weights being the latents.

> **Definition (Linear factor model).** A generative model in which a latent
> $h \sim p(h)$ is mapped to the data by an affine decoder $x = Wh + b$ plus
> independent noise. The choice of $p(h)$ and of the noise distribution fixes the
> specific model.

The dependency structure is directed: the latent causes generate the observations,
never the reverse. Each $h_j$ fans out through column $W_{:,j}$ to influence every
coordinate of $x$.

$$
% caption: The directed graph of a linear factor model: latent factors $h_1,\dots,h_k$ fan out through $W$ to the observed $x_1,\dots,x_d$.
\begin{tikzpicture}[>=stealth, font=\footnotesize,
  lat/.style={circle, draw, thick, minimum size=9mm, inner sep=0pt},
  obs/.style={circle, draw, fill=black!8, minimum size=9mm, inner sep=0pt}]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  % latent layer
  \node[lat, draw=acc, text=acc] (h1) at (0,2.0)  {$h_1$};
  \node[lat, draw=acc, text=acc] (h2) at (0,0)    {$h_2$};
  \node[lat, draw=acc, text=acc] (h3) at (0,-2.0) {$h_3$};
  \node[font=\footnotesize, text=acc] at (0,3.0) {latent factors};
  % observed layer
  \node[obs] (x1) at (5,2.7)  {$x_1$};
  \node[obs] (x2) at (5,1.35) {$x_2$};
  \node[obs] (x3) at (5,0)    {$x_3$};
  \node[obs] (x4) at (5,-1.35){$x_4$};
  \node[obs] (x5) at (5,-2.7) {$x_5$};
  \node[font=\footnotesize] at (5,3.6) {observed};
  % fan-out edges (W)
  \foreach \h in {h1,h2,h3}
    \foreach \x in {x1,x2,x3,x4,x5}
      \draw[->, black] (\h) -- (\x);
  \node[font=\footnotesize, text=acc] at (2.5,3.15) {$W$};
\end{tikzpicture}
$$

Fixing the prior to a standard Gaussian and the decoder to be linear makes the
whole joint distribution tractable, because affine maps of Gaussians stay
Gaussian. That is the regime of the first two models.

## Probabilistic PCA and factor analysis

Take the latent prior to be a unit Gaussian and the noise to be Gaussian as well.
**Factor analysis** allows the noise to have a different variance on each
coordinate; **probabilistic PCA** forces one shared variance. Both write

$$
h \sim \mathcal{N}(0, I_k),
\qquad
x \mid h \sim \mathcal{N}(W h + b,\; \Psi),
$$

with the noise covariance $\Psi$ diagonal. The two models are the two choices of
$\Psi$:

$$
\underbrace{\Psi = \diag(\psi_1, \dots, \psi_d)}_{\text{factor analysis}}
\qquad\text{vs.}\qquad
\underbrace{\Psi = \sigma^2 I_d}_{\text{probabilistic PCA}}.
$$

> **Definition (Factor analysis).** The linear-Gaussian model with $h \sim
> \mathcal{N}(0,I)$ and $x \mid h \sim \mathcal{N}(Wh+b, \Psi)$ where $\Psi$ is a
> general diagonal matrix: each observed coordinate gets its own residual
> variance $\psi_i$.

> **Definition (Probabilistic PCA).** Factor analysis with isotropic noise,
> $\Psi = \sigma^2 I$, a single residual variance shared across all coordinates.

### The marginal over the data

Because $x$ is an affine function of the jointly Gaussian pair $(h, \text{noise})$,
the marginal $p(x)$ is Gaussian, and we only need its mean and covariance. Using
$\mathbb{E}[h]=0$ and $\Cov(h)=I$, independence of $h$ and the noise:

$$
\mathbb{E}[x] = W\,\mathbb{E}[h] + b = b,
$$

$$
\Cov(x)
= \mathbb{E}\!\brackets{(Wh + \text{noise})(Wh + \text{noise})^{T}}
= W\,\underbrace{\mathbb{E}[h h^{T}]}_{I}\,W^{T} + \underbrace{\Cov(\text{noise})}_{\Psi}
= W W^{T} + \Psi.
$$

So the model says the data is a Gaussian whose covariance is **low-rank plus
diagonal**:

$$
x \sim \mathcal{N}\parens{b,\; W W^{T} + \Psi}.
$$

The $k$ columns of $W$ supply a rank-$k$ chunk of the covariance, the shared,
correlated structure, while $\Psi$ mops up the per-coordinate residual. This is
the whole content of the models: explain the **correlations** in $x$ with few
factors, and dump everything left over into independent noise.[^gf-fa]

> **Theorem (PPCA recovers PCA in the zero-noise limit).** As $\sigma^2 \to 0$ in
> probabilistic PCA, the posterior mean $\mathbb{E}[h \mid x]$ converges to the
> orthogonal projection of $x - b$ onto the subspace spanned by the top $k$
> principal directions of the data covariance, i.e. classical PCA.

> **Proof.** The Gaussian posterior has mean $\mathbb{E}[h\mid x] = (W^{T}W +
> \sigma^2 I)^{-1} W^{T}(x-b)$. The maximum-likelihood $W$ has columns aligned with
> the top-$k$ eigenvectors $U_k$ of the sample covariance, scaled by
> $(\lambda_i - \sigma^2)^{1/2}$. As $\sigma^2 \to 0$ the term $\sigma^2 I$
> vanishes and the scaling cancels in $(W^{T}W)^{-1}W^{T}$, leaving the projection
> $U_k^{T}(x-b)$ onto the principal subspace, the classical PCA encoding. $\qed$

Classical [PCA](/deep-learning/mathematical-background/linear-algebra-for-deep-learning)
is therefore the $\sigma^2 \to 0$ corner of this family: a deterministic
projection onto the leading eigenvectors, with no noise model at all. The
probabilistic version adds back a likelihood, which is what lets it generate new
samples and handle missing data.[^gf-ppca]

The geometric difference between PPCA and FA is the **shape of the noise**.
Isotropic noise (PPCA) is a round blur identical on every axis. Diagonal noise
(FA) is an axis-aligned blur that can be wider on some coordinates than others.

$$
% caption: PCA versus factor analysis on a 2D cloud: PPCA adds isotropic noise (a round blur), factor analysis per-axis noise (an axis-aligned ellipse).
\begin{tikzpicture}[>=stealth, font=\footnotesize, scale=1.0]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  \definecolor{red}{HTML}{C0392B}
  % ---- left: PPCA, isotropic noise ----
  \begin{scope}
    \draw[->, black] (-2.2,0) -- (2.4,0) node[right, black]{$x_1$};
    \draw[->, black] (0,-2.0) -- (0,2.4) node[above, black]{$x_2$};
    % data cloud points along a diagonal
    \foreach \px/\py in {-1.5/-1.2,-1.0/-0.7,-0.6/-0.6,-0.2/0.1,0.2/0.2,0.6/0.5,1.0/0.9,1.4/1.1,0.4/0.7,-0.8/-0.3}
      \fill[black] (\px,\py) circle (1.6pt);
    % principal direction
    \draw[->, acc, very thick] (0,0) -- (1.4,1.1) node[anchor=south east, text=acc]{factor};
    % isotropic noise: round circle
    \draw[green, thick] (0.6,0.5) circle (0.42);
    \node[green, anchor=east] at (-0.15,1.5) {round};
    \node[anchor=north] at (0,-2.2) {PPCA: isotropic noise};
  \end{scope}
  % ---- right: FA, per-axis noise ----
  \begin{scope}[xshift=6.4cm]
    \draw[->, black] (-2.2,0) -- (2.4,0) node[right, black]{$x_1$};
    \draw[->, black] (0,-2.0) -- (0,2.4) node[above, black]{$x_2$};
    \foreach \px/\py in {-1.5/-1.2,-1.0/-0.7,-0.6/-0.6,-0.2/0.1,0.2/0.2,0.6/0.5,1.0/0.9,1.4/1.1,0.4/0.7,-0.8/-0.3}
      \fill[black] (\px,\py) circle (1.6pt);
    \draw[->, acc, very thick] (0,0) -- (1.4,1.1) node[anchor=south east, text=acc]{factor};
    % per-axis noise: axis-aligned ellipse (wider on x1)
    \draw[red, thick] (0.6,0.5) ellipse (0.62 and 0.30);
    \node[red, anchor=west] at (1.45,0.5) {per-axis};
    \node[anchor=north] at (0,-2.2) {FA: per-axis noise};
  \end{scope}
\end{tikzpicture}
$$

Both models are fit by maximum likelihood — closed-form for PCA, by
expectation-maximization for general FA and PPCA, alternating between inferring the
posterior over $h$ and re-estimating $W, \Psi$.

### A worked covariance decomposition

For example, take $d = 3$ observed
coordinates and $k = 1$ latent factor, with a factor-loading vector and diagonal
noise

$$
W = \begin{bmatrix} 2 \\ 1 \\ -1 \end{bmatrix},
\qquad
\Psi = \diag(0.5,\; 1.0,\; 0.25).
$$

The model predicts $\Cov(x) = W W^{T} + \Psi$. The rank-one outer
product is

$$
W W^{T}
= \begin{bmatrix} 2 \\ 1 \\ -1 \end{bmatrix}
  \begin{bmatrix} 2 & 1 & -1 \end{bmatrix}
= \begin{bmatrix} 4 & 2 & -2 \\ 2 & 1 & -1 \\ -2 & -1 & 1 \end{bmatrix},
$$

and adding $\Psi$ to the diagonal gives the full model covariance

$$
\Sigma
= W W^{T} + \Psi
= \begin{bmatrix} 4.5 & 2 & -2 \\ 2 & 2 & -1 \\ -2 & -1 & 1.25 \end{bmatrix}.
$$

$$
% caption: Low-rank-plus-diagonal covariance. The full model covariance $\Sigma$ splits into a rank-$k$ block $WW^T$ carrying every correlation and a diagonal $\Psi$ of private per-coordinate noise.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{red}{HTML}{C0392B}
  % Sigma block (dense)
  \draw[acc, thick, fill=acc!10] (0,0) rectangle (1.6,1.6);
  \node at (0.8,0.8) {Cov$(x)$};
  \node[font=\footnotesize, anchor=north] at (0.8,-0.15) {dense};
  \node at (2.35,0.8) {$=$};
  % WW^T block (low rank)
  \draw[acc, thick] (3.1,0) rectangle (4.7,1.6);
  % shading a thin band to suggest rank-1
  \draw[acc, thick, fill=acc!18] (3.1,0.55) rectangle (4.7,1.05);
  \node at (3.9,0.8) {$WW^{T}$};
  \node[font=\footnotesize, anchor=north] at (3.9,-0.15) {rank $k$: correlations};
  \node at (5.45,0.8) {$+$};
  % Psi block (diagonal)
  \draw[red, thick] (6.2,0) rectangle (7.8,1.6);
  \foreach \i in {0,1,2} {
    \draw[red, thick, fill=red!15] (6.2+\i*0.533,1.6-\i*0.533-0.533) rectangle (6.2+\i*0.533+0.533,1.6-\i*0.533);
  }
  \node[font=\footnotesize, anchor=north] at (7.0,-0.15) {diagonal: noise};
\end{tikzpicture}
$$

Read the structure directly off the matrix. The **off-diagonal** entries come
entirely from $WW^{T}$: coordinates $1$ and $2$ co-vary at $+2$ because both load
positively on the factor, while coordinate $3$ anti-correlates ($-2, -1$) because
it loads negatively. The single hidden cause explains _every_ correlation among the
three observations. The **diagonal** splits into a shared part and a private part:
variance of $x_1$ is $4$ (from the factor) plus $0.5$ (its own noise). A fit with
$k$ factors would explain the covariance with $3k + 3$ numbers here ($W$ plus
$\Psi$) instead of the $6$ free entries of a full $3\times3$ covariance — the
compression that makes the model useful grows as $d$ climbs and $k$ stays small.

The posterior over the latent given an observation follows from the same Gaussian
algebra. For a point $x = (3, 1, -1)^{T}$ with $b = 0$, the posterior mean is
$\mathbb{E}[h \mid x] = (W^{T}\Psi^{-1}W + 1)^{-1} W^{T}\Psi^{-1}(x - b)$. Here
$W^{T}\Psi^{-1}W = 4/0.5 + 1/1 + 1/0.25 = 8 + 1 + 4 = 13$, and
$W^{T}\Psi^{-1}x = 2\cdot3/0.5 + 1\cdot1/1 + (-1)(-1)/0.25 = 12 + 1 + 4 = 17$, so

$$
\mathbb{E}[h \mid x] = \frac{17}{13 + 1} = \frac{17}{14} \approx 1.21.
$$

The point is best explained by a factor firing at about $1.2$; multiplying back
through $W$ gives the reconstruction $\hat x = W\cdot1.21 \approx (2.43, 1.21,
-1.21)$, and the gap $x - \hat x$ is what the noise model absorbs.

## Independent component analysis

Gaussian latents are convenient but they throw away information: a Gaussian is
rotationally symmetric, so $W$ is only identifiable up to an orthogonal rotation,
and you cannot tell which factor is which. **Independent component analysis** (ICA)
abandons the Gaussian prior precisely to recover the individual sources. It
assumes the latents are statistically **independent** and **non-Gaussian**, and
the decoder is (in the basic form) noise-free and square:

$$
x = W h,
\qquad
p(h) = \prod_{j=1}^{k} p_j(h_j),
\qquad
p_j \text{ non-Gaussian}.
$$

> **Definition (Independent component analysis).** A linear factor model with
> statistically independent, non-Gaussian latents and (typically) no noise, fit by
> finding an unmixing matrix $W^{-1}$ that makes the recovered components
> $h = W^{-1} x$ maximally independent and non-Gaussian.

The canonical setting is the **cocktail-party problem**: $k$ microphones record $k$
speakers talking at once, each microphone picking up a different linear mixture.
ICA recovers the original voices (**blind source separation**) knowing only the
mixtures, not the mixing matrix.[^gf-ica]

$$
% caption: ICA on the cocktail-party problem: two mixed signals are separated back into the original waveforms by the estimated unmixing matrix $W^{-1}$.
\begin{tikzpicture}[>=stealth, font=\footnotesize, scale=1.0]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  \definecolor{red}{HTML}{C0392B}
  % ---- mixed signals (left) ----
  \node[font=\footnotesize] at (1.6,2.3) {mixed observations $x$};
  % mixed signal 1
  \draw[acc, thick] plot[domain=0:3.2, samples=80] (\x, {1.3 + 0.32*sin(\x*180) + 0.26*sin(\x*430+40)});
  % mixed signal 2
  \draw[red, thick] plot[domain=0:3.2, samples=80] (\x, {0.0 + 0.30*sin(\x*230+60) - 0.24*sin(\x*150)});
  % arrow: unmixing
  \draw[->, very thick, black] (3.7,0.7) -- (5.3,0.7);
  \node[font=\footnotesize, anchor=south] at (4.5,0.80) {unmix};
  \node[font=\footnotesize, anchor=north, text=acc] at (4.5,0.60) {ICA};
  % ---- separated sources (right) ----
  \begin{scope}[xshift=5.8cm]
    \node[font=\footnotesize] at (1.6,2.3) {recovered sources $h$};
    % clean source 1: pure low tone
    \draw[acc, thick] plot[domain=0:3.2, samples=80] (\x, {1.3 + 0.42*sin(\x*180)});
    % clean source 2: pure high tone
    \draw[green, thick] plot[domain=0:3.2, samples=80] (\x, {0.0 + 0.38*sin(\x*430+40)});
  \end{scope}
\end{tikzpicture}
$$

### Why non-Gaussianity is required

The identifiability argument is sharp. Suppose the sources were Gaussian:
$h \sim \mathcal{N}(0, I)$. For any orthogonal $R$ ($R R^{T} = I$), the rotated
latent $h' = R h$ is **also** $\mathcal{N}(0, I)$ — the standard Gaussian is
rotation-invariant. Then the mixing $W' = W R^{T}$ produces an identical data
distribution:

$$
W' h' = W R^{T} R\, h = W h = x,
$$

so $W$ and $WR^{T}$ are indistinguishable from the data alone, for every rotation
$R$. The sources are recoverable only up to an arbitrary rotation, exactly the
non-identifiability that FA inherits. **Non-Gaussian** sources break this
symmetry: a non-Gaussian distribution is _not_ rotation-invariant, so only the true
unmixing makes the recovered components independent.

> **Theorem (ICA identifiability).** If at most one source $h_j$ is Gaussian and
> the sources are independent, the mixing matrix $W$ is identifiable up to
> permutation and per-component scaling of its columns.

> **Proof sketch.** Independence plus non-Gaussianity forces a unique rotation. The
> Darmois–Skitovich theorem states that if two independent variables have a linear
> combination that is independent of another, the components involved must be
> Gaussian. Contrapositively, with the non-Gaussian sources independent, no
> nontrivial rotation can keep the recovered components independent, leaving only
> the trivial freedoms of relabeling the sources (permutation) and rescaling each
> (you cannot separate $h_j$ from $2 h_j$ with $W_{:,j}/2$). $\qed$

$$
% caption: Why non-Gaussianity is needed. A Gaussian source cloud (left) is rotationally symmetric, so no rotation is distinguished; a non-Gaussian cloud (right) has corners that pin the axes to the true sources.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  % --- left: Gaussian, round cloud ---
  \begin{scope}
    \draw[->, black] (-2.1,0) -- (2.1,0) node[right, font=\scriptsize]{$h_1$};
    \draw[->, black] (0,-2.1) -- (0,2.1) node[above, font=\scriptsize]{$h_2$};
    \foreach \px/\py in {0.1/0.2,-0.4/0.5,0.6/-0.3,-0.7/-0.6,0.9/0.4,-0.2/0.9,0.4/-0.8,-0.9/0.3,0.2/-0.4,-0.5/-0.2,0.7/0.7,-0.6/0.8,0.3/0.5,-0.3/-0.7,0.8/-0.5}
      \fill[black] (\px,\py) circle (1.5pt);
    \draw[green, thick] (0,0) circle (1.25);
    % ambiguous rotated axes
    \draw[green, dashed, thick] (-1.3,-1.3) -- (1.3,1.3);
    \draw[green, dashed, thick] (-1.3,1.3) -- (1.3,-1.3);
    \node[font=\scriptsize, anchor=north] at (0,-2.35) {Gaussian: any rotation f\/its};
  \end{scope}
  % --- right: non-Gaussian, square/heavy-tailed cloud ---
  \begin{scope}[xshift=6.2cm]
    \draw[->, black] (-2.1,0) -- (2.1,0) node[right, font=\scriptsize]{$h_1$};
    \draw[->, black] (0,-2.1) -- (0,2.1) node[above, font=\scriptsize]{$h_2$};
    \foreach \px/\py in {1.3/1.2,-1.3/1.1,1.2/-1.3,-1.2/-1.2,1.4/0.2,-1.4/0.3,0.2/1.3,0.1/-1.4,1.1/1.4,-1.4/-1.0,1.3/-0.9,-1.0/1.4}
      \fill[black] (\px,\py) circle (1.5pt);
    \draw[acc, thick] (-1.5,-1.5) rectangle (1.5,1.5);
    \draw[acc, very thick, ->] (0,0) -- (1.5,0);
    \draw[acc, very thick, ->] (0,0) -- (0,1.5);
    \node[font=\scriptsize, anchor=north] at (0,-2.35) {non-Gaussian: axes pinned};
  \end{scope}
\end{tikzpicture}
$$

Practically, ICA maximizes a non-Gaussianity contrast (kurtosis or negentropy)
of the recovered components $h = W^{-1}x$, since by the central limit theorem a
mixture of independent signals is _more_ Gaussian than its parts, so the most
non-Gaussian projections are the un-mixed sources.[^gf-ica-id]

For example, suppose two independent sources are mixed by

$$
W = \begin{bmatrix} 1 & 1 \\ 0.5 & 2 \end{bmatrix},
\qquad
x = W h.
$$

Given the mixtures $x$, ICA searches for the unmixing matrix $W^{-1}$; inverting the
$2\times2$ (determinant $1\cdot2 - 1\cdot0.5 = 1.5$) gives

$$
W^{-1} = \frac{1}{1.5}\begin{bmatrix} 2 & -1 \\ -0.5 & 1 \end{bmatrix}
= \begin{bmatrix} 1.33 & -0.67 \\ -0.33 & 0.67 \end{bmatrix}.
$$

For a specific sample $x = (3, 3)^{T}$, applying the estimated unmixing recovers
$h = W^{-1}x = (1.33\cdot3 - 0.67\cdot3,\; -0.33\cdot3 + 0.67\cdot3) = (2, 1)$ —
the original source amplitudes, up to the scale and sign freedoms in the theorem.
In practice $W^{-1}$ is not handed to the algorithm; it is _found_ by
rotating the whitened data until the two output coordinates are maximally
non-Gaussian, at which point their scatter looks like the pinned square on the
right above rather than the ambiguous round blob on the left.

## Sparse coding

The previous models keep the latent dimension small. **Sparse coding** goes the
other way: it uses an **overcomplete** dictionary ($k > d$, more factors than
data dimensions) but insists that each observation uses only a _few_ of them. The
prior on $h$ is sharply peaked at zero (a Laplace, equivalently an $L_1$ penalty),
so most components are exactly zero.

$$
x = W h + \text{noise},
\qquad
p(h) \propto \prod_{j=1}^{k} e^{-\lambda |h_j|}
\quad(\text{Laplace, } L_1),
\qquad k > d.
$$

> **Definition (Sparse coding).** A linear factor model with an overcomplete
> dictionary $W$ and a sparsity-inducing prior on $h$ (Laplace / $L_1$): each
> input is reconstructed as a sparse combination of a few dictionary atoms.

### Inference is an optimization

There is no closed-form posterior. For a fixed dictionary $W$, the code $h$ for an
input $x$ is the MAP estimate: minimize reconstruction error plus the $L_1$
penalty that comes from the Laplace prior:

$$
h^\star = \arg\min_{h}\;
\underbrace{\tfrac{1}{2}\,\lVert x - W h \rVert_2^2}_{\text{reconstruction}}
\;+\;
\underbrace{\lambda \,\lVert h \rVert_1}_{\text{sparsity}}.
$$

This is the **LASSO**: convex in $h$, with a non-smooth $L_1$ term that drives most
coordinates to exactly zero. Learning alternates between this sparse inference and
a dictionary update on $W$.

```algorithm
caption: $\textsc{SparseCoding}(X, k, \lambda)$ — dictionary learning by alternating minimization
initialize dictionary $W \in \mathbb{R}^{d \times k}$ randomly, columns unit-norm
repeat
  for each input $x$ in the batch do // sparse inference (LASSO)
    $h \gets \arg\min_{h}\ \tfrac{1}{2}\lVert x - W h\rVert_2^2 + \lambda\lVert h\rVert_1$
  $W \gets \arg\min_{W}\ \sum_x \tfrac{1}{2}\lVert x - W h_x\rVert_2^2$ // dictionary update
  renormalize each column of $W$ to unit length // prevent $W$ blow-up
until converged
return $W$
```

Geometrically, the input is reconstructed as a sparse weighted sum: out of a large
dictionary, only a handful of atoms switch on, each scaled by its (nonzero) code.

$$
% caption: Sparse coding reconstructs an input as a weighted sum of dictionary atoms, with only a few of the $k$ codes $h_j$ nonzero.
\begin{tikzpicture}[>=stealth, font=\footnotesize,
  atom/.style={draw, minimum width=9mm, minimum height=9mm, inner sep=0pt},
  on/.style={draw=acc, thick, fill=acc!12},
  off/.style={draw=black, fill=black!8}]
  \definecolor{acc}{HTML}{2348F2}
  % dictionary atoms row
  \node[atom, on]  (a1) at (0,0)   {};
  \node[atom, off] (a2) at (1.1,0) {};
  \node[atom, on]  (a3) at (2.2,0) {};
  \node[atom, off] (a4) at (3.3,0) {};
  \node[atom, off] (a5) at (4.4,0) {};
  \node[atom, on]  (a6) at (5.5,0) {};
  \node[font=\footnotesize] at (2.75,1.0) {dictionary $W$ ($k$ atoms)};
  % coefficients under active atoms
  \node[text=acc, font=\footnotesize] at (0,-0.85)   {$h_1$};
  \node[text=black, font=\footnotesize] at (1.1,-0.85) {$0$};
  \node[text=acc, font=\footnotesize] at (2.2,-0.85)  {$h_3$};
  \node[text=black, font=\footnotesize] at (3.3,-0.85) {$0$};
  \node[text=black, font=\footnotesize] at (4.4,-0.85) {$0$};
  \node[text=acc, font=\footnotesize] at (5.5,-0.85)  {$h_6$};
  % sum to input
  \node[font=\footnotesize] at (2.75,-1.7) {only the \texttt{active} set (atoms 1, 3, 6) contributes};
  \node[atom, draw=acc, thick] (xr) at (8.0,0) {$x$};
  \node[font=\footnotesize] at (8.0,1.0) {input};
  \draw[->, acc, thick] (a6) -- (xr) node[midway, above, font=\footnotesize]{sum};
\end{tikzpicture}
$$

### Learned dictionaries are Gabor-like

Trained on natural image patches, the learned atoms reliably become **localized,
oriented, band-pass edge detectors**, the discrete analogue of **Gabor filters**
and a close match to the receptive fields of simple cells in the primary visual
cortex. Sparsity is doing the work: forcing each patch to be explained by few atoms
pushes the dictionary toward a basis of elementary edges that compose efficiently.[^gf-sparse-code]

$$
% caption: A learned sparse-coding dictionary on natural-image patches: each tile is one atom, a localized oriented edge or grating, the discrete analogue of Gabor filters.
\begin{tikzpicture}[>=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{2348F2}
  % a 4x4 grid of patches, each an oriented edge / grating icon
  \foreach \col in {0,1,2,3}
    \foreach \row in {0,1,2,3} {
      \begin{scope}[shift={(\col*1.25, -\row*1.25)}]
        \draw[black] (0,0) rectangle (1.05,1.05);
      \end{scope}
    }
  % atom 0,0: vertical edge
  \draw[very thick] (0.52,0.12) -- (0.52,0.93);
  % atom 1,0: horizontal edge
  \draw[very thick] (1.37,0.52) -- (2.18,0.52);
  % atom 2,0: diagonal /
  \draw[very thick] (2.62,0.15) -- (3.43,0.9);
  % atom 3,0: diagonal \
  \draw[very thick] (3.87,0.9) -- (4.68,0.15);
  % row 1 (y offset -1.25): gratings
  \draw[very thick] (0.30,-1.13) -- (0.30,-0.32);
  \draw[very thick] (0.52,-1.13) -- (0.52,-0.32);
  \draw[very thick] (0.74,-1.13) -- (0.74,-0.32);
  \draw[very thick] (1.30,-0.74) -- (2.11,-0.74);
  \draw[very thick] (1.30,-0.52) -- (2.11,-0.52);
  \draw[very thick] (2.66,-1.1) -- (3.40,-0.36);
  \draw[very thick] (2.86,-1.1) -- (3.50,-0.46);
  \draw[very thick] (3.95,-0.40) -- (4.62,-1.05);
  \draw[very thick] (3.80,-0.52) -- (4.50,-1.13);
  % row 2 (y -2.5): mixed orientations
  \draw[very thick] (0.18,-1.9) -- (0.92,-2.25);
  \draw[very thick] (1.55,-2.4) -- (1.95,-1.7);
  \draw[very thick] (2.62,-2.35) -- (3.43,-1.65);
  \draw[very thick] (3.92,-1.65) -- (4.58,-2.35);
  % row 3 (y -3.75)
  \draw[very thick] (0.52,-3.6) -- (0.52,-2.85);
  \draw[very thick] (1.30,-3.25) -- (2.11,-3.05);
  \draw[very thick] (2.55,-3.05) -- (3.50,-3.05);
  \draw[very thick] (3.85,-2.85) -- (4.65,-3.6);
\end{tikzpicture}
$$

## The family, side by side

All four models are one template, $x = Wh + b + \text{noise}$, distinguished by
the prior on $h$, the noise, and what structure they recover.

| model | prior $p(h)$ | noise | dimension | recovers |
| --- | --- | --- | --- | --- |
| probabilistic PCA | $\mathcal{N}(0, I)$ | isotropic $\sigma^2 I$ | $k < d$ | principal subspace |
| factor analysis | $\mathcal{N}(0, I)$ | diagonal $\diag(\psi_i)$ | $k < d$ | shared correlations |
| ICA | independent, non-Gaussian | none (or small) | $k = d$ | independent sources |
| sparse coding | Laplace / $L_1$ (sparse) | Gaussian | $k > d$ (overcomplete) | sparse edge dictionary |

Gaussian priors (PPCA, FA) give tractability at the
cost of identifiability; ICA drops the Gaussian to make individual sources
identifiable; sparse coding drops the small latent dimension for an overcomplete,
interpretable basis. Each is a different answer to the same question: what simple
hidden causes, linearly combined, produced this data?

## Linear factors in the deep era

The linear decoder looks primitive next to a neural network, yet each model in
this lesson has a direct descendant in current practice, and the reasons they were
introduced still hold.

**Sparse coding survives as dictionary learning and $L_1$ regularization.**
Olshausen and Field's original result — that a sparsity prior on natural image
patches yields Gabor-like edge filters — was the first computational account of why
V1 simple cells look the way they do, and it remains the standard reference for
sparse dictionary methods.[^olshausen] The same $L_1$ penalty reappears everywhere a
model should use few active units: LASSO regression, sparse attention, and the
sparse **autoencoders** now used to pull interpretable features out of large
language models are the same objective with a learned nonlinear encoder in place of
the LASSO solve.

**ICA is still the default for blind source separation.** The FastICA algorithm,
which maximizes negentropy by a fixed-point iteration, made ICA fast enough for
routine use and is the tool of choice for removing eye-blink and heartbeat
artifacts from EEG and MEG recordings, where the independent sources are literally
distinct physical generators.[^fastica] The non-Gaussianity argument from this
lesson explains why it works: neural and muscular sources are heavy-tailed, so
their independent directions are identifiable.

**The template itself became the encoder–decoder.** Replacing the linear map
$Wh + b$ with a deep network, and the closed-form Gaussian posterior with an
_amortized_ neural approximation, turns factor analysis into the **variational
autoencoder** of Kingma and Welling — same latent-variable-and-decoder graph, same
Gaussian prior, now nonlinear and trained by stochastic gradients.[^vae-kw] Read
this way, the rest of the chapter is one long answer to a question the linear
models raised and could not fully answer: how to keep the latent-variable structure
when the decoder is an arbitrary neural network.

[^gf-lfm]: **Goodfellow**, _Deep Learning_, Ch. 13 — linear factor models as the common template $x = Wh + b + \text{noise}$, with the prior $p(h)$ and noise model distinguishing each member.
[^gf-fa]: **Goodfellow**, _Deep Learning_, §13.1 — factor analysis: a Gaussian latent prior and diagonal noise give the data a low-rank-plus-diagonal covariance $WW^{T}+\Psi$.
[^gf-ppca]: **Goodfellow**, _Deep Learning_, §13.1 — probabilistic PCA as factor analysis with isotropic noise $\sigma^2 I$, recovering classical PCA in the $\sigma^2\to0$ limit.
[^gf-ica]: **Goodfellow**, _Deep Learning_, §13.2 — independent component analysis and blind source separation: independent, non-Gaussian latents and a (near-)noiseless mixing.
[^gf-ica-id]: **Goodfellow**, _Deep Learning_, §13.2 — why a Gaussian prior leaves $W$ identifiable only up to rotation, and how non-Gaussianity (kurtosis/negentropy contrasts) pins the sources down.
[^gf-sparse-code]: **Goodfellow**, _Deep Learning_, §13.4 — sparse coding with an overcomplete dictionary and an $L^1$/Laplace prior, and the Gabor-like atoms it learns on natural image patches.
[^olshausen]: **Olshausen & Field**, "Emergence of simple-cell receptive field properties by learning a sparse code for natural images," _Nature_ 381 (1996) — sparse coding on natural image patches yields localized, oriented, band-pass filters matching V1 simple cells.
[^fastica]: **Hyvärinen & Oja**, "Independent Component Analysis: Algorithms and Applications," _Neural Networks_ 13 (2000) — the FastICA fixed-point algorithm maximizing negentropy, and its use in EEG/MEG artifact removal.
[^vae-kw]: **Kingma & Welling**, "Auto-Encoding Variational Bayes," ICLR 2014 — the variational autoencoder: factor analysis's latent-variable template with a deep decoder and an amortized approximate posterior.