Optimal estimation

In applied statistics, optimal estimation is a regularized matrix inverse method based on Bayes' theorem. It is used very commonly in the geosciences, particularly for atmospheric sounding. A matrix inverse problem looks like this:

A x → = y → {\displaystyle \mathbf {A} {\vec {x}}={\vec {y}}}

The essential concept is to transform the matrix, A, into a conditional probability and the variables, x → {\displaystyle {\vec {x}}} and y → {\displaystyle {\vec {y}}} into probability distributions by assuming Gaussian statistics and using empirically determined covariance matrices.

Typically, one expects the statistics of most measurements to be Gaussian. So for example for P ( y → | x → ) {\displaystyle P({\vec {y}}|{\vec {x}})} , we can write:

P ( y → | x → ) = 1 ( 2 π ) m n / 2 | S y | exp ⁡ [ − 1 2 ( A x → − y → ) T S y − 1 ( A x → − y → ) ] {\displaystyle P({\vec {y}}|{\vec {x}})={\frac {1}{(2\pi )^{mn/2}|{\boldsymbol {S_{y}}}|}}\exp \left[-{\frac {1}{2}}({\boldsymbol {A}}{\vec {x}}-{\vec {y}})^{T}{\boldsymbol {S_{y}}}^{-1}({\boldsymbol {A}}{\vec {x}}-{\vec {y}})\right]}

where m and n are the numbers of elements in x → {\displaystyle {\vec {x}}} and y → {\displaystyle {\vec {y}}} respectively A {\displaystyle {\boldsymbol {A}}} is the matrix to be solved (the linear or linearised forward model) and S y {\displaystyle {\boldsymbol {S_{y}}}} is the covariance matrix of the vector y → {\displaystyle {\vec {y}}} . This can be similarly done for x → {\displaystyle {\vec {x}}} :

P ( x → ) = 1 ( 2 π ) m / 2 | S x a | exp ⁡ [ − 1 2 ( x → − x a ^ ) T S x a − 1 ( x → − x a ^ ) ] {\displaystyle P({\vec {x}})={\frac {1}{(2\pi )^{m/2}|{\boldsymbol {S_{x_{a}}}}|}}\exp \left[-{\frac {1}{2}}({\vec {x}}-{\widehat {x_{a}}})^{T}{\boldsymbol {S_{x_{a}}}}^{-1}({\vec {x}}-{\widehat {x_{a}}})\right]}

Here P ( x → ) {\displaystyle P({\vec {x}})} is taken to be the so-called "a-priori" distribution: x a ^ {\displaystyle {\widehat {x_{a}}}} denotes the a-priori values for x → {\displaystyle {\vec {x}}} while S x a {\displaystyle {\boldsymbol {S_{x_{a}}}}} is its covariance matrix.

The nice thing about the Gaussian distributions is that only two parameters are needed to describe them and so the whole problem can be converted once again to matrices. Assuming that P ( x → | y → ) {\displaystyle P({\vec {x}}|{\vec {y}})} takes the following form:

P ( x → | y → ) = 1 ( 2 π ) m n / 2 | S x | exp ⁡ [ − 1 2 ( x → − x ^ ) T S x − 1 ( x → − x ^ ) ] {\displaystyle P({\vec {x}}|{\vec {y}})={\frac {1}{(2\pi )^{mn/2}|{\boldsymbol {S_{x}}}|}}\exp \left[-{\frac {1}{2}}({\vec {x}}-{\widehat {x}})^{T}{\boldsymbol {S_{x}}}^{-1}({\vec {x}}-{\widehat {x}})\right]}

P ( y → ) {\displaystyle P({\vec {y}})} may be neglected since, for a given value of x → {\displaystyle {\vec {x}}} , it is simply a constant scaling term. Now it is possible to solve for both the expectation value of x → {\displaystyle {\vec {x}}} , x ^ {\displaystyle {\widehat {x}}} , and for its covariance matrix by equating P ( x → | y → ) {\displaystyle P({\vec {x}}|{\vec {y}})} and P ( y → | x → ) P ( x → ) {\displaystyle P({\vec {y}}|{\vec {x}})P({\vec {x}})} . This produces the following equations:

S x = ( A T S y − 1 A + S x a − 1 ) − 1 {\displaystyle {\boldsymbol {S_{x}}}=({\boldsymbol {A}}^{T}{\boldsymbol {S_{y}^{-1}}}{\boldsymbol {A}}+{\boldsymbol {S_{x_{a}}^{-1}}})^{-1}}

x ^ = x a ^ + S x A T S y − 1 ( y → − A x a ^ ) {\displaystyle {\widehat {x}}={\widehat {x_{a}}}+{\boldsymbol {S_{x}}}{\boldsymbol {A}}^{T}{\boldsymbol {S_{y}}}^{-1}({\vec {y}}-{\boldsymbol {A}}{\widehat {x_{a}}})}

Because we are using Gaussians, the expected value is equivalent to the maximum likely value, and so this is also a form of maximum likelihood estimation.

Typically with optimal estimation, in addition to the vector of retrieved quantities, one extra matrix is returned along with the covariance matrix. This is sometimes called the resolution matrix or the averaging kernel and is calculated as follows:

R = ( A T S y − 1 A + S x a − 1 ) − 1 A T S y − 1 A {\displaystyle {\boldsymbol {R}}=({\boldsymbol {A}}^{T}{\boldsymbol {S_{y}}}^{-1}{\boldsymbol {A}}+{\boldsymbol {S_{x_{a}}}}^{-1})^{-1}{\boldsymbol {A}}^{T}{\boldsymbol {S_{y}}}^{-1}{\boldsymbol {A}}}

This tells us, for a given element of the retrieved vector, how much of the other elements of the vector are mixed in. In the case of a retrieval of profile information, it typical indicates the altitude resolution for a given altitude. For instance if the resolution vectors for all the altitudes contain non-zero elements (to a numerical tolerance) in their four nearest neighbours, then the altitude resolution is only one fourth that of the actual grid size.

• Clive D. Rodgers (1976). "Retrieval of Atmospheric Temperature and Composition From Remote Measurements of Thermal Radiation". Reviews of Geophysics and Space Physics. 14 (4): 609. Bibcode:1976RvGSP..14..609R. doi:10.1029/RG014i004p00609.

• Clive D. Rodgers (2000). Inverse Methods for Atmospheric Sounding: Theory and Practice. World Scientific.

• Clive D. Rodgers (2002). "Atmospheric Remote Sensing: The Inverse Problem". Proceedings of the Fourth Oxford/RAL Spring School in Quantitative Earth Observation. University of Oxford.