Gaussian Linear Model

A Gaussian linear model is a linear model with Gaussian noise:

$$Y_i=X_i^T\theta +{ϵ}_i,{ϵ}_i\stackrel{\text{iid}}{\sim }N(0,{\sigma }^2).$$

Suppose we have a sample size $n$, $X_i\in {\mathbb{R}}^d$, $\theta \in {\mathbb{R}}^d$, $Y_i\in \mathbb{R}$. This note focuses on the regression task with a fixed design matrix $X=(X_1,\dots ,X_n{)}^T\in {\mathbb{R}}^{n\times d}$, which reduces to the Estimation of $\theta$.

Least Squares and Maximum Likelihood Estimation

Ordinary Least Squares gives the Maximum Likelihood Estimation of $\theta$:

$${\hat{\theta }}^{(\mathrm{MLE})}={\hat{\theta }}^{(\mathrm{LS})}=(X^{T}X{)}^{-1}{X}^{T}Y.$$

We can show that

$${\hat{\theta }}^{(\mathrm{LS})}\sim \mathcal{N}(\theta ,{\sigma }^2(X^{T}X{)}^{-1}{I}_n).$$

See A Probabilistic View Maximum Likelihood Estimation for the derivation.

Bayes Estimator

For prior $Q\sim \mathcal{N}(0,{\tau }^2{I}_d)$ and a Bowl-Shaped Loss $L(a,\theta )=\ell (a-\theta )$, the Bayes Optimal Estimator is

$${\hat{\theta }}^{(\mathrm{Bayes})}={\left(X^{T}X+\frac{{\sigma }^2}{{\tau }^2}I\right)}^{-1}{X}^{T}Y.$$

1st Proof

We note that this is essentially Bayesian Linear Regression. Specifically, the posterior is also a Gaussian (Conjugate Prior) that satisfies

$$\begin{array}{rl}p(\theta |Y)& \propto p(\theta ,Y)\propto \exp \left(-\frac{1}{2}((Y-X\theta {)}^T{\sigma }^{-2}(Y-X\theta )+{\theta }^T{\tau }^{-2}\theta )\right)\\ & \propto \exp \left(-\frac{1}{2}(\theta -{\hat{\theta }}_{\mathrm{post}}{)}^T{\hat{\Sigma }}_{\mathrm{post}}^{-1}(\theta -{\hat{\theta }}_{\mathrm{post}})\right)\\ & \sim \mathcal{N}({\hat{\theta }}_{\mathrm{post}},{\hat{\Sigma }}_{\mathrm{post}})\end{array}$$

where ${\hat{\theta }}_{\mathrm{post}}=(X^{T}X+\frac{{\sigma }^2}{{\tau }^2}I{)}^{-1}{X}^{T}Y$, ${\hat{\Sigma }}_{\mathrm{post}}=({\sigma }^{-2}{X}^{T}X+{\tau }^{-2}I{)}^{-1}$

By Anderson’s Lemma, ${\hat{\theta }}^{(\mathrm{Bayes})}={\hat{\theta }}_{\mathrm{post}}=(X^{T}X+\frac{{\sigma }^2}{{\tau }^2}I{)}^{-1}{X}^{T}Y$, Moreover, the Bayes risk is $R_B^{\ast }(\Theta )=\mathbb{E}[\ell (X_{\mathrm{post}}^{1/2}W)]$ where $W\sim \mathcal{N}(0,I)$.

2nd Proof

Similarly, by Anderson’s Lemma, we know that ${\hat{\theta }}^{(\mathrm{Bayes})}$ is the posterior mean. Since for normal distribution, mean and mode coincide, we have

$$\begin{array}{rl}{\hat{\theta }}_{\mathrm{post}}^{\mathrm{mean}}={\hat{\theta }}_{\mathrm{post}}^{\mathrm{mode}}& =\arg \underset{\theta }{\max }\exp \left(-\frac{1}{2}((Y-X\theta {)}^T{\sigma }^{-2}(Y-X\theta )+{\theta }^T{\tau }^{-2}\theta )\right)\\ & =\arg \underset{\theta }{\min }({\sigma }^{-2}||Y-X\theta |{|}^2+{\tau }^{-2}||\theta |{|}^2)\\ & =\arg \underset{\theta }{\min }(||Y-X\theta |{|}^2+\frac{{\sigma }^2}{{\tau }^2}||\theta |{|}^2).\end{array}$$

This corresponds to a Ridge Regression, whose solution is

$${\hat{\theta }}^{(\mathrm{Bayes})}={\hat{\theta }}_{\mathrm{post}}^{\mathrm{mean}}={\hat{\theta }}_{\mathrm{post}}^{\mathrm{mode}}={\hat{\theta }}^{(\mathrm{ridge})}={\left(X^{T}X+\frac{{\sigma }^2}{{\tau }^2}I\right)}^{-1}{X}^{T}Y.$$

Minimax Estimator

For a Bowl-Shaped Loss $L(a,\theta )=\ell (a-\theta )$ and a full column rank design matrix $X\in {\mathbb{R}}^{n\times d}$, Ordinary Least Squares also gives the Minimax Optimal Estimator of $\theta$:

$${\hat{\theta }}^{(\mathrm{Minimax})}={\hat{\theta }}^{(\mathrm{LS})}.$$

Proof

Recall that ${\hat{\theta }}^{(\mathrm{LS})}\sim \mathcal{N}(\theta ,(X^{T}X{)}^{-1}{\sigma }^2)$. Thus, the risk is

$$R({\hat{\theta }}^{(\mathrm{LS})})=\mathbb{E}[\ell ((X^{T}X/{\sigma }^2{)}^{-1}W)],W\sim \mathcal{N}(0,I_d).$$

Case I. $d=n$ and $X=I$. Note that the risk of the least squares estimator is independent of $\theta$, and thus to show it’s minimax optimal, we aim to find a prior $Q$ whose Bayes risk matches $R_M({\hat{\theta }}^{(\mathrm{LS})})=R({\hat{\theta }}^{(\mathrm{LS})})$ (see Minimax via Bayes). Recall (see Bayes Estimator) that given a normal prior $Q_{\tau }=\mathcal{N}(0,{\tau }^2{I}_d)$, the posterior is $\mathcal{N}({\hat{\theta }}_{\mathrm{post}},(\frac{X^{T}X}{{\sigma }^2}+\frac{I}{{\tau }^2}{)}^{-1})$, and the Bayes risk w.r.t $Q_{\tau }$ is

$$R_B^{\ast }(Q_{\tau })=\mathbb{E}\left[\ell \left({\left(\frac{X^{T}X}{{\sigma }^2}+\frac{I}{{\tau }^2}\right)}^{-1/2}W\right)\right]=\mathbb{E}\left[\ell \left({\left({\sigma }^{-2}+{\tau }^{-2}\right)}^{-1/2}W\right)\right]$$

Since $\ell$ is convex hence continuous, by the continuous mapping theorem,

$$\begin{array}{rl}\underset{\tau \to \infty }{\lim }{R}_B^{\ast }(Q_{\tau })=& \mathbb{E}\left[\underset{\tau \to \infty }{\lim }\ell (({\sigma }^{-2}+{\tau }^{-2}{)}^{-1/2}W)\right]\\ =& \mathbb{E}\left[\ell (\underset{\tau \to \infty }{\lim }({\sigma }^{-2}+{\tau }^{-2}{)}^{-1/2}W)\right]\\ =& \mathbb{E}\left[\ell (({\sigma }^{-2}{)}^{-1/2}W)\right]\\ =& \mathbb{E}\left[\ell \left({\left(\frac{X^{T}X}{{\sigma }^2}\right)}^{-1/2}W\right)\right]\\ =& R({\hat{\theta }}^{(\mathrm{LS})}).\end{array}$$

Case II. $d=n$ and $X\ne I$. We define a new loss function by $\tilde{\ell }(a)=\ell (X^{-1}a)$. One can check that $\tilde{\ell }$ is also bowl-shaped. For a fixed design matrix $X$, we consider a general estimator $\hat{\theta }(Y)$ determined by $Y$. We have the equivalence:

$$\begin{array}{rl}R(\hat{\theta }(Y))=& {\mathbb{E}}_Y\ell (\hat{\theta }(Y)-\theta )\\ =& {\mathbb{E}}_Y\tilde{\ell }(X\hat{\theta }(Y)-X\theta )\\ =& {\mathbb{E}}_Y\tilde{\ell }(X\hat{\theta }(X\theta +ϵ)-X\theta )\\ =& {\mathbb{E}}_Y\tilde{\ell }(\hat{\eta }(S)-\eta ),\end{array}$$

where the last equation applies the notation change:

$$\eta =X\theta ,S=I\eta +ϵ,\hat{\eta }=X\hat{\theta }.$$

In words, $\hat{\eta }$ is an estimator of the new parameter $\eta$ based on data $S$. One can see that

$${\hat{\eta }}^{(\mathrm{LS})}=S=Y=X^{-T}{X}^{T}X(X^{T}X{)}^{-1}{X}^{T}Y=X{\hat{\theta }}^{(\mathrm{LS})}.$$

We denote $\tilde{R}$ the risk of an estimator of $\eta$ w.r.t the new loss $\tilde{\ell }$. Then, applying Case I to ${\hat{\eta }}^{(\mathrm{LS})}$ gives

$$\tilde{R}({\hat{\eta }}^{(\mathrm{LS})})={\tilde{R}}_M({\hat{\eta }}^{(\mathrm{LS})})\le {\tilde{R}}_M(X\hat{\theta })=\underset{\eta }{\sup }R(X\hat{\theta },\eta ),$$

where $\hat{\theta }$ is a general estimator of $\theta$. Since $X$ is invertible, we have

$$\underset{\eta }{\sup }\tilde{R}(X\hat{\theta },\eta )=\underset{\theta }{\sup }\tilde{R}(X\hat{\theta },X\theta )=\underset{\theta }{\sup }R(\hat{\theta },\theta )=R_M(\hat{\theta }).$$

On the other hand, we have

$$\tilde{R}({\hat{\eta }}^{(\mathrm{LS})})=\mathbb{E}\tilde{\ell }(\sigma W)=\mathbb{E}\ell (\sigma X^{-1}W)=\mathbb{E}\left[\ell \left({\left(\frac{X^{T}X}{{\sigma }^2}\right)}^{-1/2}W\right)\right],$$

where the last inequality uses the fact that ${\sigma }^{-1}{X}^{-1}W\sim \mathcal{N}(0,(X^{T}X{)}^{-1}{\sigma }^2)$. Combining the above three inequalities gives

$$R({\hat{\theta }}^{(\mathrm{LS})})=\tilde{R}({\hat{\eta }}^{(\mathrm{LS})})\le R_M(\hat{\theta }),\forall \hat{\theta }.$$

Thus, ${\hat{\theta }}^{(\mathrm{LS})}$ is minimax optimal.

Case III. $d<n$. Let $U=X(X^{T}X{)}^{-1/2}$. Then $U{U}^T=X(X^{T}X{)}^{-1}{X}^T$ is the orthogonal projection onto the column space of $X$. Then, we have

$$U^{T}Y\sim \mathcal{N}((X^{T}X{)}^{1/2}\theta ,{\sigma }^2{U}^{T}U)=\mathcal{N}((X^{T}X{)}^{1/2}\theta ,{\sigma }^2{I}_d).$$

Equivalently, the original data gives a new Gaussian linear model:

$$U^{T}Y=(X^{T}X{)}^{1/2}\theta +{ϵ}^{\prime },{ϵ}^{\prime }\sim \mathcal{N}(0,{\sigma }^2{I}_d).$$

This Gaussian linear model reduces to Case II. Specifically, with a fixed design matrix $(X^{T}X{)}^{1/2}$, let $A:U^{T}Y\mapsto \hat{\theta }$ be a general estimator of $\theta$. And we denote $R_d$ the risk w.r.t this new Gaussian linear model. Then, Case II gives

$$\underset{\theta }{\sup }{R}_d(A,\theta )\ge \mathbb{E}\left[\ell \left({\left(\frac{((X^{T}X{)}^{1/2}{)}^T(X^{T}X{)}^{1/2}}{{\sigma }^2}\right)}^{-1/2}W\right)\right]=\mathbb{E}\left[\ell \left({\left(\frac{X^{T}X}{{\sigma }^2}\right)}^{-1/2}W\right)\right]=R({\hat{\theta }}^{(\mathrm{LS})}).$$

Therefore, we are left to show that for any general estimator $\hat{\theta }(Y)$ corresponding to the original Gaussian linear model, there exists a induced estimator $A(U^{T}Y)$ such that $R_d(A,\theta )=R(\hat{\theta },\theta )$. We claim this is true with the following randomized induced estimator:

$$A(T)=\hat{\theta }(Y),Y\sim \cdot \mid T,$$

where $T=U^{T}Y$.

The induced estimator works as follows: upon observing $T$, we simulate $Y$ from the conditional distribution $\cdot \mid T$, and then apply the original estimator $\hat{\theta }(Y)$ to obtain $A(T)$. In our case, since we actually observe $Y$, we directly have $A(T)=\hat{\theta }(Y)$. However, we can still interpret $Y$ as being generated from the conditional distribution $\cdot \mid T$.

Suppose $T$ is a Sufficient Statistic of $Y$, the above claim is true:

$$\begin{array}{rl}R(\hat{\theta },\theta )=& {\mathbb{E}}_{Y\sim P_{\theta }}\ell (\hat{\theta }-\theta )\\ =& {\mathbb{E}}_{Y\sim \cdot \mid T}\ell (\hat{\theta }-\theta )\\ =& {\mathbb{E}}_T[{\mathbb{E}}_Y[\ell (\hat{\theta }-\theta )\mid T]]\\ =& {\mathbb{E}}_T[{\mathbb{E}}_Y[\ell (A(T)-\theta )\mid T]]\\ =& \mathbb{E}\ell (A(T)-\theta )\\ =& R_d(A,\theta ).\end{array}$$

Therefore, we are left to show that $T$ is indeed a sufficient statistic of $Y$. We provide three methods.

Method I. Intuition

Let $P:=U{U}^T\in \mathrm{OrthBasis}(\mathrm{col}(X))$. Then,

$$X\theta +ϵ=\underset{\in \mathrm{col}(X)}{\underbrace{X\theta +Pϵ}}+\underset{\in \mathrm{col}(X{)}^{\perp }}{\underbrace{(I-P)ϵ}}$$

Note that for a Gaussian vector, zero correlation implies independence. And since perpendicular vectors have zero correlation, we have $(I-P)ϵ\perp \perp U^{T}Y=X\theta +Pϵ$. Moreover, since $(I-P)ϵ\perp \perp \theta$, we can simulate $Y=U{U}^{T}Y+(I-P)ϵ$ without $\theta$. Thus $U^{T}Y$ is sufficient.

In words, $U^{T}Y$ captures all the information of $\theta$ left-applied by $X$. The remaining part is pure noise perpendicular to the column space of $X$ and does not depend on $\theta$. This is illustrated in the following plot, where ${ϵ}_1=Pϵ$ and ${ϵ}_2=(I-P)ϵ$.

Decomposition of Z\theta and \epsilon_{2}.

Method II. Fisher-Neyman Factorization

We first show that $X^{T}Y$ is sufficient:

$$\begin{array}{rl}Y\sim \mathcal{N}(X\theta ,{\sigma }^{2}I)⟹p(y)& =(2\pi {\sigma }^2{)}^{-\frac{n}{2}}\exp \left(-\frac{1}{2{\sigma }^2}(y-X\theta {)}^T(y-X\theta )\right)\\ & =\underset{g(T(y),\theta )}{\underbrace{(2\pi {\sigma }^2{)}^{-\frac{n}{2}}\exp (-\frac{1}{2{\sigma }^2}({\theta }^T{X}^{T}X\theta -2{\theta }^T{X}^{T}y}}+y^{T}y))\end{array}$$

By Fisher-Neyman Factorization Theorem, $X^{T}Y\in {\mathbb{R}}^d$ is sufficient Since $X$ is full column-rank, $(X^{T}X{)}^{-1/2}$ exists and $(X^{T}X{)}^{-1/2}{X}^{T}Y$ is also sufficient.

Method III. Conditional Simulation

We first show that $Y\mid U{U}^{T}Y\perp \perp \theta$.

Note that

$$\begin{array}{rl}U{U}^{T}Y& =X(X^{T}X{)}^{-1}{X}^T(X\theta +ϵ)\\ & =X\theta +(U{U}^Tϵ)\\ & =Y-(I-U{U}^T)ϵ\end{array}$$

Therefore,

$$Y\mid U{U}^{T}Y=t\sim \mathcal{N}(t,(I-U{U}^T){\sigma }^2(I-U{U}^T{)}^T).$$

Finally,

$$Y\mid U^{T}Y=t\sim \mathcal{N}(Ut,(I-U{U}^T){\sigma }^{2}I)$$

Thus, $U^{T}Y$ is sufficient.