Deriving Weighted Least Squares from the Gaussian Score

From Gaussian Precision Weights to Least Squares

In Gaussian regression, observations need not contribute equally to the fit. If an observation has a larger precision weight, its variance is smaller under the model, and a given residual contributes more to the log-likelihood. This tutorial makes that connection precise: differentiating the Gaussian log-likelihood gives the weighted normal equations, and solving those equations minimizes the weighted residual sum of squares.

Let $Y_1,\ldots,Y_n$ be conditionally independent observations, with observed values $y_1,\ldots,y_n$. Let $X$ be the $n\times p$ design matrix, with row $i$ equal to $x_i^\mathsf{T}$, and let $\beta\in\mathbb{R}^p$ be the coefficient vector. Assume the identity-link Gaussian model and positive fixed precision weights $w_i$:

Y_i\mid X \sim N(x_i^\mathsf{T}\beta,\phi/w_i), \qquad w_i>0,\quad \phi>0.

Here $\phi$ is the common dispersion, and $w_i$ multiplies the precision $1/\phi$. Define $W=\operatorname{diag}(w_1,\ldots,w_n)$ and the weighted residual sum of squares

Q(\beta)=(y-X\beta)^\mathsf{T}W(y-X\beta) =\sum_{i=1}^n w_i(y_i-x_i^\mathsf{T}\beta)^2.

Definition: The weighted normal equations for this model are $$ X^\mathsf{T}WX\hat{\beta}=X^\mathsf{T}Wy. $$ They are the score equations for the Gaussian regression coefficients. Any solution minimizes $Q(\beta)$; if $X$ has full column rank, the solution is unique.

Deriving the Score and the Normal Equations

For fixed $\phi$ and weights, the independent Gaussian likelihood is the product of the observation densities. Taking its logarithm gives

\ell(\beta,\phi) =-\frac{n}{2}\log(2\pi\phi) +\frac{1}{2}\sum_{i=1}^n\log w_i -\frac{1}{2\phi}\sum_{i=1}^n w_i(y_i-x_i^\mathsf{T}\beta)^2.

The first two terms do not depend on $\beta$. Differentiating the final term with respect to the vector $\beta$, and writing $r=y-X\beta$, yields the score

U(\beta)=\frac{\partial \ell}{\partial\beta} =\frac{1}{\phi}X^\mathsf{T}W(y-X\beta).

The likelihood score equation $U(\beta)=0$ is therefore equivalent to $X^\mathsf{T}W(y-X\beta)=0$, since $\phi>0$. Rearranging gives the weighted normal equations. The derivation does not require estimating $\phi$: for any fixed positive $\phi$, the factor $1/\phi$ does not change the equation for $\beta$.

Result: Under the independent Gaussian model with variances $\phi/w_i$, the coefficient score vanishes exactly when $$ X^\mathsf{T}WX\beta=X^\mathsf{T}Wy. $$ Thus the Gaussian score weights each residual by its precision weight before setting the coefficient score to zero.

Worked Example: A Weighted Intercept

Suppose three independent observations have means equal to a common intercept $\beta$, with observed values $2,4,10$ and precision weights $1,2,1$, respectively. Find the score estimate.

For an intercept-only model, $X$ is a column of ones. The normal equation reduces to $\sum_i w_i(y_i-\beta)=0$. Equivalently,

(1+2+1)\hat{\beta}=1(2)+2(4)+1(10)=20, \qquad \hat{\beta}=5.

Check the score directly: the fitted residuals are $2-5=-3$, $4-5=-1$, and $10-5=5$. Their weighted sum is

$$ 1(-3)+2(-1)+1(5)=0. $$

The weighted residual sum of squares at this estimate is $1(9)+2(1)+1(25)=36$. The weight of 2 makes the middle observation count twice as much as an observation with weight 1 in the score and in $Q$.

Worked Example: A Weighted Line with an Intercept

Consider observations $(x_i,y_i)=(0,1),(1,2),(2,2)$, with weights $1,2,1$. Fit the mean $\beta_0+\beta_1x_i$. The design matrix has rows $(1,x_i)$. The entries of $X^\mathsf{T}WX$ are the weighted sums of $1,x_i,x_i^2$:

X^\mathsf{T}WX= \begin{pmatrix} 4&4\\ 4&6 \end{pmatrix}, \qquad X^\mathsf{T}Wy= \begin{pmatrix} 1(1)+2(2)+1(2)\\ 1(0)(1)+2(1)(2)+1(2)(2) \end{pmatrix} = \begin{pmatrix} 7\\8 \end{pmatrix}.

The equations are $4\beta_0+4\beta_1=7$ and $4\beta_0+6\beta_1=8$. Subtracting the first equation from the second gives $2\beta_1=1$, so $\hat{\beta}_1=0.5$, and then $4\hat{\beta}_0+2=7$, giving $\hat{\beta}_0=1.25$.

The fitted values are $1.25,1.75,2.25$, and the residuals are $-0.25,0.25,-0.25$. The two score components are proportional to the weighted residual sum and the weighted $x$-times-residual sum:

1(-0.25)+2(0.25)+1(-0.25)=0, \qquad 1(0)(-0.25)+2(1)(0.25)+1(2)(-0.25)=0.

The weighted residual sum of squares is $1(0.25^2)+2(0.25^2)+1(0.25^2)=0.25$, consistent with the score equations.

Why the Equations Minimize Weighted Residual Sum of Squares

A zero score gives a stationary point, but a stationary point need not minimize every possible objective. For this Gaussian model, the quadratic form of $Q$ allows a direct global argument.

Theorem: Suppose $W$ is diagonal with strictly positive entries. If $\hat{\beta}$ solves $X^\mathsf{T}WX\hat{\beta}=X^\mathsf{T}Wy$, then $Q(\hat{\beta})\leq Q(\beta)$ for every $\beta\in\mathbb{R}^p$. If, in addition, $X$ has full column rank, this minimizer is unique and equals $(X^\mathsf{T}WX)^{-1}X^\mathsf{T}Wy$.

Proof. Put $\delta=\beta-\hat{\beta}$. Since $y-X\beta=(y-X\hat{\beta})-X\delta$, expansion gives

Q(\beta)=Q(\hat{\beta}) -2\delta^\mathsf{T}X^\mathsf{T}W(y-X\hat{\beta}) +\delta^\mathsf{T}X^\mathsf{T}WX\delta.

The middle term is zero by the normal equations. The remaining difference is

Q(\beta)-Q(\hat{\beta}) =\delta^\mathsf{T}X^\mathsf{T}WX\delta =(X\delta)^\mathsf{T}W(X\delta) =\sum_{i=1}^n w_i(x_i^\mathsf{T}\delta)^2\geq 0.

This proves global minimization. If $X$ has full column rank and $\delta\ne0$, then $X\delta\ne0$. Because every $w_i>0$, the displayed sum is strictly positive. Hence no other $\beta$ attains the same minimum. Also, for any nonzero vector $v$, $v^\mathsf{T}X^\mathsf{T}WXv=\sum_i w_i(x_i^\mathsf{T}v)^2>0$; thus $X^\mathsf{T}WX$ is invertible, and solving the normal equations gives the stated formula. $\square$

The rank condition matters. If columns of $X$ are linearly dependent, different coefficient vectors can give the same fitted means and the same value of $Q$. The normal equations still describe the minimizers, but the coefficient vector need not be unique. This is why the inverse formula must not be used without checking full column rank.

Worked Example: A Weighted Regression Through the Origin

Suppose a model has no intercept, so its mean is $\beta x_i$. The observations are $(x_i,y_i)=(1,2),(2,3),(3,5)$, with weights $1,3,2$. There is one coefficient, and the normal equation is

\left(\sum_i w_ix_i^2\right)\hat{\beta} =\sum_i w_ix_iy_i.

The two sums are

\sum_i w_ix_i^2=1(1^2)+3(2^2)+2(3^2)=1+12+18=31, $$ $$ \sum_i w_ix_iy_i=1(1)(2)+3(2)(3)+2(3)(5)=2+18+30=50.

Therefore $\hat{\beta}=50/31$. The weighted score is proportional to $50-31\beta$, which is zero at this estimate. The fitted values are $50/31,100/31,150/31$, so the residuals are $12/31,-7/31,5/31$. Their weighted residual sum of squares is

Q(\hat{\beta}) =1(12/31)^2+3(-7/31)^2+2(5/31)^2 =\frac{144+147+50}{961} =\frac{341}{961}.

Since the design column is nonzero, it has full column rank in this one-coefficient model. The theorem therefore guarantees that this estimate is the unique global minimizer.

What the Result Does—and Does Not—Say

The derivation gives two equivalent ways to find the Gaussian coefficient estimate: set the coefficient score to zero, or minimize $Q(\beta)$. The equivalence follows because the log-likelihood depends on $\beta$ through $-Q(\beta)/(2\phi)$, with $\phi>0$. It is not necessary to expand the likelihood further to locate the coefficient estimate.

A common mistake is to assume that a larger precision weight divides the residual's contribution. Under the model used here, the variance is $\phi/w_i$, and the likelihood contains $w_i(y_i-x_i^\mathsf{T}\beta)^2$, not that residual squared divided by $w_i$. For example, replacing $W$ by $W^{-1}$ would change the normal equations and would correspond to a different objective. Always read the variance specification first: the weights in this derivation multiply the residual squares because they are precision weights.

Another error is to report $(X^\mathsf{T}WX)^{-1}X^\mathsf{T}Wy$ when $X^\mathsf{T}WX$ is singular. The normal equations remain the right condition, but an ordinary inverse then does not exist, and a unique coefficient vector is not guaranteed. For a full-column-rank design, positive weights ensure the inverse exists, as shown in the proof.

Key takeaway: With independent Gaussian responses having means $X\beta$ and variances $\phi/w_i$, the coefficient score is $X^\mathsf{T}W(y-X\beta)/\phi$. Its zero gives the weighted normal equations, and the quadratic expansion proves that their solutions minimize the weighted residual sum of squares.

Check Your Understanding

Use the model and arguments developed above to answer the following.

For independent Gaussian responses with variances $\phi/w_i$, write the coefficient score in matrix form and state the resulting normal equations.
In an intercept-only model, show how the normal equation expresses the estimate as a weighted average of the observed responses.
Why does the quadratic expansion of $Q(\beta)$ establish global minimization rather than only stationarity?
What condition on $X$ makes the coefficient estimate unique, and why do strictly positive weights matter?
If a stated variance is $\phi/w_i$, explain why the residual sum of squares uses $w_i$, rather than $1/w_i$.

Deriving Weighted Least Squares from the Gaussian Score

What You'll Learn