General Bayesian Estimators (MMSE, MAP, LMMSE)

Reference:
Kay S M. Fundamentals of statistical signal processing[M]. Prentice Hall PTR, 1993. (Ch. 11 - 11.5, Ch. 12 - 12.5)
Slides of ET4386, TUD

Motivation:

• The optimal Bayesian estimators discussed in the previous chapter are difficult to determine in closed form, and in practice too computationally intensive to implement.

• We now study more general Bayesian estimators and their properties. To do so, the concept of the Bayesian risk function is discussed.

• Minimization of this criterion results in a variety of estimators. The ones that we will concentrate on are the MMSE estimator and the maximum a posteriori estimator.

• They involve multidimensional integration for the MMSE estimator and multidimensional maximization for the MAP estimator. Although under the jointly Gaussian assumption these estimators are easily found, in general, they are not. When we are unable to make the Gaussian assumption, another approach must be used. To fill this gap we can choose to retain the MMSE criterion but constrain the estimator to be linear. Then, an explicit form for the estimator may be determined which depends only on the first two moments of the PDF. In many ways this approach is analogous to the BLUE in classical estimation, and some parallels will become apparent.

Risk Functions

Previously, we had derived the MMSE estimator by minimizing $E[(\theta -\hat{\theta }{)}^2]$, where the expectation is with respect to the PDF $p(\mathbf{x},\theta )$. If we let $ϵ=\theta -\hat{\theta }$ denote the error of the estimator for a particular realization of $\mathbf{x}$ and $\theta$, and also let $\mathcal{C}(ϵ)={ϵ}^2$, then the MSE criterion minimizes $E[\mathcal{C}(ϵ)]$. $\mathcal{C}(ϵ)$ is termed as cost function, and the average cost or $E[\mathcal{C}(ϵ)]$ is termed the Bayes risk $\mathcal{R}$ or
$$\mathcal{R}=E[\mathcal{C}(ϵ)]\text{\hspace{1em}(RF.1)}$$
If $\mathcal{C}(ϵ)={ϵ}^2$, then the cost function is quadratic and the Bayes risk is just the MSE. Of course, there is no need to restrict ourselves to quadratic cost functions. In Figure 11.1b we have the “absolute” error
$$\mathcal{C}(ϵ)=|ϵ|\text{\hspace{1em}(RF.2)}$$
In Figure 11.1c the “hit-or-miss” cost function is displayed
$$\mathcal{C}(ϵ)=\{\begin{array}{cc}0& |ϵ|<\delta \\ 1& |ϵ|>\delta \end{array}\text{\hspace{1em}(RF.3)}$$

We already have seen that the Bayes risk is minimized for a quadratic cost function by the MMSE estimator $\hat{\theta }=E(\theta |\mathbf{x})$. We now determine the optimal estimators for the other cost functions.

The Bayes risk $\mathcal{R}$ is
$$\begin{array}{rl}\mathcal{R}& =E[\mathcal{C}(ϵ)]\\ & =\int \int \mathcal{C}(\theta -\hat{\theta })p(\mathbf{x},\theta )d\mathbf{x}d\theta \\ & =\int \left[\int \mathcal{C}(\theta -\hat{\theta })p(\theta |\mathbf{x})d\theta \right]p(\mathbf{x})d\mathbf{x}\end{array}\text{\hspace{1em}(RF.4)}$$
We attempt to minimize the inner integral $g(\hat{\theta })=\int \mathcal{C}(\theta -\hat{\theta })p(\theta |\mathbf{x})d\theta$ for each $\mathbf{x}$.

• absolute error cost function
  $$\begin{array}{rl}g(\hat{\theta })& =\int |\theta -\hat{\theta }|p(\theta |\mathbf{x})d\theta \\ & ={\int }_{-\infty }^{\hat{\theta }}(\hat{\theta }-\theta )p(\theta |\mathbf{x})d\theta +{\int }_{\hat{\theta }}^{\infty }(\theta -\hat{\theta })p(\theta |\mathbf{x})d\theta \end{array}$$
  Differentiating w.r.t. $\hat{\theta }$ using Leibnitz’s rule
  $$\frac{\partial }{\partial u}{\int }_{{\varphi }_1(u)}^{{\varphi }_2(u)}h(u,v)dv={\int }_{{\varphi }_1(u)}^{{\varphi }_2(u)}\frac{\partial h(u,v)}{\partial u}dv+\frac{d{\varphi }_2(u)}{du}h(u,{\varphi }_2(u))-\frac{d{\varphi }_1(u)}{du}h(u,{\varphi }_1(u))$$
  we have
  $$\frac{dg(\hat{\theta })}{d\hat{\theta }}={\int }_{-\infty }^{\hat{\theta }}p(\theta |\mathbf{x})d\theta -{\int }_{\hat{\theta }}^{\infty }p(\theta |\mathbf{x})d\theta =0$$
  or
  $${\int }_{-\infty }^{\hat{\theta }}p(\theta |\mathbf{x})d\theta ={\int }_{\hat{\theta }}^{\infty }p(\theta |\mathbf{x})d\theta$$
  By definition $\hat{\theta }$ is the median of the posterior PDF or the point for which Pr ⁡ { θ ≤ θ ^ ∣ x } = 1 / 2 \Pr \{ \theta\le \hat \theta|\mathbf x\}=1/2 Pr{θ≤θ^∣x}=1/2.

• hit or miss cost function
  $$\begin{array}{rl}g(\hat{\theta })=& {\int }_{-\infty }^{\hat{\theta }-\delta }1\cdot p(\theta |\mathbf{x})d\theta +{\int }_{\hat{\theta }+\delta }^{\infty }1\cdot p(\theta |\mathbf{x})d\theta \\ =& 1-{\int }_{\hat{\theta }-\delta }^{\hat{\theta }+\delta }p(\theta |\mathbf{x})d\theta \end{array}$$
  This is minimized by maximizing
  $${\int }_{\dot{\theta }-\delta }^{\hat{\theta }+\delta }p(\theta |\mathbf{x})d\theta$$
  For $\delta$ arbitrarily small this is maximized by choosing $\hat{\theta }$ to correspond to the location of the maximum of $p(\theta |\mathbf{x})$. The estimator that minimizes the Bayes risk for the “hit-or-miss” cost function is therefore the mode (location of the maximum) of the posterior PDF. It is termed the maximum a posteriori (MAP) estimator and will be described in more detail later.

In summary, the estimators that minimize the Bayes risk for the cost functions of Figure 11.1 are the mean, median, and mode of the posterior PDF. This is illustrated in Figure 11.2a. For some posterior PDFs these three estimators are identical. A notable example is the Gaussian posterior PDF
$$p(\theta |\mathbf{x})=\frac{1}{\sqrt{2\pi {\sigma }_{\theta |x}^2}}\exp \left[-\frac{1}{2{\sigma }_{\theta |x}^2}(\theta -{\mu }_{\theta |x}{)}^2\right]$$
The mean ${\mu }_{\theta |x}$ is identical to the median (due to the symmetry) and the mode, as illustrated in Figure 11.2b.

Minimum Mean Square Error Estimators (MMSE)

In The Bayesian Philosophy the MMSE estimator was determined to be $E(\theta |\mathbf{x})$ or the mean of the posterior PDF. We continue our discussion of this important estimator by first extending it to the vector parameter case and then studying some of its properties.

If $\mathbf{\theta }$ is a vector parameter of dimension $p\times 1$, then to estimate ${\theta }_1$, for example, we may view the remaining parameters as nuisance parameters
$$p({\theta }_1|\mathbf{x})=\int \cdots \int p(\mathbf{\theta }|\mathbf{x})d{\theta }_2\cdots d{\theta }_p\text{\hspace{1em}(MMSE.1)}$$
where
$$p(\mathbf{\theta }|\mathbf{x})=\frac{p(\mathbf{x}|\mathbf{\theta })p(\mathbf{\theta })}{\int p(\mathbf{x}|\mathbf{\theta })p(\mathbf{\theta })d\mathbf{\theta }}\text{\hspace{1em}(MMSE.2)}$$
Then we have
$${\hat{\theta }}_i=E({\theta }_i|\mathbf{x})=\int {\theta }_{i}p({\theta }_i|\mathbf{x})d{\theta }_i,i=1,2,\cdots ,p\text{\hspace{1em}(MMSE.3)}$$
This is the MMSE estimator that minimizes
$$E[({\theta }_i-{\hat{\theta }}_i{)}^2]=\int \int ({\theta }_i-{\hat{\theta }}_i{)}^{2}p(\mathbf{x},{\theta }_i)d\mathbf{x}d{\theta }_i\text{\hspace{1em}(MMSE.4)}$$
or the squared error when averaged with respect to the marginal PDF $p(\mathbf{x},{\theta }_i)$. Alternatively, we can express the MMSE estimator for the first parameter from $(MMSE.1)$
$$\begin{array}{rl}{\hat{\theta }}_1& =\int {\theta }_{1}p({\theta }_1|\mathbf{x})d{\theta }_1\\ & =\int {\theta }_1\left[\int \cdots \int p(\mathbf{\theta }|\mathbf{x})d{\theta }_2\cdots d{\theta }_p\right]d{\theta }_1\\ & =\int {\theta }_{1}p(\mathbf{\theta }|\mathbf{x})d\mathbf{\theta }\end{array}$$
or in general
$${\hat{\theta }}_i=\int {\theta }_{i}p(\mathbf{\theta }|\mathbf{x})d\mathbf{\theta },i=1,2,\cdots ,p$$
In vector form we have
$$\hat{\mathbf{\theta }}=\left[\begin{array}{c}\int {\theta }_{1}p(\mathbf{\theta }|\mathbf{x})d\mathbf{\theta }\\ \int {\theta }_{2}p(\mathbf{\theta }|\mathbf{x})d\mathbf{\theta }\\ ⋮\\ \int {\theta }_{p}p(\mathbf{\theta }|\mathbf{x})d\mathbf{\theta }\end{array}\right]=\int \mathbf{\theta }p(\mathbf{\theta }|\mathbf{x})d\mathbf{\theta }=E(\mathbf{\theta }|\mathbf{x})\text{\hspace{1em}(MMSE.5)}$$
where the expectation is with respect to the posterior PDF of the vector parameter or $p(\mathbf{\theta }|\mathbf{x})$.

According to the definition, the minimum Bayesian MSE for ${\hat{\theta }}_i$
$$\mathrm{Bmse}({\hat{\theta }}_1)=E[({\theta }_1-{\hat{\theta }}_1{)}^2]=\int \int ({\theta }_1-{\hat{\theta }}_1{)}^{2}p(\mathbf{x},{\theta }_1)d{\theta }_{1}d\mathbf{x}$$
and since ${\hat{\theta }}_1=E({\theta }_1|\mathbf{x})$, we have
$$\begin{array}{rl}\mathrm{Bmse}({\hat{\theta }}_1)& =\int \int [{\theta }_1-E({\theta }_1|\mathbf{x}){]}^{2}p(\mathbf{x},{\theta }_1)d{\theta }_{1}d\mathbf{x}\\ & =\int \left[\int [{\theta }_1-E({\theta }_1|\mathbf{x}){]}^{2}p({\theta }_1|\mathbf{x})d{\theta }_1\right]p(\mathbf{x})d\mathbf{x}\end{array}$$
The posterior PDF can be written as
$$p({\theta }_1|\mathbf{x})=\int \cdots \int p(\mathbf{\theta }|\mathbf{x})d{\theta }_2\cdots d{\theta }_p$$
so that
$$\mathrm{Bmse}({\hat{\theta }}_1)=\int \left[\int [{\theta }_1-E({\theta }_1|\mathbf{x}){]}^{2}p(\mathbf{\theta }|\mathbf{x})d\mathbf{\theta }\right]p(\mathbf{x})d\mathbf{x}\text{\hspace{1em}(MMSE.6)}$$
which is just the $[1,1]$ element of ${\mathbf{C}}_{\theta |x}$, the covariance matrix of the posterior PDF. Hence, in general we have the minimum Bayesian MSE is
$$\mathrm{Bmse}({\hat{\theta }}_i)=\int [{\mathbf{C}}_{\theta |x}{]}_{ii}p(\mathbf{x})d\mathbf{x}i=1,2,\cdots ,p\text{\hspace{1em}(MMSE.7)}$$
where
$${\mathbf{C}}_{\theta |x}=E_{\theta |x}[(\mathbf{\theta }-E(\mathbf{\theta }|\mathbf{x}))(\mathbf{\theta }-E(\mathbf{\theta }|\mathbf{x}){)}^T]\text{\hspace{1em}(MMSE.8)}$$
We can obtain $\hat{\mathbf{\theta }}$ and ${\mathbf{C}}_{\theta |x}$ from Theorem 10.2 mentioned in The Bayesian Philosophy.

Theorem 10.2 (Conditional PDF of Multivariate Gaussian) If $\mathbf{x}$ and $\mathbf{y}$ are jointly Gaussian, where $\mathbf{x}$ is $k\times 1$ and $\mathbf{y}$ is $l\times 1,$ with mean vector ${\left[E(\mathbf{x}{)}^{T}E(\mathbf{y}{)}^T\right]}^T$ and
partitioned covariance matrix
$$\mathbf{C}=\left[\begin{array}{ll}{\mathbf{C}}_{xx}& {\mathbf{C}}_{xy}\\ {\mathbf{C}}_{yx}& {\mathbf{C}}_{yy}\end{array}\right]=\left[\begin{array}{ll}k\times k& k\times l\\ l\times k& l\times l\end{array}\right]\text{\hspace{1em}(BF.25)}$$
so that
$$p(\mathbf{x},\mathbf{y})=\frac{1}{(2\pi {)}^{\frac{k+1}{2}}{\det }^{\frac{1}{2}}(\mathbf{C})}\exp \left[-\frac{1}{2}{\left(\left[\begin{array}{l}\mathbf{x}-E(\mathbf{x})\\ \mathbf{y}-E(\mathbf{y})\end{array}\right]\right)}^T{\mathbf{C}}^{-1}\left(\left[\begin{array}{l}\mathbf{x}-E(\mathbf{x})\\ \mathbf{y}-E(\mathbf{y})\end{array}\right]\right)\right]$$
Then the conditional PDF $p(\mathbf{y}|\mathbf{x})$ is also Gaussian and
$$E(\mathbf{y}|\mathbf{x})=E(\mathbf{y})+{\mathbf{C}}_{yx}{\mathbf{C}}_{xx}^{-1}(\mathbf{x}-E(\mathbf{x}))\text{\hspace{1em}(BF.26)}$$

$${\mathbf{C}}_{y|x}={\mathbf{C}}_{yy}-{\mathbf{C}}_{yx}{\mathbf{C}}_{xx}^{-1}{\mathbf{C}}_{xy}\text{\hspace{1em}(BF.27)}$$

For general linear model, we can apply the result in Theorem 10.3.

Theorem 10.3 (Posterior PDF for the Bayesian General Linear Model) If the observed data $\mathbf{x}$ can be modeled as
$$\mathbf{x}=\mathbf{H}\mathbf{\theta }+\mathbf{w}$$
where $\mathbf{x}$ is an $N\times 1$ data vector, $\mathbf{H}$ is a known $N\times p$ matrix, $\mathbf{\theta }$ is a $p\times 1$ random vector with prior PDF $\mathcal{N}\left({\mu }_{\theta },{\mathbf{C}}_{\theta }\right),$ and $\mathbf{w}$ is an $N\times 1$ noise vector with PDF $\mathcal{N}\left(0,{\mathbf{C}}_w\right)$ and independent of $\mathbf{\theta },$ then the posterior PDF $p(\mathbf{\theta }|\mathbf{x})$ is Gaussian with mean
$$E(\mathbf{\theta }|\mathbf{x})={\mathbf{\mu }}_{\theta }+{\mathbf{C}}_{\theta }{\mathbf{H}}^T{\left(\mathbf{H}{\mathbf{C}}_{\theta }{\mathbf{H}}^T+{\mathbf{C}}_w\right)}^{-1}\left(\mathbf{x}-\mathbf{H}{\mathbf{\mu }}_{\theta }\right)\text{\hspace{1em}(BF.29)}$$
and covariance
$${\mathbf{C}}_{\theta \mid x}={\mathbf{C}}_{\theta }-{\mathbf{C}}_{\theta }{\mathbf{H}}^T{\left(\mathbf{H}{\mathbf{C}}_{\theta }{\mathbf{H}}^T+{\mathbf{C}}_w\right)}^{-1}\mathbf{H}{\mathbf{C}}_{\theta }\text{\hspace{1em}(BF.30)}$$
In contrast to the classical general linear model, $\mathbf{H}$ need not be full rank to ensure the invertibility of ${\mathbf{H}\mathbf{C}}_{\theta }{\mathbf{H}}^T+{\mathbf{C}}_w.$

Alternative formulation using Matrix inversion lemma:
$$E(\mathbf{\theta }|\mathbf{x})={\mathbf{\mu }}_{\theta }+{\left({\mathbf{C}}_{\theta }^{-1}+{\mathbf{H}}^T{\mathbf{C}}_w^{-1}\mathbf{H}\right)}^{-1}{\mathbf{H}}^T{\mathbf{C}}_w^{-1}\left(\mathbf{x}-\mathbf{H}{\mathbf{\mu }}_{\theta }\right)\text{\hspace{1em}(BF.31)}$$

$${\mathbf{C}}_{\theta \mid x}={\left({\mathbf{C}}_{\theta }^{-1}+{\mathbf{H}}^T{\mathbf{C}}_w^{-1}\mathbf{H}\right)}^{-1}\text{\hspace{1em}(BF.32)}$$

It is interesting to note that in the absence of prior knowledge in the Bayesian linear model the MMSE estimator yields the same form as the MVU estimator for the classical linear model. To verify this result note that from $(BF.31)$, for no prior knowledge ${\mathbf{C}}_{\theta }^{-1}\to 0$, therefore
$$\hat{\mathbf{\theta }}\to {\left({\mathbf{H}}^T{\mathbf{C}}_w^{-1}\mathbf{H}\right)}^{-1}{\mathbf{H}}^T{\mathbf{C}}_w^{-1}\mathbf{x}$$
which is recognized as the MVU estimator for the general linear model.

Properties

• The MMSE estimator commutes over linear transformations. Assume that we wish to estimate $\mathbf{\alpha }$ for
  $$\mathbf{\alpha }=\mathbf{A}\mathbf{\theta }+\mathbf{b}$$
  $\mathbf{\alpha }$ is a random vector for which the MMSE estimator is $\hat{\mathbf{\alpha }}=E(\mathbf{\alpha }|\mathbf{x})$. Because the linearity of the expectation operator
  $$\hat{\mathbf{\alpha }}=E(\mathbf{A}\mathbf{\theta }+\mathbf{b}|\mathbf{x})=\mathbf{A}E(\mathbf{\theta }|\mathbf{x})+\mathbf{b}=\mathbf{A}\hat{\mathbf{\theta }}+\mathbf{b}$$
  This holds regardless of the joint PDF $p(\mathbf{x},\mathbf{\theta })$.

• The MMSE has an additivity property for independent data sets. Assume ${\mathbf{x}}_1,{\mathbf{x}}_2$ are independent data vectors and $\mathbf{\theta },{\mathbf{x}}_1,{\mathbf{x}}_2$ are jointly Gaussian. Letting $\mathbf{x}=[{\mathbf{x}}_1^T{\mathbf{x}}_2^T{]}^T$, we have from Theorem 10.2
  $$\hat{\mathbf{\theta }}=E(\mathbf{\theta })+{\mathbf{C}}_{\theta x}{\mathbf{C}}_{xx}^{-1}(\mathbf{x}-E(\mathbf{x}))$$
  Since ${\mathbf{x}}_1,{\mathbf{x}}_2$ are independent,
  $${\mathbf{C}}_{xx}^{-1}={\left[\begin{array}{cc}{\mathbf{C}}_{x_1{x}_2}& 0\\ 0& {\mathbf{C}}_{x_2{x}_2}\end{array}\right]}^{-1}=\left[\begin{array}{cc}{\mathbf{C}}_{x_1{x}_2}^{-1}& 0\\ 0& {\mathbf{C}}_{x_2{x}_2}^{-1}\end{array}\right]$$

  $${\mathbf{C}}_{\theta x}=E\left[\mathbf{\theta }{\left[\begin{array}{c}{\mathbf{x}}_1\\ {\mathbf{x}}_2\end{array}\right]}^T\right]=[{\mathbf{C}}_{\theta x_1}{\mathbf{C}}_{\theta x_2}]$$

  Therefore
  $$\hat{\mathbf{\theta }}=E(\mathbf{\theta })+{\mathbf{C}}_{\theta x_1}{\mathbf{C}}_{x_1{x}_1}^{-1}({\mathbf{x}}_1-E({\mathbf{x}}_1))+{\mathbf{C}}_{\theta x_2}{\mathbf{C}}_{x_2{x}_2}^{-1}({\mathbf{x}}_2-E({\mathbf{x}}_2))$$

• In the jointly Gaussian case the MMSE is linear in the data, as can be seen from $(BF.26)$.

Example: Bayesian Fourier Analysis

Data model:
$$x[n]=a\cos 2\pi f_{0}n+b\sin 2\pi f_{0}n+w[n]n=0,1,\cdots ,N-1$$
where $f_0$ is a multiple of $1/N$, excepting $0$ or $1/2$ (for which $\sin 2\pi f_{0}n$ is identically zero), and $w[n]$ is WGN with variance ${\sigma }^2$. It is desired to estimate $\mathbf{\theta }=[ab{]}^T$.

To find the MMSE estimator, the data model is rewritten as
$$\mathbf{x}=\mathbf{H}\mathbf{\theta }+\mathbf{w}$$
where
$$\mathbf{H}=\left[\begin{array}{cc}1& 0\\ \cos 2\pi f_0& \sin 2\pi f_0\\ ⋮& ⋮\\ \cos \left[2\pi f_0(N-1)\right]& \sin \left[2\pi f_0(N-1)\right]\end{array}\right]$$
Because the columns of $\mathbf{H}$ are orthogonal, we have
$${\mathbf{H}}^T\mathbf{H}=\frac{N}{2}\mathbf{I}$$
Therefore it is more convenient to apply $(BF.31)$ and $(BF.32)$
$$\begin{array}{rl}\hat{\mathbf{\theta }}& ={\left(\frac{1}{{\sigma }_{\theta }^2}\mathbf{I}+{\mathbf{H}}^T\frac{1}{{\sigma }^2}\mathbf{H}\right)}^{-1}{\mathbf{H}}^T\frac{1}{{\sigma }^2}\mathbf{x}=\frac{1/{\sigma }^2}{1/{\sigma }_{\theta }^2+N/2{\sigma }^2}{\mathbf{H}}^T\mathbf{x}\\ {\mathbf{C}}_{\theta \mid x}& ={\left(\frac{1}{{\sigma }_{\theta }^2}\mathbf{I}+{\mathbf{H}}^T\frac{1}{{\sigma }^2}\mathbf{H}\right)}^{-1}=\frac{1}{1/{\sigma }_{\theta }^2+N/2{\sigma }^2}\mathbf{I}\end{array}$$
The MMSE estimator is
$$\begin{array}{rl}\hat{a}& =\frac{1}{1+\frac{2{\sigma }^2/N}{{\sigma }_{\theta }^2}}\left[\frac{2}{N}\sum _{n=0}^{N-1}x[n]\cos 2\pi f_{0}n\right]\\ \hat{b}& =\frac{1}{1+\frac{2{\sigma }^2/N}{{\sigma }_{\theta }^2}}\left[\frac{2}{N}\sum _{n=0}^{N-1}x[n]\sin 2\pi f_{0}n\right]\end{array}$$
with
$$\mathrm{B}\mathrm{m}\mathrm{s}\mathrm{e}(\hat{a})=\mathrm{B}\mathrm{m}\mathrm{s}\mathrm{e}(\hat{b})=\frac{1}{\frac{1}{{\sigma }_{\theta }^2}+\frac{1}{2{\sigma }^2/N}}$$

Maximum A Posteriori Estimators (MAP)

In the MAP estimation approach we choose $\hat{\theta }$ to maximize the posterior PDF or
$$\hat{\theta }=\arg \underset{\theta }{\max }p(\theta |\mathbf{x})$$
This was shown to minimize the Bayes risk for a “hit-or-miss” cost function $E[|{\theta }_i-{\hat{\theta }}_i|]$ . Since
$$p(\theta |\mathbf{x})=\frac{p(\mathbf{x}|\theta )p(\theta )}{p(\mathbf{x})}$$
It is equivalent to
$$\hat{\theta }=\arg \underset{\theta }{\max }p(\mathbf{x}|\theta )p(\theta )\text{\hspace{1em}(MAP.1)}$$
or
$$\hat{\theta }=\arg \underset{\theta }{\max }[\ln p(\mathbf{x}|\theta )+\ln p(\theta )]\text{\hspace{1em}(MAP.2)}$$
To extend the MAP estimator to the vector parameter case in which the posterior PDF is now $p(\mathbf{\theta }|\mathbf{x})$, we again employ the result
$$p({\theta }_1|\mathbf{x})=\int \cdots \int p(\mathbf{\theta }|\mathbf{x})d{\theta }_2\cdots d{\theta }_p\text{\hspace{1em}(MAP.3)}$$
Then the MAP estimator is given by
$${\hat{\theta }}_1=\arg \underset{{\theta }_1}{\max }p({\theta }_1|\mathbf{x})$$
or in general
$${\hat{\theta }}_i=\arg \underset{{\theta }_i}{\max }p({\theta }_i|\mathbf{x})i=1,2,\cdots ,p\text{\hspace{1em}(MAP.4)}$$
One of the advantages of the MAP estimator for a scalar parameter is that to numerically determine it we need only maximize $p(\mathbf{x}|\theta )p(\theta )$. No integration is required. This desirable property of the MAP estimator for a scalar parameter does not carry over to the vector parameter case due to the need to obtain $p({\theta }_i|\mathbf{x})$ as per $(MAP.3)$. However, we might propose the following vector MAP estimator
$$\hat{\mathbf{\theta }}=\arg \underset{\mathbf{\theta }}{\max }p(\mathbf{\theta }|\mathbf{x})=\arg \underset{\mathbf{\theta }}{\max }p(\mathbf{x}|\mathbf{\theta })p(\mathbf{\theta })\text{\hspace{1em}(MAP.5)}$$
in which the posterior PDF for the vector parameter $\mathbf{\theta }$ is maximized to find the estimator. The estimator in $(MAP.5)$ is not in general the same as $(MAP.4)$.

Properties

• If $p(\theta )$ is uniform and $p(\mathbf{x}|\theta )$ falls within this interval, then
  $$\hat{\theta }=\arg \underset{\theta }{\max }p(\mathbf{x}|\theta )p(\theta )=\arg \underset{\theta }{\max }p(\mathbf{x}|\theta )$$
  which is essentially the Bayesian MLE.

• Generally, it holds that if $N\to \infty$, the PDF $p(\mathbf{x}|\theta )$ becomes dominant over $p(\theta )$ and the MAP becomes thus identical to the Bayesian MLE.

• If $\mathbf{x}$ and $\mathbf{\theta }$ are jointly Gaussian, the posterior PDF is Gaussian and thus the MAP estimator is identical to the MMSE estimator.

• The invariance property encountered in maximum likelihood theory (MLE) does not hold for the MAP estimator.

Example: DC Level in WGN - Uniform Prior PDF

Recall the introductory example in The Bayesian Philosophy. There we discussed the MMSE estimator of $A$ for a DC level in WGN with a uniform prior PDF. The MMSE estimator could not be obtained in explicit form due to the need to evaluate the integrals. The posterior PDF was given as
$$p(A|\mathbf{x})=\{\begin{array}{rlrl}& \frac{\frac{1}{\sqrt{2\pi {\sigma }^2/N}}\exp \left[-\frac{1}{2{\sigma }^2/N}(A-\bar{x}{)}^2\right]}{{\int }_{-A_0}^{A_0}\frac{1}{\sqrt{2\pi {\sigma }^2/N}}\exp \left[-\frac{1}{2{\sigma }^2/N}(A-\bar{x}{)}^2\right]dA}& & |A|\le A_0\\ & 0,& & |A|>A_0\end{array}\text{\hspace{1em}(BP.9)}$$

The MAP estimator is explicitly given as
$$\hat{A}=\{\begin{array}{cccc}& -A_0& & \bar{x}<-A_0\\ & \bar{x}& & -A_0\le \bar{x}\le A_0\\ & A_0& & \bar{x}>A_0\end{array}$$
It can be seen that the MAP estimator is usually easier to determine since it does not involve any integration but only a maximization.

Linear MMSE Estimation (LMMSE)

We begin our discussion by assuming a scalar parameter $\theta$ is to be estimated based on the data set $\mathbf{x}=[x[0]x[1]\cdots x[N-1]{]}^T$. We do not assume any specific form for the joint PDF $p(\mathbf{x},\theta )$, but as we shall see shortly, only a knowledge of the first two moments.

We now consider the class of all linear estimators of the form
$$\hat{\theta }=\sum _{n=0}^{N-1}{a}_{n}x[n]+a_N\text{\hspace{1em}(LMMSE.1)}$$
and choosing the weighting coefficients $a_n$'s to minimize the Bayesian MSE
$$\mathrm{Bmse}(\hat{\theta })=E\left[(\theta -\hat{\theta }{)}^2\right]=\int \int (\theta -\hat{\theta }{)}^{2}p(\mathbf{x},\theta )d\theta d\mathbf{x}\text{\hspace{1em}(LMMSE.2)}$$
The resultant estimator is termed the linear minimum mean square error (LMMSE) estimator. Note that we included the $a_N$ coefficient to allow for nonzero means of $\mathbf{x}$ and $\mathbf{\theta }$.

Before determining the LMMSE estimator we should keep in mind that the estimator will be suboptimal unless the MMSE estimator happens to be linear. Such would be the case, for example, if the Bayesian linear model applied. Since the LMMSE estimator relies on the correlation between random variables, a parameter uncorrelated with the data cannot be linearly estimated.

We now derive the optimal weighting coefficients for use in $(LMMSE.1)$. Substituting $(LMMSE.1)$ into $(LMMSE.2)$ and differentiating w.r.t. $a_N$
$$\frac{\partial }{\partial a_N}E\left[{\left(\theta -\sum _{n=0}^{N-1}{a}_{n}x[n]-a_N\right)}^2\right]=-2E\left[\theta -\sum _{n=0}^{N-1}{a}_{n}x[n]-a_N\right]$$
Setting this equal to zero produces
$$a_N=E(\theta )-\sum _{n=0}^{N-1}{a}_{n}E(x[n])\text{\hspace{1em}(LMMSE.3)}$$
which is zero if the means are zero. Continuing, we need to minimize
$$\mathrm{B}\mathrm{m}\mathrm{s}\mathrm{e}(\hat{\theta })=E\left\{{\left[\sum _{n=0}^{N-1}{a}_n(x[n]-E(x[n]))-(\theta -E(\theta ))\right]}^2\right\}$$
over the remaining $a_n$'s, where $a_N$ has been replaced by $(LMMSE.3)$. Letting $\mathbf{a}={\left[a_0{a}_1\cdots a_{N-1}\right]}^T,$ we have
$$\begin{array}{rl}\mathrm{Bmse}(\hat{\theta })=& E\left\{{\left[{\mathbf{a}}^T(\mathbf{x}-E(\mathbf{x}))-(\theta -E(\theta ))\right]}^2\right\}\\ =& E\left[{\mathbf{a}}^T(\mathbf{x}-E(\mathbf{x}))(\mathbf{x}-E(\mathbf{x}){)}^T\mathbf{a}\right]-E\left[{\mathbf{a}}^T(\mathbf{x}-E(\mathbf{x}))(\theta -E(\theta ))\right]\\ & -E\left[(\theta -E(\theta ))(\mathbf{x}-E(\mathbf{x}){)}^T\mathbf{a}\right]+E\left[(\theta -E(\theta ){)}^2\right]\\ =& {\mathbf{a}}^T{\mathbf{C}}_{xx}\mathbf{a}-{\mathbf{a}}^T{\mathbf{C}}_{x\theta }-{\mathbf{C}}_{\theta x}\mathbf{a}+C_{\theta \theta }\end{array}\text{\hspace{1em}(LMMSE.4)}$$
where ${\mathbf{C}}_{xx}$ is the $N\times N$ covariance matrix of $\mathbf{x},$ and ${\mathbf{C}}_{\theta x}$ is the $1\times N$ cross-covariance vector having the property that ${\mathbf{C}}_{\theta x}^T={\mathbf{C}}_{x\theta },$ and $C_{\theta \theta }$ is the variance of $\theta .$ Taking the gradient to yield
$$\frac{\partial \mathrm{B}\mathrm{m}\mathrm{s}\mathrm{e}(\hat{\theta })}{\partial \mathbf{a}}=2{\mathbf{C}}_{xx}\mathbf{a}-2{\mathbf{C}}_{x\theta }$$
which when set to zero results in
$$\mathbf{a}={\mathbf{C}}_{xx}^{-1}{\mathbf{C}}_{x\theta }\text{\hspace{1em}(LMMSE.5)}$$
Using $(LMMSE.3)$ and $(LMMSE.5)$ in $(LMMSE.1)$ produces
$$\begin{array}{rl}\hat{\theta }& ={\mathbf{a}}^T\mathbf{x}+a_N\\ & ={\mathbf{C}}_{x\theta }^T{\mathbf{C}}_{xx}^{-1}\mathbf{x}+E(\theta )-{\mathbf{C}}_{x\theta }^T{\mathbf{C}}_{xx}^{-1}E(\mathbf{x})\end{array}$$
or finally the LMMSE estimator is
$$\hat{\theta }=E(\theta )+{\mathbf{C}}_{\theta x}{\mathbf{C}}_{xx}^{-1}(\mathbf{x}-E(\mathbf{x}))\text{\hspace{1em}(LMMSE.6)}$$
Note that it is identical in form to the MMSE estimator for jointly Gaussian $\mathbf{x}$ and $\theta$, as can be verified from $(BF.26)$. This is because in the Gaussian case the MMSE estimator happens to be linear, and hence our constraint is automatically satisfied.

The minimum Bayesian MSE is obtained by substituting $(LMMSE.5)$ into $(LMMSE.4)$ to yield
$$\begin{array}{rl}\mathrm{B}\mathrm{m}\mathrm{s}\mathrm{e}(\hat{\theta })=& {\mathbf{C}}_{x\theta }^T{\mathbf{C}}_{xx}^{-1}{\mathbf{C}}_{xx}{\mathbf{C}}_{xx}^{-1}{\mathbf{C}}_{x\theta }-{\mathbf{C}}_{x\theta }^T{\mathbf{C}}_{xx}^{-1}{\mathbf{C}}_{x\theta }-{\mathbf{C}}_{\theta x}{\mathbf{C}}_{xx}^{-1}{\mathbf{C}}_{x\theta }+C_{\theta \theta }\\ =& {\mathbf{C}}_{\theta x}{\mathbf{C}}_{xx}^{-1}{\mathbf{C}}_{x\theta }-2{\mathbf{C}}_{\theta x}{\mathbf{C}}_{xx}^{-1}{\mathbf{C}}_{x\theta }+C_{\theta \theta }\\ =& C_{\theta \theta }-{\mathbf{C}}_{\theta x}{\mathbf{C}}_{xx}^{-1}{\mathbf{C}}_{x\theta }\end{array}\text{\hspace{1em}(LMMSE.7)}$$
In general, all that is required to determine the LMMSE estimator are the first two moments of $p(\mathbf{x},\theta )$ or
$$\left[\begin{array}{c}E(\theta )\\ E(\mathbf{x})\end{array}\right],\left[\begin{array}{cc}{C}_{\theta \theta }& {\mathbf{C}}_{\theta x}\\ {\mathbf{C}}_{x\theta }& {\mathbf{C}}_{xx}\end{array}\right]$$

Geometrical Interpretations

In this formulation, we assume that $\theta$ and $\mathbf{x}$ are zero mean. If not, we can always define the zero mean random variables ${\theta }^{\prime }=\theta -E(\theta )$ and ${\mathbf{x}}^{\prime }=\mathbf{x}-E(\mathbf{x})$, and consider estimation of ${\theta }^{\prime }$ by a linear function of ${\mathbf{x}}^{\prime }$. Now we wish to find the $a_n$'s so that
$$\hat{\theta }=\sum _{n=0}^{N-1}{a}_{n}x[n]$$
minimizes
$$\mathrm{Bmse}(\hat{\theta })=E\left[(\theta -\hat{\theta }{)}^2\right]=\int \int (\theta -\hat{\theta }{)}^{2}p(\mathbf{x},\theta )d\theta d\mathbf{x}$$
It can be shown that an appropriate definition, i.e., one that satisfies the properties of an inner product between the vectors $x$ and $y$, is
$$(x,y)=E(xy)\text{\hspace{1em}(LMMSE.8)}$$
With this definition we have that
$$(x,x)=E(x^2)=\Vert x{\Vert }^2\text{\hspace{1em}(LMMSE.9)}$$
Also, we can now define two vectors to be orthogonal if
$$(x,y)=E(xy)=0\text{\hspace{1em}(LMMSE.10)}$$
Since the vectors are zero mean, this is equivalent to saying that two vectors are orthogonal if and only if they are uncorrelated.

Then the Bayesian MSE can be presented as
$$E\left[(\theta -\hat{\theta }{)}^2\right]={\Vert \theta -\sum _{n=0}^{N-1}{a}_{n}x[n]\Vert }^2$$
This means that minimization of the MSE is equivalent to a minimization of the squared length of the error vector $ϵ=\theta -\hat{\theta }$. Clearly, the length of the error vector is minimized when $ϵ$ is orthogonal to the subspace spanned by { x [ 0 ] , x [ 1 ] , . . . , x [ N − 1 ] } \{x[0], x[1], ... , x[N - 1]\} {x[0],x[1],...,x[N−1]}. Hence, we require
$$ϵ\perp x[0],x[1],\cdots ,x[N-1]\text{\hspace{1em}(LMMSE.11)}$$
or by using our definition of orthogonality
$$E\left[(\theta -\hat{\theta })x[n]\right]=0n=0,1,\cdots ,N-1\text{\hspace{1em}(LMMSE.12)}$$
This is the important orthogonality principle. It says that in estimating the realization of a random variable by a linear combination of data samples, the optimal estimator is obtained when the error is orthogonal to each data sample. Using the orthogonality principle, the weighting coefficients are easily found as
$$\sum _{n=0}^{N-1}{a}_{m}E(x[m]x[n])=E(\theta x[n])n=0,1,\cdots ,N-1$$
In matrix form this is
$$\left[\begin{array}{cccc}E\left(x^2[0]\right)& E(x[0]x[1])& \dots & E(x[0]x[N-1])\\ E(x[1]x[0])& E\left(x^2[1]\right)& \dots & E(x[1]x[N-1])\\ ⋮& ⋮& \ddots & ⋮\\ E(x[N-1]x[0])& E(x[N-1]x[1])& \dots & E\left(x^2[N-1]\right)\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{N-1}\end{array}\right]=\left[\begin{array}{c}E(\theta x[0])\\ E(\theta x[1])\\ ⋮\\ E(\theta x[N-1])\end{array}\right]$$
These are normal equations. The matrix is recognized as ${\mathbf{C}}_{xx}$, and the right-hand vector as ${\mathbf{C}}_{x\theta }$. Therefore,
$${\mathbf{C}}_{xx}\mathbf{a}={\mathbf{C}}_{x\theta }$$
and
$$\mathbf{a}={\mathbf{C}}_{xx}^{-1}{\mathbf{C}}_{x\theta }$$
The LMMSE estimator of $\theta$ is
$$\hat{\theta }={\mathbf{a}}^T\mathbf{x}={\mathbf{C}}_{\theta x}{\mathbf{C}}_{xx}^{-1}\mathbf{x}$$
in agreement with $(LMMSE.5)$ and $(LMMSE.6)$. The minimum Bayesian MSE is the squared length of the error vector or
$$\begin{array}{rl}E\left[(\theta -\hat{\theta }{)}^2\right]& =E\left[(\theta -\sum _{n=0}^{N-1}{a}_{n}x[n])ϵ\right]\\ & =E\left[\theta ϵ\right]\\ & =E({\theta }^2)-\sum _{n=0}^{N-1}{a}_{n}E(x[n]\theta )\\ & =C_{\theta \theta }-{\mathbf{a}}^T{\mathbf{C}}_{x\theta }\\ & =C_{\theta \theta }-{\mathbf{C}}_{\theta x}{\mathbf{C}}_{xx}^{-1}{\mathbf{C}}_{x\theta }\end{array}$$
in agreement with $(LMMSE.7)$.

Example: DC Level in WGN - Uniform Prior PDF

The data model is
$$x[n]=A+w[n]n=0,1,\dots ,N-1$$
where $A\sim \mathcal{U}\left[-A_0,A_0\right],w[n]$ is WGN with variance ${\sigma }^2,$ and $A$ and $w[n]$ are independent. We wish to estimate $A$. The MMSE estimator cannot be obtained in closed form due to the integration required. Applying the LMMSE estimator, we first note that $E(A)=0,$ and hence $E(x[n])=0.$ since $E(\mathbf{x})=0,$ the covariances are
$$\begin{array}{rl}{\mathbf{C}}_{xx}& =E\left(\mathbf{x}{\mathbf{x}}^T\right)\\ & =E\left[(A1+\mathbf{w})(A1+\mathbf{w}{)}^T\right]\\ & =E\left(A^2\right)11^T+{\sigma }^2\mathbf{I}\\ {\mathbf{C}}_{\theta x}& =E\left(A{\mathbf{x}}^T\right)\\ & =E\left[A(A1+\mathbf{w}{)}^T\right]\\ & =E\left(A^2\right)1^T\end{array}$$
Hence, from $(LMMSE.6)$
$$\begin{array}{rl}\hat{A}& ={\mathbf{C}}_{\theta x}{\mathbf{C}}_{xx}^{-1}\mathbf{x}\\ & ={\sigma }_A^2{1}^T{\left({\sigma }_A^{2}11^T+{\sigma }^2\mathbf{I}\right)}^{-1}\mathbf{x}\\ & =\frac{\frac{A_0^2}{3}}{\frac{A_0^2}{3}+\frac{{\sigma }^2}{N}}\bar{x}\end{array}$$
where we have let ${\sigma }_A^2=E\left(A^2\right).$

The Vector LMMSE Estimator

Now we wish to find the linear estimator that minimizes the Bayesian MSE for each element. We assume that
$${\hat{\theta }}_i=\sum _{n=0}^{N-1}{a}_{in}x[n]+a_{iN}\text{\hspace{1em}(LMMSE.13)}$$
for $i=1,2,\dots ,p$ and choose the weighting coefficients to minimize
$$\mathrm{Bmse}\left({\hat{\theta }}_i\right)=E\left[{\left({\theta }_i-{\hat{\theta }}_i\right)}^2\right]i=1,2,\dots ,p$$
where the expectation is with respect to $p\left(\mathbf{x},{\theta }_i\right)$. Since we are actually determining $p$ separate estimators, the scalar solution can be applied, and we obtain from $(LMMSE.6)$
$${\hat{\theta }}_i=E\left({\theta }_i\right)+{\mathbf{C}}_{{\theta }_{i}x}{\mathbf{C}}_{xx}^{-1}(\mathbf{x}-E(\mathbf{x}))i=1,2,\dots ,p$$
and the minimum Bayesian MSE is from $(LMMSE.7)$
$$\mathrm{B}\mathrm{m}\mathrm{s}\mathrm{e}\left({\hat{\theta }}_i\right)=C_{{\theta }_i{\theta }_i}-{\mathbf{C}}_{{\theta }_i}{\mathbf{C}}_{xx}^{-1}{\mathbf{C}}_{x{\theta }_i}i=1,2,\dots ,p$$
The scalar LMMSE estimators can be combined into a vector estimator as
$$\begin{array}{rl}\hat{\mathbf{\theta }}& =\left[\begin{array}{c}E\left({\theta }_1\right)\\ E\left({\theta }_2\right)\\ ⋮\\ E\left({\theta }_p\right)\end{array}\right]+\left[\begin{array}{c}{\mathbf{C}}_{{\theta }_{1}x}{\mathbf{C}}_{xx}^{-1}(\mathbf{x}-E(\mathbf{x}))\\ {\mathbf{C}}_{{\theta }_{2}x}{\mathbf{C}}_{xx}^{-1}(\mathbf{x}-E(\mathbf{x}))\\ ⋮\\ {\mathbf{C}}_{{\theta }_{p}x}{\mathbf{C}}_{xx}^{-1}(\mathbf{x}-E(\mathbf{x}))\end{array}\right]\\ & =\left[\begin{array}{c}E\left({\theta }_1\right)\\ E\left({\theta }_2\right)\\ ⋮\\ E\left({\theta }_p\right)\end{array}\right]+\left[\begin{array}{c}{\mathbf{C}}_{{\theta }_{1}x}\\ {\mathbf{C}}_{{\theta }_{2}x}\\ ⋮\\ {\mathbf{C}}_{{\theta }_{p}x}\end{array}\right]{\mathbf{C}}_{xx}^{-1}(\mathbf{x}-E(\mathbf{x}))\\ & =E(\mathbf{\theta })+{\mathbf{C}}_{\theta x}{\mathbf{C}}_{xx}^{-1}(\mathbf{x}-E(\mathbf{x}))\end{array}\text{\hspace{1em}(LMMSE.14)}$$
where now ${\mathbf{C}}_{\theta x}$ is a $p\times N$ matrix.

By a similar approach we find that the Bayesian MSE matrix is (see Problem 12.7)
$$\begin{array}{rl}{\mathbf{M}}_{\hat{\theta }}& =E\left[(\mathbf{\theta }-\hat{\mathbf{\theta }})(\mathbf{\theta }-\hat{\mathbf{\theta }}{)}^T\right]\\ & ={\mathbf{C}}_{\theta \theta }-{\mathbf{C}}_{\theta x}{\mathbf{C}}_{xx}^{-1}{\mathbf{C}}_{x\theta }\end{array}\text{\hspace{1em}(LMMSE.15)}$$
where ${\mathbf{C}}_{\theta \theta }$ is the $p\times p$ covariance matrix. Consequently, the minimum Bayesian $\mathrm{M}\mathrm{S}\mathrm{E}$ is
$$\mathrm{Bmse}({\hat{\theta }}_i)=[{\mathbf{M}}_{\hat{\theta }}{]}_{ii}\text{\hspace{1em}(LMMSE.16)}$$
Of course, these results are identical to those for the MMSE estimator in the Gaussian case, for which the estimator is linear.

Properties

• LMMSE estimators are identical in form to the MMSE estimator for jointly Gaussian $\mathbf{x}$ and $\mathbf{\theta }$

• LMMSE estimator commutes over linear transformations. This is to say that if
  $$\mathbf{\alpha }=\mathbf{A}\mathbf{\theta }+\mathbf{b}$$
  then the LMMSE estimator of $\alpha$ is

$$\hat{\mathbf{\alpha }}=\mathbf{A}\hat{\mathbf{\theta }}+\mathbf{b}$$
​ with $\hat{\theta }$ given by $(LMMSE.14)$.

• The LMMSE estimator of a sum of unknown parameters is the sum of the individual estimators. Specifically, if we wish to estimate $\alpha ={\theta }_1+{\theta }_2,$ then

$$\hat{\mathbf{\alpha }}={\hat{\mathbf{\theta }}}_1+{\hat{\mathbf{\theta }}}_2$$

• In analogy with the BLUE there is a corresponding Gauss-Markov theorem for the Bayesian case

  Theorem 12.1 (Bayesian Gauss-Markov Theorem) If the data are described by the Bayesian linear model form
  $$\mathbf{x}=\mathbf{H}\mathbf{\theta }+\mathbf{w}$$
  where $\mathbf{x}$ is an $N\times 1$ data vector, $\mathbf{H}$ is a known $N\times p$ observation matrix, $\mathbf{\theta }$ is a $p\times 1$ random vector of parameters whose realization is to be estimated and has mean $E(\mathbf{\theta })$ and covariance matrix ${\mathrm{C}}_{\theta \theta }$, and $\mathbf{w}$ is an $N\times 1$ random vector with zero mean and covariance matrix ${\mathrm{C}}_w$ and is uncorrelated with $\mathbf{\theta }$, then the LMMSE estimator of $\mathbf{\theta }$ is
  $$\begin{array}{rl}\hat{\mathbf{\theta }}& =E(\mathbf{\theta })+{\mathbf{C}}_{\theta \theta }{\mathbf{H}}^T{\left(\mathbf{H}{\mathbf{C}}_{\theta \theta }{\mathbf{H}}^T+{\mathbf{C}}_w\right)}^{-1}(\mathbf{x}-\mathbf{H}E(\mathbf{\theta }))\\ & =E(\mathbf{\theta })+{\left({\mathbf{C}}_{\theta \theta }^{-1}+{\mathbf{H}}^T{\mathbf{C}}_w^{-1}\mathbf{H}\right)}^{-1}{\mathbf{H}}^T{\mathbf{C}}_w^{-1}(\mathbf{x}-\mathbf{H}E(\mathbf{\theta }))\end{array}$$
  The performance of the estimator is measured by the error $\mathbf{ϵ}=\mathbf{\theta }-\hat{\mathbf{\theta }}$ whose mean is zero and whose covariance matrix is
  $$\begin{array}{rl}{\mathbf{C}}_{ϵ}& =E_{x,\theta }\left(\mathbf{ϵ}{\mathbf{ϵ}}^T\right)\\ & ={\mathbf{C}}_{\theta \theta }-{\mathbf{C}}_{\theta \theta }{\mathbf{H}}^T{\left(\mathbf{H}{\mathbf{C}}_{\theta \theta }{\mathbf{H}}^T+{\mathbf{C}}_w\right)}^{-1}\mathbf{H}{\mathbf{C}}_{\theta \theta }\\ & ={\left({\mathbf{C}}_{\theta \theta }^{-1}+{\mathbf{H}}^T{\mathbf{C}}_w^{-1}\mathbf{H}\right)}^{-1}\end{array}$$
  The error covariance matrix is also the minimum MSE matrix ${\mathbf{M}}_{\hat{\theta }}$ whose diagonal elements yield the minimum Bayesian MSE
  $$\begin{array}{rl}{\left[{\mathbf{M}}_{\hat{\theta }}\right]}_{ii}& ={\left[{\mathbf{C}}_{ϵ}\right]}_{ii}=\mathrm{Bmse}\left({\hat{\theta }}_i\right)\end{array}$$

These results are identical to those in Theorem 10.3 (Posterior PDF for the Bayesian General Linear Model) except that the error vector is not necessarily Gaussian.