Predicting one random variable from another
\[ % Number sets \newcommand{\Nat}{\mathbb{N}} \newcommand{\Integer}{\mathbb{Z}} \newcommand{\Real}{\mathbb{R}} \newcommand{\Complex}{\mathbb{C}} % Operators \newcommand{\SNR}{\operatorname{SNR}} \newcommand{\tr}{\operatorname{tr}} \newcommand{\re}{\operatorname{Re}} \newcommand{\im}{\operatorname{Im}} \newcommand{\sign}{\operatorname{sign}} \newcommand{\rect}{\operatorname{rect}} \newcommand{\sinc}{\operatorname{si}} \newcommand{\erf}{\operatorname{erf}} \newcommand{\erfc}{\operatorname{erfc}} \newcommand{\kurt}{\operatorname{kurt}} \newcommand{\diag}{\operatorname{diag}} \newcommand{\MSE}{\operatorname{MSE}} \newcommand{\argmin}{\operatorname*{argmin}} \newcommand{\argmax}{\operatorname*{argmax}} \newcommand{\bias}{\operatorname{bias}} \newcommand{\Bin}{\operatorname{Bin}} % Delimiters \newcommand{\paren}[1]{\left(#1\right)} \newcommand{\bracket}[1]{\left[#1\right]} \newcommand{\set}[1]{\left\{#1\right\}} \newcommand{\abs}[1]{\left\lvert #1 \right\rvert} \newcommand{\norm}[1]{\left\lVert #1 \right\rVert} \newcommand{\ip}[1]{\left\langle #1 \right\rangle} % Notation \newcommand{\dd}{\,\mathrm{d}} \newcommand{\tran}{^{\top}} \newcommand{\conv}{\ast} \newcommand{\conj}{^{*}} \newcommand{\herm}{^{\mathsf{H}}} \newcommand{\opt}{_{\star}} % Expectation and variance % \newcommand{\E}[1]{\mathbb{E}\!\left[#1\right]} \newcommand{\E}[2][]{\mathbb{E}_{#1}\!\left[#2\right]} \newcommand{\Var}[1]{\mathbb{V}\!\left[#1\right]} % Distributions \newcommand{\Gauss}{\mathcal{N}} \newcommand{\Uniform}{\mathcal{U}} \newcommand{\FT}{\mathcal{F}} % Probability \newcommand{\Prob}[1]{\mathbb{P}\!\left(#1\right)} % Misc \newcommand{\bsy}[1]{\boldsymbol{#1}} \newcommand{\calR}{\mathcal{R}} \newcommand{\calC}{\mathcal{C}} \]
\[ \text{Observations }\mathbf{X}\ \overset{\textstyle\text{predict}}{\longrightarrow}\ \text{Estimand } Y \]
Find the best mapping \(g\)
\[ \widehat{Y} = g(\mathbf{X}) \]
We want \(\widehat{Y}\) close to \(Y\).
Measure it with a cost \(\mathcal{C}(Y,\widehat{Y})\) and average:
\[ g^\star = \argmin_{g}\ \E[\mathbf{X},Y]{\mathcal{C}\big(Y,\, g(\mathbf{X})\big)} \]
where the expectation is over the joint distribution of \((\mathbf{X}, Y)\).
Split the joint expectation — outer over \(\mathbf{X}\), inner over \(Y\mid\mathbf{X}\):
\[ \E[\mathbf{X},Y]{\mathcal{C}\big(Y,\, g(\mathbf{X})\big)} = \E[\mathbf{X}]{\ \underbrace{\E[Y\mid \mathbf{X}]{\mathcal{C}\big(Y, g(\mathbf{X})\big)}}_{\text{conditional risk}}\ } \]
Minimize inside, separately for each \(\mathbf{x}\).
For fixed \(\mathbf{x}\), the value \(g(\mathbf{x}) = \widehat{y}\) is just a constant.
\[ g^\star(\mathbf{x}) = \argmin_{\widehat{y}} \E[Y\mid \mathbf{X}=\mathbf{x}]{\mathcal{C}\big(Y, \widehat{y}\big)} \]
Squared error \(\mathcal{C} = (Y-\widehat{y})^2 \;\Rightarrow\; \widehat{y} = \E{Y\mid\mathbf{x}}\) — conditional mean (MMSE)
Hit-or-miss \(\mathcal{C} = \mathbb{1}[\,Y\neq\widehat{y}\,] \;\Rightarrow\; \widehat{y} = \argmax_{y} f_{Y\mid X}(y\mid\mathbf{x})\) — conditional mode (MAP)
These are the best unrestricted predictors. But they require knowledge of the full conditional distribution \(f_{Y\mid X}(y\mid\mathbf{x})\).
| mapping | |
|---|---|
| General — any function | \(g(\mathbf{X})\) |
| Parametric — e.g. a neural network | \(g_{\bsy{\phi}}(\mathbf{X})\) |
| Linear | \(g_{\bsy{\phi}}(\mathbf{X}) = \bsy{\phi}\tran\mathbf{X}\) |
| Affine | \(g_{\bsy{\phi}}(\mathbf{X}) = \bsy{\phi}\tran\mathbf{X} + \phi_0\) |
Optimizing over functions becomes optimizing over parameters:
\[ \argmin_{g}\ \E{\mathcal{C}(Y, g(\mathbf{X}))} \quad\longrightarrow\quad \argmin_{\bsy{\phi}}\ \E{\mathcal{C}(Y, g_{\bsy{\phi}}(\mathbf{X}))} \]
Allow a feature map \(\bsy{\psi}(\mathbf{X})\in\Real^K\):
\[ \widehat{Y} = \bsy{\phi}\tran\bsy{\psi}(\mathbf{X}) \]
Substitute \(\mathbf{Z} = \bsy{\psi}(\mathbf{X})\):
\[ \widehat{Y} = \bsy{\phi}\tran\mathbf{Z} \]
Again a linear problem — now in the features \(\mathbf{Z}\).
E.g., polynomial regression with \(\bsy{\psi}(X) = [1, X, X^2, \ldots, X^K]\tran\).
\[ \widehat{Y} = \phi_0 + \phi_1 X + \phi_2 X^2 + \ldots + \phi_K X^K \]
Affine mapping \(\;\widehat{Y} = aX + b\;\) — find \(a,b\) minimizing
\[ \min_{a,b} J(a,b) = \min_{a,b} \E{(Y - aX - b)^2} \]
Solve for the intercept \(b\)
\[ \frac{\partial J}{\partial b} = \E{-2\,(Y - aX - b)} = 0 \]
\[ \E{Y} - a\,\E{X} - b = 0 \]
\[ \boxed{\,b = m_Y - a\,m_X\,} \]
The line passes through the means \((m_X, m_Y)\).
Insert \(b = m_Y - a\,m_X\), so the residual is centered:
\[ Y - aX - b = (Y-m_Y) - a\,(X-m_X) \]
\[ \frac{\partial J}{\partial a} = \E{-2\,(X-m_X)\big[(Y-m_Y) - a\,(X-m_X)\big]} = 0 \]
\[ \E{(X-m_X)(Y-m_Y)} = a\,\E{(X-m_X)^2} \]
\[ \boxed{\,a_{\mathrm{MMSE}} = \frac{C_{XY}}{\sigma_X^2}\,} \qquad \boxed{\,b_{\mathrm{MMSE}} = m_Y - a_{\mathrm{MMSE}}\,m_X\,} \]
Only second-order moments are needed — the affine predictor is fixed by \(m_X, m_Y, \sigma_X^2, C_{XY}\).
Now \(\mathbf{X}\in\Real^K\) — a vector slope \(\mathbf{a}\) and intercept \(b\):
\[ \widehat{Y} = \mathbf{a}\tran\mathbf{X} + b, \qquad J(\mathbf{a},b) = \E{(Y - \mathbf{a}\tran\mathbf{X} - b)^2} \]
Intercept first — \(\partial J/\partial b = 0\), exactly as before:
\[ \boxed{\,b = m_Y - \mathbf{a}\tran\mathbf{m}_X\,} \]
The fit passes through the means \((\mathbf{m}_X, m_Y)\).
Center the variables and set \(\nabla_{\mathbf{a}} J = \mathbf{0}\):
\[ \E{(\mathbf{X}-\mathbf{m}_X)\big[(Y-m_Y) - \mathbf{a}\tran(\mathbf{X}-\mathbf{m}_X)\big]} = \mathbf{0} \quad\text{(orthogonality)} \]
\[ \underbrace{\E{(\mathbf{X}-\mathbf{m}_X)(\mathbf{X}-\mathbf{m}_X)\tran}}_{\mathbf{C}_{XX}}\,\mathbf{a} = \underbrace{\E{(\mathbf{X}-\mathbf{m}_X)(Y-m_Y)}}_{\mathbf{c}_{XY}} \]
\[ \boxed{\,\mathbf{a}_{\mathrm{MMSE}} = \mathbf{C}_{XX}^{-1}\,\mathbf{c}_{XY}\,} \]
The matrix version of \(a = C_{XY}/\sigma_X^2\).
We rarely know the joint distribution — only \(N\) data pairs \((\mathbf{x}_i, y_i)\). Replace the expectation by an empirical sum (linear model \(\widehat{Y} = \bsy{\phi}\tran\mathbf{X}\)):
\[ \bsy{\phi}_{\mathrm{LS}} = \argmin_{\bsy{\phi}} \sum_{i=1}^N (y_i - \bsy{\phi}\tran\mathbf{x}_i)^2 = \argmin_{\bsy{\phi}} \norm{\mathbf{y} - \bsy{\chi}\,\bsy{\phi}}^2, \qquad \bsy{\chi} = \begin{bmatrix} \mathbf{x}_1\tran \\ \vdots \\ \mathbf{x}_N\tran \end{bmatrix} \]
\[ \boxed{\,\bsy{\phi}_{\mathrm{LS}} = (\bsy{\chi}\tran\bsy{\chi})^{-1}\bsy{\chi}\tran\mathbf{y}\,} \]
By the law of large numbers, empirical moments converge to population moments:
\[ \frac{1}{N}\bsy{\chi}\tran\bsy{\chi} \;\to\; \mathbf{R}_{XX} = \E{\mathbf{X}\mathbf{X}\tran}, \qquad \frac{1}{N}\bsy{\chi}\tran\mathbf{y} \;\to\; \mathbf{R}_{XY} = \E{Y\mathbf{X}} \]
\[ \boxed{\,\bsy{\phi}_{\mathrm{LS}} \;\xrightarrow{\ N\to\infty\ }\; \mathbf{R}_{XX}^{-1}\mathbf{R}_{XY} = \bsy{\phi}_{\mathrm{MMSE}}\,} \]
Least squares is sample-based MMSE.
Both solutions invert \(\mathbf{R}_{XX}\). This is fragile when \(\mathbf{R}_{XX}\) is ill-conditioned:
Eigendecomposition \(\;\mathbf{R}_{XX} = \sum_i \lambda_i\,\mathbf{u}_i\mathbf{u}_i\tran\;\) (eigenvalues \(\lambda_i\)):
\[ \mathbf{R}_{XX}^{-1} = \sum_i \frac{1}{\lambda_i}\,\mathbf{u}_i\mathbf{u}_i\tran \]
A small \(\lambda_i \to 0\) blows up the inverse — a huge, noise-sensitive \(\bsy{\phi}\).
Add a small constant \(\lambda > 0\) to the diagonal:
\[ \boxed{\,\bsy{\phi}_\lambda = (\mathbf{R}_{XX} + \lambda\mathbf{I})^{-1}\mathbf{R}_{XY}\,} \]
Each eigenvalue is lifted away from zero, \(\;\lambda_i \mapsto \lambda_i + \lambda\):
\[ (\mathbf{R}_{XX} + \lambda\mathbf{I})^{-1} = \sum_i \frac{1}{\lambda_i + \lambda}\,\mathbf{u}_i\mathbf{u}_i\tran \]
The \(1/\lambda_i\) blow-up is capped at \(1/\lambda\). But where does \(\lambda\) come from?
So far: fit a mapping to data. Now: model how the data is generated.
Suppose the truth is linear with Gaussian observation noise:
\[ Y = \bsy{\phi}\tran\mathbf{X} + N, \qquad N \sim \Gauss(0, \sigma_n^2), \qquad N \perp \mathbf{X} \]
Then the conditional law of \(Y\) is Gaussian — symmetric and unimodal:
\[ f_{Y\mid X}(y\mid\mathbf{x}) = \Gauss\!\big(\bsy{\phi}\tran\mathbf{x},\ \sigma_n^2\big) \;\Rightarrow\; \E{Y\mid\mathbf{x}} = \argmax_y f_{Y\mid X}(y\mid\mathbf{x}) = \bsy{\phi}\tran\mathbf{x} \]
The linear predictor is the best predictor — MMSE and MAP coincide. Restricting to linear loses nothing.
Now treat the slope as random, with a Gaussian prior:
\[ \bsy{\phi} \sim \Gauss(\mathbf{0},\ \sigma_\phi^2\mathbf{I}) \]
MAP maximizes posterior \(\propto\) likelihood \(\times\) prior, i.e. minimizes
\[ \frac{1}{\sigma_n^2}\,\E{(Y - \bsy{\phi}\tran\mathbf{X})^2} \;+\; \frac{1}{\sigma_\phi^2}\,\norm{\bsy{\phi}}^2 \]
\[ \boxed{\,\bsy{\phi}_{\mathrm{MAP}} = (\mathbf{R}_{XX} + \lambda\mathbf{I})^{-1}\mathbf{R}_{XY}, \qquad \lambda = \frac{\sigma_n^2}{\sigma_\phi^2}\,} \]
Tikhonov is a Gaussian prior on \(\bsy{\phi}\). Noisy data (large \(\sigma_n^2\)) or a tight prior (small \(\sigma_\phi^2\)) \(\Rightarrow\) more shrinkage. → Demo: Conditioning & Ridge