3.1.

Metropolis Algorithm Search for MAP Solution

Suppose that we have in our possession prior knowledge that a target patch, forming part of the lower boundary of the atmosphere, is composed of an intimate mixture of $m+1$ spectral endmembers $\{{\epsilon }_i(k)\}$ at temperature $T$. The label $k$ may, depending upon context, refer to the wavenumber or to a finite number of wavenumber-averaged spectral bands.

For the development that follows, it is desirable to restate some of the main definitions used in spectral mixing theory using language borrowed from topology.^21,22 Let $m+1$ distinct points ${\mathbf{y}}_0,{\mathbf{y}}_1,\cdots {\mathbf{y}}_m$ in ${\mathbf{R}}^{\mathbf{m}}$ be chosen so that the vectors ${\mathbf{y}}_1-{\mathbf{y}}_0,{\mathbf{y}}_2-{\mathbf{y}}_0\cdots {\mathbf{y}}_m-{\mathbf{y}}_0$ are linearly independent. Then the set

with and is an $m$-simplex.

A spectral mixture amounts to a mapping into a geometric $m$-simplex, whose vertices have spectral endmembers at a stipulated temperature $T$ for “coefficients.” A mixture with stipulated weights ${\lambda }_i$ corresponds to the vector

The “interior” of $K_m$ is the subset of $K_m$ for which ${\lambda }_i>0$, i.e., the closure of its interior. The “polyhedron” of $K_m$, denoted as $|K_m|$, is the set comprised of the points of $\mathbf{x}\in K_m$ considered as a subset of ${\mathbf{R}}^m$ and is a convex compact subset of ${\mathbf{R}}^{\mathbf{m}}$.

In the case $m+1=3$, a familiar example of a 2-simplex is the ternary diagram used to classify phreatic igneous rocks. The double three-component diagram used in the quartz, alkali feldspar, plagioclase, feldspathoid classification²³ scheme is a union of two 2-simplices, and is an example of a “simplicial complex.”

For present purposes, the $i$’th pure endmember for the $n$’th trial is assigned to the $i$’th vertex of $K_m$

with the spectral mixture corresponding to a point in the polyhedron of $K_m$.

It is necessary to account for surface temperature in a somewhat different way. Let the minimum and maximum physically admissible surface temperatures be $T_{\min }$ and $T_{\max }$, respectively. Then the temperature of our target patch is given as

with Corresponding to $\mathbf{x}$ introduced already, we have from Eq. (6) the unit interval, with The quantity that appears in the forward model for the $n$’th trial is where $B_k(T_n)$ is the (band-integrated, as needed) Planck function at temperature $T_n$. The parametrization of the choice $\{T_n,\epsilon (k)\}$ in terms of the vector $\mathbf{x}$ is a mapping into the topological product of $I^1$ and $|K_m|$. The set $H^{m+1}$ is not a simplex nor is it necessarily a simplicial complex. It is, however, a convex polytope and is the convex hull of its vertices ${\mathbf{x}}_i$, $0\le i\le m+1$. Although we will not need it in what follows, $H^{m+1}$ can be decomposed into either a simplicial complex or a union of simplices.

We score trial mixtures by that we most wish to maximize; the posterior probability for the observed spectral radiance to originate from a surface patch with temperature $T$ and spectral emissivity $\epsilon (k)$. A standard argument^3,5 gives the posterior probability in terms of a MAXENT estimator

in terms of a forward model that is some function of the $n$’th trial, in each spectral bin $k$. We note that while the equation of transfer is linear, the dependence of its solution $I_{\mathrm{FM}}$ upon ${\epsilon }_i(k)B_k(T_n)$ need not be. The assumed noise variance ${\sigma }^2$ is shown as having a formal dependence upon a parameter, the “annealing temperature” $T_{\mathrm{a}}$, which governs the annealing schedule for the search for an MAP solution. The joint posterior probability in $J$ spectral bands is proportional to If radiance $I_k$ in each of $J$ bands originating from a patch on the Earth’s surface has been detected at the top of the atmosphere, the posterior probability that the surface patch is at a temperature $T$ given prior knowledge $K$ is given by Bayes’ theorem as The noise variance is assumed to be known and the functional dependence of probabilities upon ${\sigma }_i$ is omitted. The prior probability $P(\{I_i\}|K)$ for the radiances $\{I_k\}$ has no dependence upon $T$ and for our purposes may be absorbed into an overall normalization.²⁴ Equation (15) is evaluated with the aid of the prior probability for the surface to be at temperature $T$ and has spectral emissivity $\epsilon (k)$, given available knowledge $K$⁵ where $P[T,{\epsilon }_i(k)|\{I_i\},K]$ is the conditional probability for the hypothesis that the surface temperature is $T$, and the spectral emissivity ${\epsilon }_k$, given observed radiances $\{I_i\}$ and prior knowledge $K$.

Each trial is thus scored according to the joint posterior probability for observed spectral radiance $I_i$ to result from surface temperature $T$ and spectral emissivity ${\epsilon }_k$

where $x$ stands for $\{\mathbf{x},x_{m+1}\}$. Thus, in going from the $(n-1)$’th to the $n$’th trial, the $n$’th candidate mixture is selected by Monte Carlo draw and $p_n$ for the new trial is compared to $p_{n-1}$ for the last one. The probability that it is accepted is^18,20 where the probability $P(T_{\mathrm{a}})$ of taking a downward step in $p_n$ is determined by the annealing schedule. The dependence of $P(T_{\mathrm{a}})$ on the annealing schedule is symbolized by the annealing temperature $T_{\mathrm{a}}$, which is taken to decrease systematically during the MAP search. The actual form $P(T_{\mathrm{a}})$ takes in practical calculations and is determined empirically.

3.2.

Convergence

We now examine the question of convergence. Corresponding to the sequence of $m$-simplices, $K_m$, as the number of trials $n$ increases without bound is a sequence of trials $\{T_n,{\epsilon }_i(k)\}$ with associated loci $\{x\}\in H^{m+1}$.

As a closed bounded subset of ${\mathbf{R}}^{m+1}$, $H^{m+1}$ is a compactum. Therefore, as $n\to \infty$, the sequence of trials $x$ contains a convergent subsequence, whatever the value of $m$. Correspondingly, the sequence of posterior probabilities likewise has a convergent subsequence that, by construction, tends to the maximum value of the posterior probability, i.e., to an MAP solution for $T$ and $\epsilon (k)$.

Consider the map $x^{\prime }=\mathrm{\Phi }(x)$ given by

The mapping of Eq. (19) gives the action of the Metropolis algorithm according to Eq. (18) at sufficiently late times in the annealing schedule so that a transition to a state of decreased posterior probability occurs rarely, in the limit, almost never. We have noted that at a sufficiently late point in the annealing schedule, trials that decrease the posterior probability [Eq. (17)] will become infrequent. We may elide any such trials without affecting the convergence of the subsequence, which then takes the form For all $n$ greater than some $M$, convergence of the subsequence implies the Cauchy condition (with the Euclidean norm supplying a suitable metric for finite $m$) so that

The mapping of Eq. (19) generates a sequence of trials $x$ for which $p_n$ is nondecreasing. By Zorn’s lemma, the set comprised of all admissible trials $x$ has at least one element with a maximal value of $p_n$. We note that maximizing $p_n$ also maximizes the information-theoretic entropy by Eq. (12). According to the usual statement of the second law, the state of maximum entropy is one of thermodynamic equilibrium. We may, in view of the thermodynamic analogy underlying the Metropolis algorithm, therefore call the limit Eq. (22) an equilibrium point.

This nomenclature is attractive for another reason. In the limit, Eq. (22) amounts to a fixed point of Eq. (19). Ordinary fixed-point theorems are inapplicable to Eq. (19) because it is neither continuous nor semicontinuous: It can map an open set $\in H^{m+1}$ to a singleton $x^{\prime }$. We can, however, adapt the celebrated construction introduced by Nash²⁵ to prove the existence of a fixed (equilibrium) point of an equivalent self-mapping.

In fact, we shall prove a somewhat stronger result. Consider

for stipulated $x$. The function $\phi$ is continuous in the mixture $x_{\alpha }$. Define the mapping $N:x\to x^{\prime }$ by where the index $\alpha$ is taken to run over members of any finite set of admissible trials $x_{\alpha }$ in the execution of the Metropolis algorithm. (One may think of the collection of all sequences $x_{\alpha }$ in ensemble-theoretic terms.) Suppose that $x^{\prime }$ is a fixed point under Eq. (24). In Eq. (24), some values of $\alpha$ correspond to choices for $\{T_i,{\epsilon }_i(k)\}$ for which the posterior probability does not increase For these values of $\alpha$ If the choice $x$ is fixed under the mapping $N$ in Eq. (24), then the contribution to $x^{\prime }$ from any $x_{\beta }$ must not decrease; therefore, ${\phi }_{\beta }=0$, $\forall \beta$, lest the denominator in $\mathrm{\Phi }$ exceed unity. Put another way, no other choice of $\{T_i,{\epsilon }_i(k)\}$ can increase the posterior probability, but that is the definition of an equilibrium point. If, on the other hand, an equilibrium point $x$ maximizes the posterior probability [Eq. (17)], every ${\phi }_{\alpha }$ vanishes, so that $x$ is a fixed point.

Equation (24) is continuous and maps points $x$ into a convex compactum $\subset {\mathbf{R}}^{m+1}$. A fixed point

Therefore, exists according to the Brouwer fixed-point theorem that, by construction, maximizes the a-posteriori probability of $x$.

The mapping Eq. (19) generates a sequence of trials $x$ for which $p_n$ is nondecreasing and gives the maximal value of $p_n$ in the limit, while Eq. (27) demonstrates the existence of a trial $x^{*}$ for which $p_n$ cannot be made greater. In view of the ensemble-theoretic freedom to choose $x_{\alpha }$, we may identify the limit in Eq. (22) with the fixed point in Eq. (27). Therefore, a convergent subsequence of annealing trials exists that tends to an equilibrium point. Moreover, Eq. (27) demonstrates that the annealing search can, in principle, find $x^{*}$ in a finite number of trials. We conclude that, granted a suitable annealing schedule, there exists at least one convergent sequence of trials that tends to MAP surface temperature and spectral emissivity estimates consistent with observed spectral radiances $I_k$.

3.3.

Uniqueness

Sequential compactness guarantees the existence of a convergent subsequence of trials. In practice, we must expect that there will be more than one such sequence. The nonuniqueness of solutions to the TES problem suggests that there will be a continuous infinity of possible trials $\{T_n,{\epsilon }_i(k)\}$ that yield any stipulated value for the posterior probability. In any realizable search strategy, however, we need only contend with a countable set of convergent subsequences. Among these, there will be one for which the posterior probability is greatest. Its limit will be the closest approach to the MAP solution achieved by simulated annealing. In the nature of things, more than one convergent subsequence may be expected to exist that yields this same maximal estimate, with the same asymptotic annealing temperature $T_{\mathrm{a}}^{\infty }$. We ignore all subsequences except these maximal ones.

In Refs. 5 and 3, expectation values for $T$ and $\{\epsilon (k)\}$ over the the posterior probability [Eq. (15)] were shown to give good estimates for physical surface temperatures and emissivities. We claim that the mean of a large number of subsequences that converge to the limiting MAP value will tend to the expectation values for $T$ and $\{\epsilon (k)\}$ with respect to Eq. (15).

The MAXENT estimator is constructed from the posterior probability of noise power in a spectral bin. For the sake of simplicity, we assume identical noise power in each bin. (This assumption is inessential and may be relaxed.) A fully annealed MAP estimate may be thought of as an individual Bernoulli trial drawn from the likelihood function for $\{T_n,{\epsilon }_i(k)\}$. By construction, all such trials are independent and identically distributed with bounded expectation values. Moments over Eq. (15) are bounded despite bad behavior of the Jeffreys prior at $T=0$, because of the rapid decay of the exponentials away from the MAP solution, as L’Hôpital’s rule demonstrates.

Let

and be the means of MAP surface temperature and spectral emissivity taken over $N$ convergent subsequences. Suppose the covariance matrix $\mathbf{\Sigma }$ of the trials to be nonsingular. We invoke the multivariate central limit theorem to conclude the mean values converge weakly to the multivariate Gaussian distribution

Reliance on the mixing hypothesis in the form given by Eq. (10), however, brings with it the concern that the relevant covariance matrix might be singular. In that event, the strong law of large numbers^26,27 ensures

and but without giving estimated variances of the mean values, such as with Eq. (30).

To the extent that the estimator used in the simulated annealing search is zero-mean error, we conclude the estimates yield accurate values for the physical values of $T$ and $\{\epsilon (k)\}$. As the spectral weights $\mathbf{x}$, lying as they do between zero and unity, possess bounded moments, this conclusion applies to the limiting mean values of $\{T_n,\{x_m^n\}\}$ as well.