Constrained Gaussian Mean Estimation Under Affine Restrictions

Why Constraints Change the Estimate

In a multivariate Gaussian model, an unrestricted estimate of the mean is the sample mean. But scientific knowledge may require its components to satisfy exact linear relationships: two population means may be equal, a weighted sum may be fixed, or several means may obey conservation equations. The estimate must then stay within the set of means allowed by those restrictions.

The best constrained value is not generally the closest allowed vector in ordinary Euclidean distance. The likelihood measures discrepancies using the covariance matrix, so the appropriate projection depends on the scale and correlation of the measurements. This tutorial derives that covariance-weighted projection for a known covariance matrix.

Model and definition: Let $X_1,\ldots,X_n$ be independent random vectors with $X_i\sim N_p(\mu,\Sigma)$, where $n\geq 1$, $\mu\in\mathbb{R}^p$, and the known covariance matrix $\Sigma$ is symmetric positive definite. Suppose the parameter is restricted by $A\mu=b$, where $A$ is an $r\times p$ matrix of full row rank $r$, with $1\leq r\leq p$, and $b\in\mathbb{R}^r$. This system is feasible because $A$ has full row rank. The constrained maximum-likelihood estimate is the value of $\mu$ satisfying $A\mu=b$ that maximizes the likelihood.

From the Likelihood to a Constrained Distance

Write $\bar X=n^{-1}\sum_{i=1}^n X_i$ for the sample mean, and let $\bar x$ denote its observed value. With $\Sigma$ fixed, the part of the log-likelihood that depends on $\mu$ is

\ell(\mu)=C-\frac{1}{2}\sum_{i=1}^n (x_i-\mu)^{\mathsf T}\Sigma^{-1}(x_i-\mu) =C'-\frac{n}{2}(\bar x-\mu)^{\mathsf T}\Sigma^{-1}(\bar x-\mu),

where $C$ and $C'$ do not depend on $\mu$. The second expression follows by writing $x_i-\mu=(x_i-\bar x)+(\bar x-\mu)$: the cross terms sum to zero because $\sum_i(x_i-\bar x)=0$. Thus maximizing the likelihood over $A\mu=b$ is equivalent to minimizing the squared distance from $\bar x$ to $\mu$, measured by $\Sigma^{-1}$:

\text{minimize}_{\mu:\,A\mu=b}\quad (\bar x-\mu)^{\mathsf T}\Sigma^{-1}(\bar x-\mu).

Since $\Sigma^{-1}$ is positive definite, this objective is strictly convex in $\mu$. The constraint set is an affine subspace, so a minimizer, if found, is unique. The following result gives it explicitly.

Theorem: Under the model and rank conditions above, the constrained maximum-likelihood estimate is $$ \hat\mu_c=\bar x-\Sigma A^{\mathsf T}(A\Sigma A^{\mathsf T})^{-1}(A\bar x-b). $$ The matrix $A\Sigma A^{\mathsf T}$ is invertible, and the estimate satisfies $A\hat\mu_c=b$.

To see why the inverse exists, take any nonzero $v\in\mathbb{R}^r$. Full row rank of $A$ implies $A^{\mathsf T}v\neq 0$, and positive definiteness of $\Sigma$ gives

v^{\mathsf T}A\Sigma A^{\mathsf T}v =(A^{\mathsf T}v)^{\mathsf T}\Sigma(A^{\mathsf T}v)>0.

Therefore $A\Sigma A^{\mathsf T}$ is positive definite and invertible.

Deriving and Verifying the Estimate

Introduce a Lagrange multiplier $\lambda\in\mathbb{R}^r$ for the constraint. It is convenient to minimize one half of the quadratic objective; this does not change its minimizer. The Lagrangian is

Q(\mu,\lambda)=\frac12(\mu-\bar x)^{\mathsf T}\Sigma^{-1}(\mu-\bar x) +\lambda^{\mathsf T}(A\mu-b).

Differentiating with respect to $\mu$ and setting the gradient equal to zero gives

\Sigma^{-1}(\mu-\bar x)+A^{\mathsf T}\lambda=0, \qquad \mu=\bar x-\Sigma A^{\mathsf T}\lambda.

Substitute this expression into $A\mu=b$. Solving the resulting $r$-dimensional system determines $\lambda$, and hence $\mu$:

A\bar x-A\Sigma A^{\mathsf T}\lambda=b, \qquad \lambda=(A\Sigma A^{\mathsf T})^{-1}(A\bar x-b), $$ $$ \mu=\bar x-\Sigma A^{\mathsf T}(A\Sigma A^{\mathsf T})^{-1}(A\bar x-b).

This candidate satisfies the constraint, since

A\hat\mu_c =A\bar x-A\Sigma A^{\mathsf T}(A\Sigma A^{\mathsf T})^{-1}(A\bar x-b) =b.

It is the unique minimizer because the objective has positive-definite Hessian $\Sigma^{-1}$. More explicitly, for any nonzero $h\in\mathbb{R}^p$, $h^{\mathsf T}\Sigma^{-1}h>0$; thus the objective is strictly convex, and its stationary point on the feasible affine set is the unique constrained minimum. This proves the theorem. $\square$

Worked Example: One Restriction With Correlated Components

Suppose $X_i\sim N_2(\mu,\Sigma)$, where $\Sigma=\begin{pmatrix}4&1\\1&2\end{pmatrix}$ is known, and the observed sample mean is $\bar x=(3,1)^{\mathsf T}$. Impose the restriction $\mu_1+\mu_2=2$. Find the constrained estimate.

Here $A=(1,1)$ and $b=2$. First calculate the two quantities in the correction:

A\bar x-b=3+1-2=2, \qquad A\Sigma A^{\mathsf T}=4+1+1+2=8, \qquad \Sigma A^{\mathsf T}= \begin{pmatrix}5\\3\end{pmatrix}.

The formula therefore gives

\hat\mu_c= \begin{pmatrix}3\\1\end{pmatrix} -\begin{pmatrix}5\\3\end{pmatrix}\frac{2}{8} = \begin{pmatrix}3\\1\end{pmatrix} -\begin{pmatrix}5/4\\3/4\end{pmatrix} = \begin{pmatrix}7/4\\1/4\end{pmatrix}.

The estimate obeys the restriction: $7/4+1/4=2$. The two components are adjusted by different amounts because the covariance-weighted geometry is not ordinary Euclidean geometry.

Worked Example: A Diagonal Covariance and a Weighted Equality

Let $\Sigma=\operatorname{diag}(1,9)$, let $\bar x=(4,1)^{\mathsf T}$, and require $\mu_1-2\mu_2=0$. Thus $A=(1,-2)$, $b=0$, and the restriction is $\mu_1=2\mu_2$.

Compute the discrepancy and the covariance terms:

A\bar x-b=4-2(1)=2, \qquad A\Sigma A^{\mathsf T}=1+4(9)=37, \qquad \Sigma A^{\mathsf T}= \begin{pmatrix}1\\-18\end{pmatrix}.

Hence

\hat\mu_c= \begin{pmatrix}4\\1\end{pmatrix} -\begin{pmatrix}1\\-18\end{pmatrix}\frac{2}{37} = \begin{pmatrix}146/37\\73/37\end{pmatrix}.

The first component is twice the second, since $146/37=2(73/37)$. The adjustment is $(-2/37,36/37)^{\mathsf T}$. In particular, the component with larger variance moves more in this example; the restriction is met at the likelihood-optimal cost, not by minimizing the unweighted size of the adjustment.

Projection Geometry and the Cost of the Restriction

The estimate is a projection in the inner product defined by $\Sigma^{-1}$: the squared length of a vector $u$ is $u^{\mathsf T}\Sigma^{-1}u$. This choice of metric follows directly from the Gaussian likelihood. A useful way to confirm the projection interpretation is to compare the constrained estimate with any other feasible value.

Result: For every $\mu$ satisfying $A\mu=b$, let $d=\mu-\hat\mu_c$. Then $$ (\bar x-\mu)^{\mathsf T}\Sigma^{-1}(\bar x-\mu) = (\bar x-\hat\mu_c)^{\mathsf T}\Sigma^{-1}(\bar x-\hat\mu_c) +d^{\mathsf T}\Sigma^{-1}d. $$ Consequently, every other feasible value has an objective at least as large, with equality only at $\mu=\hat\mu_c$.

To prove the result, both $\mu$ and $\hat\mu_c$ satisfy the restriction, so $Ad=0$. The formula for the estimate also gives

\bar x-\hat\mu_c =\Sigma A^{\mathsf T}(A\Sigma A^{\mathsf T})^{-1}(A\bar x-b).

Expanding $\bar x-\mu=(\bar x-\hat\mu_c)-d$, the difference between the two sides' quadratic terms is the extra term $d^{\mathsf T}\Sigma^{-1}d$ plus a cross term. That cross term vanishes:

(\bar x-\hat\mu_c)^{\mathsf T}\Sigma^{-1}d =(A\bar x-b)^{\mathsf T}(A\Sigma A^{\mathsf T})^{-1}Ad =0.

The displayed identity follows. Since $\Sigma^{-1}$ is positive definite, $d^{\mathsf T}\Sigma^{-1}d\geq 0$, and it is zero only when $d=0$. This proves both minimality and uniqueness. $\square$

The same identity quantifies the increase in the minimized quadratic distance. The unrestricted estimate is $\bar x$, so the constrained estimate pays the distance

(\bar x-\hat\mu_c)^{\mathsf T}\Sigma^{-1}(\bar x-\hat\mu_c) =(A\bar x-b)^{\mathsf T}(A\Sigma A^{\mathsf T})^{-1}(A\bar x-b).

The expression is nonnegative and is zero exactly when the sample mean already satisfies the restriction. Because the log-likelihood contains $-n/2$ times the quadratic distance, the corresponding loss in maximized log-likelihood relative to the unrestricted fit is $n/2$ times this quantity.

Worked Example: Two Restrictions in Three Dimensions

Let $\Sigma=\operatorname{diag}(1,2,1)$, $\bar x=(2,0,1)^{\mathsf T}$, and impose $\mu_1+\mu_3=0$ and $\mu_2+\mu_3=0$. Use $A=\begin{pmatrix}1&0&1\\0&1&1\end{pmatrix}$ and $b=(0,0)^{\mathsf T}$.

The discrepancy is $A\bar x-b=(3,1)^{\mathsf T}$. Direct multiplication gives

A\Sigma A^{\mathsf T} =\begin{pmatrix}2&1\\1&3\end{pmatrix}, \qquad (A\Sigma A^{\mathsf T})^{-1} =\frac15\begin{pmatrix}3&-1\\-1&2\end{pmatrix}, $$ $$ (A\Sigma A^{\mathsf T})^{-1}(A\bar x-b) =\frac15\begin{pmatrix}8\\-1\end{pmatrix}, \qquad \Sigma A^{\mathsf T}= \begin{pmatrix}1&0\\0&2\\1&1\end{pmatrix}.

Multiplying the last two matrices gives the adjustment $(8/5,-2/5,7/5)^{\mathsf T}$. Therefore

\hat\mu_c= \begin{pmatrix}2\\0\\1\end{pmatrix} -\begin{pmatrix}8/5\\-2/5\\7/5\end{pmatrix} = \begin{pmatrix}2/5\\2/5\\-2/5\end{pmatrix}.

Both restrictions hold: $2/5-2/5=0$ and $2/5-2/5=0$. This example also shows that the constraints can be handled together using the $r\times r$ matrix $A\Sigma A^{\mathsf T}$, even when the mean vector has more components.

Common Pitfall and Interpretation

A common mistake is to use the ordinary Euclidean projection of $\bar x$ onto $A\mu=b$, effectively replacing $\Sigma$ by the identity matrix. That works only when the covariance metric is proportional to the Euclidean metric, or in special cases where both projections happen to coincide. In general, ignoring $\Sigma$ ignores how the Gaussian likelihood scales discrepancies and accounts for correlations.

Another error is to apply the restriction to the observed sample mean and then adjust components by intuition. The formula shows that the adjustment has a specific direction, $\Sigma A^{\mathsf T}$, and a size determined by $A\Sigma A^{\mathsf T}$. It both enforces the restriction and minimizes the likelihood-based distance. The assumptions matter: this result treats $\Sigma$ as known and positive definite. If the covariance is also unknown, the likelihood optimization changes; this formula should not be used without checking the model.

Key takeaway: With known positive-definite covariance and full-row-rank affine restrictions, the constrained Gaussian mean MLE is the covariance-weighted projection of $\bar x$ onto the feasible set. Its adjustment is $\Sigma A^{\mathsf T}(A\Sigma A^{\mathsf T})^{-1}(A\bar x-b)$, and its uniqueness follows from strict convexity.

Check Your Understanding

Use the model assumptions stated above to answer these questions.

Why must $A\Sigma A^{\mathsf T}$ be invertible under the stated rank and covariance conditions?
For a single restriction $a^{\mathsf T}\mu=b$, identify the direction in which the sample mean is adjusted.
When is the constrained estimate equal to the unrestricted sample mean?
Explain why a Euclidean projection need not maximize the Gaussian likelihood when $\Sigma$ is not proportional to the identity matrix.
What is the increase in the minimized quadratic distance caused by imposing the restriction?