Testing that a Local Optimum of the Likelihood is Globally Optimum using Reparameterized Embeddings

TL;DR

This paper introduces a statistical test to verify if a local maximum of the likelihood function is actually the global maximum, using reparameterization to improve accuracy and computational efficiency in non-convex imaging problems.

Contribution

It proposes a novel reparameterization approach for likelihood functions to enhance global maximum testing in non-convex optimization problems.

Findings

Improved accuracy of global maximum testing through reparameterization.

Reduced computational cost in practical imaging applications.

Validated method on camera-blur estimation problem.

Abstract

Many mathematical imaging problems are posed as non-convex optimization problems. When numerically tractable global optimization procedures are not available, one is often interested in testing ex post facto whether or not a locally convergent algorithm has found the globally optimal solution. When the problem is formulated in terms of maximizing the likelihood function under a statistical model for the measurements, one can construct a statistical test that a local maximum is in fact the global maximum. A one-sided test is proposed for the case that the statistical model is a member of the generalized location family of probability distributions, a condition often satisfied in imaging and other inverse problems. We propose a general method for improving the accuracy of the test by reparameterizing the likelihood function to embed its domain into a higher dimensional parameter space. We show that the proposed global maximum testing method results in improved accuracy and reduced computation for a physically-motivated joint-inverse problem arising in camera-blur estimation.

Tables

Table 1. Table 1: Empirical performance of global convergence tests for an iterative quasi-Newton maximum likelihood estimator of the camera point spread function (100 Monte-Carlo trials). The probability of detection (PD) achieved by the proposed global maximum test ( 32 ) is perfect, while Biernacki's test [5] only attains 0.22 0.22 0.22 , when the probability of false alarm (PFA) is constrained to be 0.01 0.01 0.01 .

Test	PFA	PD
Biernacki biernacki:2005	$0.01$	$0.22$
Proposed (32)	$0.01$	$1.0$

Equations

$$\ell \left(\mathbf{d};\mathbf{\theta }\right)=\frac{1}{n}\sum _{k=1}^n\ln f\left({\mathbf{d}}_k;\mathbf{\theta }\right)$$

$${\hat{\mathbf{\theta }}}_{\mathrm{Global}}=\underset{\mathbf{\theta }\in \mathbf{\Theta }}{\mathrm{argmax}}\ell \left(\mathbf{d};\mathbf{\theta }\right).$$

$$H_0:\hat{\mathbf{\theta }}={\hat{\mathbf{\theta }}}_{\mathrm{Global}}vs.H_1:\hat{\mathbf{\theta }}\ne {\hat{\mathbf{\theta }}}_{\mathrm{Global}}.$$

$$d(x)$$

$$\mathbf{d}$$

$${\hat{\theta }}_{\mathrm{Global}}$$

$$\phi \left(\mathbf{d},\hat{\mathbf{\theta }}\right)\stackrel{\mathrm{def}}{=}\ell \left(\mathbf{d};\hat{\mathbf{\theta }}\right)-m(\hat{\mathbf{\theta }},\hat{\mathbf{\theta }}),$$

$$m({\mathbf{\theta }}_0,{\mathbf{\theta }}_1)={\mathrm{E}}_{{\mathbf{\theta }}_0}\left[\ell \left(\mathbf{D};{\mathbf{\theta }}_1\right)\right],$$

$$\frac{1}{\sqrt{n}}\sum _{k=1}^n\phi \left({\mathbf{D}}_k,\hat{\mathbf{\theta }}\right)\stackrel{D}{\to }\mathcal{N}\left(0,{\mathrm{Var}}_{{\mathbf{\theta }}_0}\left[\ell \left({\mathbf{D}}_1,{\mathbf{\theta }}_0\right)\right]\right),$$

$$\displaystyle\frac{\left(\ell\!\left({{\bm{{D}}};\widehat{{\bm{{\theta}}}}}\right)-m(\widehat{{\bm{{\theta}}}},\widehat{{\bm{{\theta}}}})\right)^{2}}{v(\widehat{{\bm{{\theta}}}})}\;\;\overset{H_{1}}{\underset{H_{0}}{\operatorname*{\begin{array}[]{c}>\vspace{-6pt}\\ <\end{array}}}}\;\;\eta,$$

$$-\sqrt{\eta v(\hat{\mathbf{\theta }})}\le \phi \left(\mathbf{D},\hat{\mathbf{\theta }}\right)\le \sqrt{\eta v(\hat{\mathbf{\theta }})}.$$

$$m(\hat{\mathbf{\theta }},\hat{\mathbf{\theta }})\ge m({\mathbf{\theta }}_0,\hat{\mathbf{\theta }}),$$

$$\displaystyle\frac{\ell\!\left({{\bm{{D}}};\widehat{{\bm{{\theta}}}}}\right)-m(\widehat{{\bm{{\theta}}}},\widehat{{\bm{{\theta}}}})}{\sqrt{v(\widehat{{\bm{{\theta}}}})}}\;\;\overset{H_{0}}{\underset{H_{1}}{\operatorname*{\begin{array}[]{c}>\vspace{-6pt}\\ <\end{array}}}}\;\;\eta_{1}.$$

$$\mathbf{D}=\mathbf{H}{\mathbf{\theta }}_0+\epsilon \text{ s.t. }\epsilon \sim \mathcal{N}\left(0,\mathbf{\Sigma }\right),$$

$$\ell \left(\mathbf{D};\mathbf{\theta }\right)=$$

$$-\frac{1}{2}\ln \left(\det \mathbf{\Sigma }\right)-\frac{q}{2}\ln \left(2\pi \right).$$

$$\lambda ={\left({\mathbf{\theta }}_0-\mathbf{\theta }\right)}^T{\mathbf{H}}^T{\mathbf{\Sigma }}^{-1}\mathbf{H}\left({\mathbf{\theta }}_0-\mathbf{\theta }\right).$$

$${\mathrm{E}}_{{\mathbf{\theta }}_0}\left[\ell \left(\mathbf{D};\mathbf{\theta }\right)\right]=-\frac{1}{2}\left(q+\lambda \right)-\frac{1}{2}\ln \left(\det \mathbf{\Sigma }\right)-\frac{q}{2}\ln \left(2\pi \right).$$

$$m(\mathbf{\theta },\mathbf{\theta })\ge m({\mathbf{\theta }}_0,\mathbf{\theta }).$$

$$m({\mathbf{\theta }}_0,\mathbf{\theta })=$$

$$=$$

$$m({\mathbf{\theta }}_0,{\mathbf{\theta }}_0)-m({\mathbf{\theta }}_0,\mathbf{\theta })=-\int f(\mathbf{d};{\mathbf{\theta }}_0)\log \frac{f(\mathbf{d};\mathbf{\theta })}{f(\mathbf{d};{\mathbf{\theta }}_0)}d\mu (\mathbf{d})$$

$$E_{{\mathbf{\theta }}_0}[\log f(\mathbf{D};{\mathbf{\theta }}_0)]\ge E_{{\mathbf{\theta }}_0}[\log f(\mathbf{D};\mathbf{\theta })].$$

$$m({\mathbf{\theta }}_0,{\mathbf{\theta }}_0)=$$

$$=$$

$$\tilde{g}(\mathbf{d};\mathbf{\theta },\tilde{\mathbf{\theta }})=\tilde{\ell }(\mathbf{d};\tilde{\mathbf{\theta }})-\ell (\mathbf{d};\mathbf{\theta }),\mathbf{\theta }\in \mathbf{\Theta },\tilde{\mathbf{\theta }}\in \mathbf{\Theta }\times {\mathbf{\Theta }}^{\prime }.$$

$$G(\mathbf{d},\mathbf{\theta })=\underset{\tilde{\mathbf{\theta }}\in \left\{\mathbf{\theta }\right\}\times {\mathbf{\Theta }}^{\prime }}{\max }\tilde{g}\left(\mathbf{d};\mathbf{\theta },\tilde{\mathbf{\theta }}\right).$$

$$m_G({\mathbf{\theta }}_0,\mathbf{\theta })=\underset{{\mathbf{\theta }}^{\prime }\in {\mathbf{\Theta }}^{\prime }}{\max }{E}_{{\mathbf{\theta }}_0}[\tilde{g}(\mathbf{D};\mathbf{\theta },(\mathbf{\theta },{\mathbf{\theta }}^{\prime }))].$$

$$\int f_0(\mathbf{d};{\mathbf{\theta }}_0)\log f_1(\mathbf{d};{\mathbf{\theta }}_0,{\mathbf{\theta }}^{\prime })d\mu (\mathbf{d}),$$

$$g(\mathbf{d},\hat{\mathbf{\theta }})=\underset{\tilde{\mathbf{\theta }}\in S_{\tilde{\ell }}\left(\hat{\mathbf{\theta }}\right)}{\max }\tilde{g}\left(\mathbf{d};\hat{\mathbf{\theta }},\tilde{\mathbf{\theta }}\right).$$

$$m_g({\mathbf{\theta }}_0,\mathbf{\theta })=\underset{\tilde{\mathbf{\theta }}\in S_{\mu }\left(\mathbf{\theta }\right)}{\max }{E}_{{\mathbf{\theta }}_0}[\tilde{g}(\mathbf{D};\mathbf{\theta },\tilde{\mathbf{\theta }})],$$

Full text

††journal: Journal of Mathematical Imaging and Vision

1 Introduction

Mathematical imaging problems are often formulated as non-convex energy minimization problems that impose desirable properties on the global optima, e.g., corresponding to a denoised, deblurred, or segmented image. Much of the work of Mila Nikolova addressed the problem of local and global optima. As stated succinctly in one of her early papers: ``The resultant … energy generally exhibits numerous local minima. Calculating its local minimum, placed in the vicinity of the maximum likelihood estimate, is inexpensive but inadequate'' nikolova1999markovian . Study of local and global optima was a recurring theme in her work, in which she addressed the nature of objective functions associated with non-convex probabilistic models, i.e., maximum likelihood (ML) and maximum a posteriori (MAP) nikolova1997estimees , nikolova1998inversion , nikolova1999markovian , nikolova2000segmentation , alberge2006blind , nikolova2007model , as well as non-linear least squares durand2006stabilityI , durand2006stabilityII . Some of the optimization algorithms she introduced were only shown to converge to one of several possible local optima. For such algorithms, an important question is whether an observed convergent limit is, in fact, the global maximum. Searching for and identifying the global maximum is the problem that we address in this paper.

We consider this problem in the general setting of maximum likelihood parameter estimation from multiple samples from a probability distribution that belongs to a parametric family. The conceptual simplicity and tractability of the Maximum Likelihood (ML) principle, along with its theoretical optimality properties, has made ML approaches prevalent in many fields. Yet, questions surrounding its practical application remain open. The pioneering statistician Sir Ronald Fisher fisher:1925 was an early advocate of the ML approach and is generally credited with its development, although similar concepts predate Fisher's work. Stigler stigler:2007 provides a historical account of the theory's maturation throughout the nineteenth and twentieth centuries. Asymptotic (large sample) characterization of local vs. global optima of the likelihood function was established by Le Cam (c.f. le-cam:2012a , ch. 6) using central limit theory, but can be challenging to apply in practice. As the sample size increases, it has long been known that, under mild smoothness conditions, all statistically consistent stationary points of the likelihood function converge with probability one to the global maximum cramer:1946 ; wald:1949 . The natural question then becomes: when there exist multiple stationary points, and the number of samples is finite, how can one identify the globally optimal one?

There exist general-purpose algorithms to address this question, e.g., simulated annealing and genetic algorithms. These algorithms, however, are rarely applied to high-dimensional problems because of high computational demands andrieu:2000 ; andrieu:2001 ; sharman:1989 . Stationary points of the likelihood function can be readily found using iterative root-finding methods such as Quasi-Newton gradient descent nocedal:1999 . Once a stationary point is found, it would be useful to have access to a simple test to determine if it is globally optimal without knowing the maximum value of the likelihood function. Several such tests have been proposed for this purpose biernacki:2005 ,blatt:2007 . In this paper, the focus is on testing local maxima of the likelihood function in the context of high dimensional inverse problems arising in signal processing and imaging.

Specifically, this paper makes the following contributions. Starting with the global maximum validation function introduced by Biernacki biernacki:2005 , we demonstrate that its mean is always less than or equal to zero when the likelihood function belongs to a generalized location family of distributions: distributions parameterized by a shift in location. This property provides the impetus for constructing a one-sided variant of the test. This generalized location family is relevant to many linear and non-linear inverse problems. Furthermore, we introduce a new approach of testing for the global maximum by expanding the parameter space to a higher dimension through a reparameterized embedding and defining an augmented validation function for testing local maxima. The augmented validation function can better discriminate between local and global maxima due to the expanded parameter space. We provide a computational procedure for identifying useful candidate embeddings that significantly improves the accuracy of the test. Significant accuracy and computational advantages are demonstrated for the application of camera blur-function estimation leblanc:2018 . In particular, when implemented as a multiple restart stopping criterion, the proposed global maximum testing procedure is shown to significantly reduce computation as compared to simulated annealing methods for finding global maxima.

The remainder of the paper is organized as follows. Section 2 describes the general global maximum testing problem and introduces a simple illustrative example used to demonstrate the key concepts. Section 3 introduces the one-sided global maximum test as a variant of the two-sided test of Biernacki biernacki:2005 , develops the reparameterized embedding method, and proposes a numerical spectral embedding procedure for identifying good embeddings. These concepts and procedures are illustrated in the context of a simple non-linear maximum likelihood estimation example. Section 4 illustrates the proposed methods for application to camera blur-function estimation.

2 Problem Description

2.1 Background

The problem setting is as follows. The observed data comes in the form of a matrix ${\bm{{d}}}=[{\bm{{d}}}_{1},\ldots,{\bm{{d}}}_{n}]\in\mathbb{R}^{m\times n}$ where the columns are independent and identically distributed (i.i.d.) realizations of an $m$-dimensional random vector ${\bm{{D}}}_{1}$ having a parametric Lebesgue density $f\!\left({{\bm{{d}}}_{1};{\bm{{\theta}}}}\right)$. The joint probability distribution $f({\bm{{d}}},{\bm{{\theta}}})=\prod_{k=1}^{n}f({\bm{{d}}}_{k},{\bm{{\theta}}})$ of ${\bm{{D}}}=[{\bm{{D}}}_{1},\ldots,{\bm{{D}}}_{n}]$ is a known function of ${\bm{{d}}}$ and ${\bm{{\theta}}}$, and is in a parametric family $\left\{f({\bm{{d}}};{\bm{{\theta}}}):{\bm{{d}}}\in\mathcal{D},{\bm{{\theta}}}\in{\bm{{\Theta}}}\right\}$. Here the vector of parameters ${\bm{{\theta}}}$ is unknown, taking values in a parameter space ${\bm{{\Theta}}}$, an open subset of $\mathbb{R}^{p}$, and the matrix of measurements ${\bm{{d}}}$ takes values in a sample space $\mathcal{D}\subset\mathbb{R}^{m\times n}$. For any mean square integrable function $X$ of ${\bm{{D}}}$ we define the statistical expectation $E_{{\bm{{\theta}}}_{0}}[X({\bm{{D}}})]=\int X({\bm{{d}}})f\!\left({{\bm{{d}}};{\bm{{\theta}}}_{0}}\right)d\mu({\bm{{d}}})$, where $d\mu({\bm{{d}}})$ indicates integration with respect to the Lebesgue measure on $\mathbb{R}^{m}$. The subscript ${\bm{{\theta}}}_{0}$ of the expectation operator $E_{{\bm{{\theta}}}_{0}}$ is called the true value of ${\bm{{\theta}}}$ as it parameterizes the underlying density $f({\bm{{d}}};{\bm{{\theta}}}_{0})$ generating the observations.

As the columns of the measurement matrix ${\bm{{d}}}$ are i.i.d. realizations, the log-likelihood function is

[TABLE]

The associated maximum likelihood estimator (MLE) $\widehat{{\bm{{\theta}}}}:\mathcal{D}\rightarrow{\bm{{\Theta}}}$ is defined as the global maximum

[TABLE]

An estimator $\widehat{{\bm{{\theta}}}}$ is said to be consistent (statistically consistent in norm) when $\lim_{n\rightarrow\infty}E_{{\bm{{\theta}}}_{0}}[\|\widehat{{\bm{{\theta}}}}-{\bm{{\theta}}}_{0}\|^{2}]\rightarrow 0$. We will assume that $f$ is continuously differentiable in ${\bm{{\theta}}}$ for all ${\bm{{d}}}$, and define the score function as ${\bm{{s}}}\!\left({{\bm{{d}}},{\bm{{\theta}}}}\right)=\nabla_{{\bm{{\theta}}}}\ell\!\left({{\bm{{d}}},{\bm{{\theta}}}}\right)$. The Fisher information matrix

${\bm{{I}}}\!\left({{\bm{{\theta}}}_{0}}\right)=E_{{\bm{{\theta}}}_{0}}[{\bm{{s}}}\!\left({{\bm{{D}}},{\bm{{\theta}}}_{0}}\right){\bm{{s}}}\!\left({{\bm{{D}}},{\bm{{\theta}}}_{0}}\right)^{T}]$ is assumed to exist and be invertible, and the notation $\overset{P}{\to}$''and $\overset{D}{\to}$'' will be used to describe convergence in probability and distribution respectively.

For the problem addressed in this paper, the global maximum of the log-likelihood is unknown, and only a local maximum $\widehat{{\bm{{\theta}}}}$ is available, which is not necessarily equal to $\widehat{{\bm{{\theta}}}}_{\mathrm{Global}}$. For example, the local maximum could be the limit of a convergent gradient descent algorithm. Given $\widehat{{\bm{{\theta}}}}$, the local maximum testing problem is to decide between the two hypotheses

[TABLE]

A test between $H_{0}$ and $H_{1}$ is defined as a binary valued function $\phi:\mathcal{D}\rightarrow\left\{0,1\right\}$ that maps the data ${\bm{{d}}}$ to [math] or $1$, indicating the decision $H_{0}$ or $H_{1}$, respectively. The accuracy of a test is measured by its probability of false alarm PFA=$E_{{\bm{{\theta}}}_{0}}[\phi|H_{0}]$ and its probability of detection PD=$E_{{\bm{{\theta}}}_{0}}[\phi|H_{1}]$. If for two tests $\phi_{1}$ and $\phi_{2}$ having identical PFA, PD of $\phi_{1}$ is greater than PD of $\phi_{2}$, then $\phi_{1}$ is said to be more powerful than $\phi_{2}$.

Many approaches to the general hypothesis testing problem (3) have been studied over the years. Blatt and Hero blatt:2007 presented a historical context, which is summarized here. The likelihood ratio test wilks:1938 , Wald test wald:1943 , and Rao score test rao:1948 are asymptotically equivalent tests as the number $n$ of samples approaches infinity. The likelihood ratio and Wald tests require the distribution under $H_{0}$ to be known, which for (3) requires knowledge of the true parameter. On the other hand, the Rao score test, later independently discovered and popularized under the name Lagrange multiplier test silvey:1959 , can be implemented when the true parameter is unknown. Rao's test measures the Euclidean norm of the score function weighted by the inverse Fisher information evaluated at a local maximum $\xi_{R}=\frac{1}{p}{\bm{{s}}}\!\left({{\bm{{d}}},\widehat{{\bm{{\theta}}}}}\right)^{T}{\bm{{I}}}^{-1}\!\left({\widehat{{\bm{{\theta}}}}}\right){\bm{{s}}}\!\left({{\bm{{d}}},\widehat{{\bm{{\theta}}}}}\right)$. Gan and Jiang gan:1999 propose a similar test for consistency of a stationary point of the log-likelihood based on White's information test white:1982 . White's original work was concerned with testing for model misspecification under the assumption that the global maximum of the likelihood function had been located, and Gan uses the same test statistic but in the converse situation.

The Rao test may be used to test for consistency of a local maximum of the log-likelihood function. Unfortunately, Monte Carlo experiments indicate that this test may not be very powerful even in the univariate setting gan:1999 ; biernacki:2005 . In biernacki:2005 an improved test was proposed for testing consistency of a stationary point following ideas presented by Cox cox:1961 ; cox:1962 . This test, called the Biernacki test, uses a bootstrap estimate to directly compare the observed value of the locally maximized log-likelihood to its statistical expectation. Both the Rao score and the Biernacki tests fall under the more general M-testing framework described by Blatt and Hero blatt:2007 , where additional types of tests of local maxima are proposed.

2.2 Motivating Example

To illustrate the difficulties in testing local maxima of the log-likelihood consider the following one dimensional statistical estimation problem. Let $x\in[0,T]$ be a time interval and $x_{i}=iT/N$, $i=0,\ldots,N-1$. The measurements $\left\{d(x_{i})\right\}_{i=1}^{N}$ are a set of time samples of a sinusoidal signal in additive Gaussian noise

[TABLE]

Here $\theta_{0}\geq 0$ is an unknown sinusoidal frequency parameter to be estimated. More generally, it will be more convenient to express the measurement model in vector form ${\bm{{d}}}=[d(x_{0}),\ldots,d(x_{N-1})]^{T}$

[TABLE]

where ${\bm{{\mu}}}(\theta_{0})=\mathrm{sin}(\theta_{0}{\bm{{x}}})$ is the mean of ${\bm{{d}}}$, a vector of time samples of the noiseless signal, and ${\bm{{\varepsilon}}}$ is an independent identically distributed (i.i.d.) zero mean $N$-dimenstional Gaussian noise with identity covariance matrix scaled by the variance parameter $\sigma^{2}$, which is assumed known. The maximum likelihood estimator $\widehat{\theta}$ of the frequency parameter is then the globally optimal solution

[TABLE]

For the sinusoidal signal in noise model (4), there will be local maxima of the log-likelihood function corresponding to the multiple stationary points of $\left\|{\bm{{\mu}}}\!\left({\theta}\right)-{\bm{{d}}}\right\|^{2}$ over $\theta$. Figure 1(a) shows two of these local maxima for the case that the noise variance $\sigma^{2}$ is zero. The solid curve in Figure 1(a) corresponds to the true signal ${\bm{{\mu}}}(\theta)$, where $\theta=\theta_{0}$ is the global maximum, and the dashed curve corresponds to another signal for which $\theta$ is a local (non-global) maximum. Also shown is noisy data ${\bm{{d}}}$ generated with each of these two signals, corresponding to blue circles and red crosses, respectively, where the Gaussian noise variance is $\sigma^{2}=1$. The difficulty in perceiving differences between these two noisy signals suggests that distinguishing a sub-optimal local maximum from the global maximum will be challenging. This becomes even more difficult for higher dimensional estimation problems occurring in imaging (c.f. leblanc:2018  Section 3c).

3 Tests for Local Optima

In the hypothesis testing problem described by (3), the null-hypothesis $H_{0}$ is that the discovered local maximum $\widehat{{\bm{{\theta}}}}$ of the log-likelihood is a global maximum. It is important to note that a failure to reject the null hypothesis is not a positive statement about the global optimality of $\widehat{{\bm{{\theta}}}}$. Instead, when a test accepts $H_{0}$, all that can be said is that it does not rule out the point as a local maximum with sufficient statistical certainty.

3.1 A Two-Sided Test

To test whether a local maximum of the likelihood function is, in fact, the global maximum one defines a suitable global maximum validation function whose statistical distribution changes depending on whether the local maximum $\widehat{{\bm{{\theta}}}}$ is global or not blatt:2007 . Define the validation function

[TABLE]

where $m:\mathbb{R}^{p}\times\mathbb{R}^{p}\rightarrow\mathbb{R}$ is the mean function

[TABLE]

This function is called the ambiguity function and is the statistical expectation under the distribution $f({\bm{{d}}};{\bm{{\theta}}}_{0})$ of the log-likelihood function $\ell({\bm{{D}}},{\bm{{\theta}}})$ evaluated at ${\bm{{\theta}}}={\bm{{\theta}}}_{1}$. Assuming that the global maximum $\widehat{{\bm{{\theta}}}}_{\mathrm{Global}}$ is near the true value ${\bm{{\theta}}}_{0}$, under the null hypothesis $H_{0}$ we have $\widehat{{\bm{{\theta}}}}=\widehat{{\bm{{\theta}}}}_{\mathrm{Global}}$, and the distribution of $\varphi\!\left({{\bm{{D}}},\widehat{{\bm{{\theta}}}}}\right)$ will have approximately zero mean. On the other hand, under the alternative hypothesis $H_{1}$ that $\widehat{{\bm{{\theta}}}}$ is a non-global local maximum the mean of the distribution of $\varphi\!\left({{\bm{{D}}},\widehat{{\bm{{\theta}}}}}\right)$ will shift away from zero. This is the key motivation for using the validation function (7) to test for a global maximum.

For an i.i.d. data sample ${\bm{{D}}}_{1},\ldots,{\bm{{D}}}_{n}$ the following asymptotic result was established in (biernacki:2005, , Theorem 2). Under $H_{0}$:

[TABLE]

where ${\mathrm{Var}}\!\left[\ell\!\left({{\bm{{D}}}_{1},{\bm{{\theta}}}_{0}}\right)\right]$ is the variance of the log-likelihood function for a single data sample ($n=1$). Recalling the definition of the random data matrix ${\bm{{D}}}=\left\{{\bm{{D}}}_{1},\ldots,{\bm{{D}}}_{n}\right\}$, this Gaussian limit motivates us to define the following test of the hypotheses $H_{0}:\widehat{{\bm{{\theta}}}}=\widehat{{\bm{{\theta}}}}_{Global}$ vs. $H_{1}:\widehat{{\bm{{\theta}}}}\not=\widehat{{\bm{{\theta}}}}_{Global}$

[TABLE]

where the function $v:\mathbb{R}^{p}\rightarrow\mathbb{R}$ is the variance $v(\theta)={\mathrm{Var}}_{{\bm{{\theta}}}}\!\left[\ell\!\left({{\bm{{D}}};{\bm{{\theta}}}}\right)\right]$ under the distribution $f({\bm{{d}}},{\bm{{\theta}}})$ of the log-likelihood evaluated at ${\bm{{\theta}}}$, and $\eta$ is a threshold selected to fix the false alarm probability equal to a suitably small number $\alpha\in[0,1]$. Under local asymptotically normal (LAN) conditions on the likelihood function lehmann:1998 $\widehat{{\bm{{\theta}}}}\overset{P}{\to}{\bm{{\theta}}}_{0}$ (a.s.) and the test statistic on the left hand side of (12) has an approximately chi-square distribution under $H_{0}$. Hence $\eta$ can be selected as the $1-\alpha$ quantile of the chi-square distribution. In biernacki:2005 this test was implemented by approximating the mean ${\mathrm{E}}_{{\bm{{\theta}}}_{0}}\!\left[\ell\!\left({{\bm{{D}}};{\bm{{\theta}}}_{0}}\right)\right]$ and the variance ${\mathrm{Var}}\!\left[\ell\!\left({{\bm{{D}}};{\bm{{\theta}}}_{0}}\right)\right]$ using a parametric bootstrap estimator.

The test (12) is called a two-sided test because the condition for which the null hypothesis is accepted can be equivalently be expressed as

[TABLE]

This is thus a test for which, as compared to the global maximum $\widehat{{\bm{{\theta}}}}_{\mathrm{Global}}$, a sub-optimal local maximum $\widehat{{\bm{{\theta}}}}$ will cause the test function to undergo a shift in mean, where the shift could either be in a positive or a negative direction.

3.2 A One-Sided Test

If it were known a priori that a sub-optimal local maximum causes a negative shift in the mean of the global maximum validation function $\varphi\!\left({{\bm{{D}}},\widehat{{\bm{{\theta}}}}}\right)$, a one-sided test would be advantageous over a two-sided test. More specifically, a one-sided test would be expected to have higher power than the two-sided test (12) when for all $\widehat{{\bm{{\theta}}}}\neq{\bm{{\theta}}}_{0}$,

[TABLE]

where $m({\bm{{\theta}}}_{0},{\bm{{\theta}}}_{1})$ is the ambiguity function defined in (8). When this condition is satisfied the two-sided test (12) can be replaced by the one-sided test

[TABLE]

The condition (13) is satisfied for many imaging and inverse problems. For example, consider the case where ${\bm{{\theta}}}$ is a clean image that one wishes to recover from samples of ${\bm{{D}}}$, the output of an imaging sensor with known point spread function (forward operator) in additive correlated noise. When the point-spread function (PSF) and the covariance are known, we will show that this model always satisfies the inequality (13), and the one-sided test might be expected to lead to a better test for a global maximum. Define ${\bm{{\theta}}}_{0}\in\mathbb{R}^{p}$ as the vectorized true image to be recovered and ${\bm{{D}}}\in\mathbb{R}^{q}$ as the random vectorized image acquired from the camera, which obeys the model:

[TABLE]

where ${\bm{{H}}}$ is a $q\times p$ matrix representing the forward operator and ${\bm{{\Sigma}}}$ is the $q\times q$ camera noise covariance matrix.

To show that (13) holds in this case, start with the log-likelihood function for the above model

[TABLE]

For any value of ${\bm{{\theta}}}$, (18) is a quadratic form in $\varepsilon$ that is distributed non-central chi-squared with non-centrality parameter

[TABLE]

The moment properties of the non-central chi-squared distribution kotz2004continuous thus specify the statistical expectation of the log-likelihood function (18) :

[TABLE]

The difference ${\mathrm{E}}_{{\bm{{\theta}}}}\!\left[\ell\!\left({{\bm{{D}}};{\bm{{\theta}}}}\right)\right]-{\mathrm{E}}_{{\bm{{\theta}}}_{0}}\!\left[\ell\!\left({{\bm{{D}}};{\bm{{\theta}}}}\right)\right]=m({\bm{{\theta}}},{\bm{{\theta}}})-m({\bm{{\theta}}}_{0},{\bm{{\theta}}})=\lambda/2$, is non-negative, establishing that (13) holds as claimed. For this example, the unconstrained maximum likelihood estimator of ${\bm{{\theta}}}$ is a solution to a convex optimization problem, which is strictly convex when ${\bm{{H}}}$ is full column-rank, and thus there will be no sub-optimal isolated local maxima of (2). As our simple example in Figure 1 illustrated, additional constraints can give rise to local maxima.

The condition (13) is satisfied for a more general class of camera models where the probability distribution of the data is in the generalized location family.

Definition 1

Let $f({\bm{{d}}};{\bm{{\theta}}})$ be a distribution defined on ${\bm{{d}}}\in\mathbb{R}^{m}$ parameterized by ${\bm{{\theta}}}\in{\bm{{\Theta}}}\subset\mathbb{R}^{p}$. The distribution belongs to the generalized location family of distributions if there exists a function ${\bm{{g}}}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{m}$ such that $f\!\left({{\bm{{d}}};{\bm{{\theta}}}}\right)=f\!\left({{\bm{{d}}}-{\bm{{g}}}\!\left({{\bm{{\theta}}}}\right)}\right)$ for all ${\bm{{\theta}}}\in\Theta$ and all ${\bm{{x}}}\in\mathbb{R}^{m}$.

Any camera model of the form ${\bm{{D}}}_{k}={\bm{{\mu}}}({\bm{{\theta}}}_{0})+{\bm{{\varepsilon}}}_{k}$, $k=1,\ldots,n$ where ${\bm{{\mu}}}(\cdot)$ is a possibly non-linear function and ${\bm{{\varepsilon}}}_{k}$ is i.i.d. but possibly non-Gaussian noise, will have a distribution that is in the generalized location family.

Theorem 3.1

Let ${\bm{{D}}}_{1},\ldots,{\bm{{D}}}_{n}$ be an i.i.d. sample and assume that ${\bm{{D}}}_{1}$ has distribution $f\!\left({{\bm{{d}}}_{1};{\bm{{\theta}}}}\right)$ belonging to a generalized location family. Then the inequality (13) holds.

Proof

The proof of the Theorem proceeds in two parts. The first part establishes that, for any parameters ${\bm{{\theta}}}$ and ${\bm{{\theta}}}_{0}$, the ambiguity function (8) satisfies: $m({\bm{{\theta}}}_{0},{\bm{{\theta}}}_{0})\geq m({\bm{{\theta}}}_{0},{\bm{{\theta}}})$ (Claim 1). The second part establishes that, for $f({\bm{{d}}};{\bm{{\theta}}})$ in a generalized location family, $m({\bm{{\theta}}}_{0},{\bm{{\theta}}}_{0})=m({\bm{{\theta}}},{\bm{{\theta}}})$ (Claim 2). Putting these two parts together implies

[TABLE]

The theorem then follows upon specialization of this inequality to ${\bm{{\theta}}}=\widehat{{\bm{{\theta}}}}$.

We recall the integral form for the mean function

[TABLE]

and the identity $E_{{\bm{{\theta}}}_{0}}[\log f({\bm{{D}}};{\bm{{\theta}}})]=nE_{{\bm{{\theta}}}_{0}}[\log f({\bm{{D}}}_{1};{\bm{{\theta}}})]$, which follows from the i.i.d. assumption.

Claim 1 in this proof follows from the non-negativity property of the Kullback-Liebler (KL) divergence, a well known result in statistics and information theory kullback1997information . For completeness, we give a self contained proof. Start with the expression:

[TABLE]

Now, using the elementary inequality $\log(u)\leq u-1$ and the fact that $\int f({\bm{{d}}},\tilde{{\bm{{\theta}}}})d\mu({\bm{{d}}})=1$ for all $\tilde{{\bm{{\theta}}}}$, the right hand side of (LABEL:eq:KL) is non-negative. Therefore, using the definition (21), this establishes the claim

[TABLE]

Claim 2 of this proof is a direct result of $f({\bm{{d}}}_{k};{\bm{{\theta}}})$ being in the generalized location family. Specifically,

[TABLE]

where the second equality comes from the generalized location family definition and the third equality follows from making the change of variable of integration $\tilde{{\bm{{d}}}}_{1}={\bm{{d}}}_{1}+(g({\bm{{\theta}}}_{0})-g({\bm{{\theta}}}))$. This establishes the Theorem. ∎

3.3 Reparameterized Embeddings

The detection performance of the one-sided and two-sided tests for the global maximum depends on how much shift a local maximum $\widehat{{\bm{{\theta}}}}\in{\bm{{\Theta}}}$ causes in the distribution of the associated validation function (7) introduced in Sec. 3.1. In this section, we explain how embedding the parameters into a higher dimensional parameter space and a modified validation function can increase this shift, leading to improved detection performance.

The proposed approach can be viewed as analogous to the advantageous use of higher dimensional embeddings in other mathematical problems. For example, in mathematical imaging, the level-set method for image segmentation osher2004level ; sethian1996fast embeds a two-dimensional curve in the plane as a level set of a higher-dimensional parameterized surface. As another example, parameter expansion is applied to stabilize numerical solutions to non-linear differential equations wang2008nonlinear . A similar approach is used in computational statistics for accelerating the convergence of parameter estimates in the iterative parameter expansion expectation-maximization (PX-EM) algorithm liu1998parameter . Furthermore, in machine learning, the support vector machine (SVM) vapnik1998support improves classification performance by representing the decision region in the native lower dimensions as a separating hyperplane in a much higher dimensional space. In analogy to the above examples, the reparameterization embedding we propose will lead to significant improvements in testing for the global optimum.

We define the reparameterized embedding as follows. As above let the log-likelihood function $\ell\!\left({{\bm{{\theta}}}}\right)$ be parameterized by a native parameter ${\bm{{\theta}}}\in{\bm{{\Theta}}}$. Let ${\bm{{\theta}}}^{\prime}\in{\bm{{\Theta}}}^{\prime}$ be a fictitious parameter and define the expanded parameterization $\tilde{{\bm{{\theta}}}}=({\bm{{\theta}}},{\bm{{\theta}}}^{\prime})$ living in $\tilde{{\bm{{\Theta}}}}={\bm{{\Theta}}}\times{\bm{{\Theta}}}^{\prime}$. The native parameter ${\bm{{\theta}}}$ is thus embedded in the cylinder set $\left\{{\bm{{\theta}}}\right\}\times{\bm{{\Theta}}}^{\prime}$ of the higher dimensional space $\tilde{{\bm{{\Theta}}}}$. Associated with this embedding, define the augmented log-likelihood function $\tilde{\ell}({\bm{{d}}};\tilde{{\bm{{\theta}}}})$ parameterized by $\tilde{{\bm{{\theta}}}}\in\tilde{{\bm{{\Theta}}}}$. We can link the embedded parameterization $\tilde{{\bm{{\theta}}}}=({\bm{{\theta}}},{\bm{{\theta}}}^{\prime})\in\tilde{{\bm{{\Theta}}}}={\bm{{\Theta}}}\times{\bm{{\Theta}}}^{\prime}$ to the native parameterization ${\bm{{\theta}}}\in{\bm{{\Theta}}}$ by fixing ${\bm{{\theta}}}^{\prime}$, which we can assume is equal to ${{\bm{{0}}}}$ without loss of generality. Thus the native parameterization is equal to a cross-section of the embedded parameterization $\left\{\tilde{{\bm{{\theta}}}}=({\bm{{\theta}}},{\bm{{0}}}):{\bm{{\theta}}}\in{\bm{{\Theta}}}\right\}$ which gives the relation between the native and augmented log-likelihood functions: $\ell({\bm{{d}}};{\bm{{\theta}}})=\tilde{\ell}({\bm{{d}}};{\bm{{\theta}}},{\bm{{0}}})$.

Define the difference between the augmented log-likelihood function and the native log-likelihood function:

[TABLE]

We first introduce a prototype for a new class of global maximum validation functions introduced in this section. This prototype, called the max gap function, is defined as the maximum of $\tilde{g}$ with respect to $\tilde{{\bm{{\theta}}}}$ over the restricted embedding space $\left\{{\bm{{\theta}}}\right\}\times{\bm{{\Theta}}}^{\prime}$

[TABLE]

When ${\bm{{\theta}}}=\widehat{{\bm{{\theta}}}}$ is a non-global maximum of the native log-likelihood function $\ell({\bm{{d}}};{\bm{{\theta}}})$, $G({\bm{{d}}},\widehat{{\bm{{\theta}}}})$ measures the gap between the augmented log-likelihood $\tilde{\ell}({\bm{{d}}};(\widehat{{\bm{{\theta}}}},{\bm{{\theta}}}^{\prime}))$ maximized over ${\bm{{\theta}}}^{\prime}\in{\bm{{\Theta}}}^{\prime}$ and $\ell({\bm{{d}}};{\bm{{\theta}}})$ evaluated at ${\bm{{\theta}}}=\widehat{{\bm{{\theta}}}}$. We would like to design the reparameterized embedding to maximize this gap. A mean gap function would provide a criterion for doing this.

For native parameters ${\bm{{\theta}}}_{0},{\bm{{\theta}}}\in{\bm{{\Theta}}}$ define the max mean gap function of ${\bm{{\theta}}}\in{\bm{{\Theta}}}$

[TABLE]

The mean of the max gap function is lower bounded by the max mean gap function.

Theorem 3.2

For any ${\bm{{\theta}}}_{0}$ and ${\bm{{\theta}}}$ in ${\bm{{\Theta}}}$ the max gap function (25) and the max mean gap function (26) are non-negative. Furthermore, $E_{{\bm{{\theta}}}_{0}}[G({\bm{{D}}};{\bm{{\theta}}})]\geq m_{G}({\bm{{\theta}}}_{0},{\bm{{\theta}}})$ and $m_{G}({\bm{{\theta}}}_{0},{\bm{{\theta}}}_{0})=0$.

Proof: Non-negativity of $G({\bm{{d}}};{\bm{{\theta}}})$ follows from the relation $\max_{{\bm{{\theta}}}^{\prime}\in{\bm{{\Theta}}}^{\prime}}\tilde{g}({\bm{{d}}};{\bm{{\theta}}},({\bm{{\theta}}},{\bm{{\theta}}}^{\prime}))\geq\tilde{g}({\bm{{d}}};{\bm{{\theta}}},({\bm{{\theta}}},\bm{0}))$, which is equal to zero since $\tilde{\ell}({\bm{{d}}};({\bm{{\theta}}},\bm{0}))=\ell({\bm{{d}}};{\bm{{\theta}}})$. An identical argument establishes non-negativity of $m_{G}({\bm{{\theta}}}_{0},{\bm{{\theta}}})$. The stated inequality follows from the fact that, for any random variable $Y(u)$ depending on a non-random parameter $u$, $\max_{u}E[Y(u)]$ is no greater than $E[\max_{u}Y(u)]$.

The proof that $m_{G}({\bm{{\theta}}}_{0},{\bm{{\theta}}}_{0})=0$ uses an argument similar to what we used to establish Claim 1 in the proof of Theorem 1. As $E_{{\bm{{\theta}}}_{0}}[\ell({\bm{{D}}};{\bm{{\theta}}}_{0})]=m({\bm{{\theta}}}_{0},{\bm{{\theta}}}_{0})$, $m_{G}({\bm{{\theta}}}_{0},{\bm{{\theta}}}_{0})=\max_{{\bm{{\theta}}}^{\prime}\in{\bm{{\Theta}}}^{\prime}}E_{{\bm{{\theta}}}_{0}}[\tilde{\ell}({\bm{{D}}};({\bm{{\theta}}}_{0},{\bm{{\theta}}}^{\prime})]-m({\bm{{\theta}}}_{0},{\bm{{\theta}}}_{0}).$ $E_{{\bm{{\theta}}}_{0}}[\tilde{\ell}({\bm{{D}}};({\bm{{\theta}}}_{0},{\bm{{\theta}}}^{\prime}))]$ is of the form

[TABLE]

where the densities $f_{0}$ and $f_{1}$ satisfy $f_{1}({\bm{{d}}};{\bm{{\theta}}}_{0},{\bm{{0}}})=f_{0}({\bm{{d}}},{\bm{{\theta}}}_{0})$. Therefore, invoking the non-negativity of the KL divergence, which led to inequality (23) in the proof of Theorem 1, the integral above takes its maximum over ${\bm{{\theta}}}^{\prime}\in{\bm{{\Theta}}}^{\prime}$ when ${\bm{{\theta}}}^{\prime}=\bm{0}$ and by definition (8) this value is equal to $m({\bm{{\theta}}}_{0},{\bm{{\theta}}}_{0})$. ∎

If the prototype validation function $G({\bm{{d}}};\widehat{{\bm{{\theta}}}})$ defined in (25) were to be used to test a local maximum $\widehat{{\bm{{\theta}}}}$, Theorem 2 suggests that the function $m_{G}({\bm{{\theta}}}_{0},{\bm{{\theta}}})$ defined in (26) can be used to design reparameterized embeddings that induce the largest possible positive shifts in $G({\bm{{\theta}}}_{0},{\bm{{\theta}}})$ when ${\bm{{\theta}}}$ deviates from the true parameter ${\bm{{\theta}}}_{0}$, which would be expected to occur if ${\bm{{\theta}}}=\widehat{{\bm{{\theta}}}}$ were a non-global maximum of $\ell({\bm{{d}}};{\bm{{\theta}}})$. We describe a spectral embedding procedure below that implements such a design strategy.

While the prototype validation function $G$ in (25) has the benefits of simplicity and the analytical lower bound in Thm 2, it has two deficiencies that motivate an alternative gap function, called $g$ and defined in (27) below. First, the maximization in $G$ is restricted to the sub-space of embedded parameter values $\tilde{{\bm{{\theta}}}}=({\bm{{\theta}}},{\bm{{\theta}}}^{\prime})\in\tilde{{\bm{{\Theta}}}}$ over which ${\bm{{\theta}}}=\widehat{{\bm{{\theta}}}}$ is fixed. This restriction deprives the prototype of extra degrees of freedom, potentially reducing its sensitivity to non-global maxima. Second, the evaluation of $G$ requires a global maximization over ${\bm{{\theta}}}^{\prime}\in{\bm{{\Theta}}}^{\prime}$, which may be challenging in practice.

We next introduce an alternative gap function that overcomes these deficiencies: it expands the maximization in (25) to the full embedding space $\tilde{{\bm{{\Theta}}}}={\bm{{\Theta}}}\times{\bm{{\Theta}}}^{\prime}$ but allows the global maximum to be replaced by a local maximum found by an iterative algorithm initialized at $\tilde{{\bm{{\theta}}}}=(\widehat{{\bm{{\theta}}}},{\bm{{0}}})$.

Definition 2

Given a point ${\bm{{\theta}}}\in{\bm{{\Theta}}}$ and an iterative algorithm for finding local maxima of a function $f(\tilde{{\bm{{\theta}}}})$, $\tilde{{\bm{{\theta}}}}\in\tilde{{\bm{{\Theta}}}}={\bm{{\Theta}}}\times{\bm{{\Theta}}}^{\prime}$, the basin of attraction $S_{f}\!\left({{\bm{{\theta}}}}\right)\subseteq\tilde{{\bm{{\Theta}}}}$ is defined as an open set containing the point $\tilde{{\bm{{\theta}}}}=({\bm{{\theta}}},{\bm{{0}}})$ such that the algorithm converges to a local maximum of $f(\tilde{{\bm{{\theta}}}})$ when initialized at that point.

Let $\widehat{{\bm{{\theta}}}}\in{\bm{{\Theta}}}$ be a local maximum of the native log-likelihood function $\ell({\bm{{d}}};{\bm{{\theta}}})$ and let $\widehat{\tilde{{\bm{{\theta}}}}}\in S_{\tilde{\ell}}\!\left({\widehat{{\bm{{\theta}}}}}\right)\subseteq\tilde{{\bm{{\Theta}}}}$ be a local maximum of the augmented log-likelihood function $\tilde{\ell}({\bm{{d}}};\tilde{{\bm{{\theta}}}})$ in a basin of attraction containing $\tilde{{\bm{{\theta}}}}=(\widehat{{\bm{{\theta}}}},{{\bm{{0}}}})$. As an alternative gap function to (25) we define the augmented validation function

[TABLE]

This measures the gap between the value of the native log-likelihood function evaluated at $\widehat{{\bm{{\theta}}}}$ and the maximum of the augmented log-likelihood function in the neighborhood of $(\widehat{{\bm{{\theta}}}},{\bm{{0}}})$.

As an analog to (26), for native parameters ${\bm{{\theta}}}_{0},{\bm{{\theta}}}\in{\bm{{\Theta}}}$, define the locally maximized mean of this gap function, called the augmented ambiguity function

[TABLE]

where for a point ${\bm{{\theta}}}\in{\bm{{\Theta}}}$, $S_{\mu}\!\left({{\bm{{\theta}}}}\right)$ denotes a basin of attraction of the mean function $\mu(\tilde{{\bm{{\theta}}}})=E_{{\bm{{\theta}}}_{0}}[\tilde{\ell}({\bm{{D}}};\tilde{{\bm{{\theta}}}})]$ containing the point $\tilde{{\bm{{\theta}}}}=({\bm{{\theta}}},{{\bm{{0}}}})$.

Theorem 3.3

*For any ${\bm{{\theta}}}_{0}$ and ${\bm{{\theta}}}$ in ${\bm{{\Theta}}}$ *

[TABLE]

The leftmost inequality is achieved with equality when ${\bm{{\theta}}}={\bm{{\theta}}}_{0}$.

Proof: The proof is similar to the proof of Theorem 2 and relies on the fact that $({\bm{{\theta}}},{\bm{{0}}})\in S_{\mu}({\bm{{\theta}}})$, by definition of the bassin of attraction $S_{\mu}({\bm{{\theta}}})$. In particular, this implies that $m_{g}({\bm{{\theta}}}_{0},{\bm{{\theta}}})\geq E_{{\bm{{\theta}}}_{0}}[\tilde{\ell}({\bm{{D}}};({\bm{{\theta}}},{\bm{{0}}}))-\ell({\bm{{D}}};{\bm{{\theta}}})]=0,$ establishing the leftmost inequality in (29). The condition for equality in the leftmost inequality follows from the fact that $E_{{\bm{{\theta}}}_{0}}[\tilde{\ell}({\bm{{D}}};\tilde{{\bm{{\theta}}}})]\leq m({\bm{{\theta}}}_{0},{\bm{{\theta}}}_{0})$, achieving equality when $\tilde{{\bm{{\theta}}}}=({\bm{{\theta}}}_{0},{{\bm{{0}}}})$. The rightmost inequality in (29) follows from the same property of expectation of maximized random variables as used in the proof of Theorem 2. ∎

Under the condition that the basins of attraction $S_{\tilde{\ell}}(\widehat{{\bm{{\theta}}}})$ and $S_{\mu}(\widehat{{\bm{{\theta}}}})$ are equal, Theorem 3 would give the lower bound $m_{g}({\bm{{\theta}}}_{0},\widehat{{\bm{{\theta}}}})$ on the mean shift in $g({\bm{{d}}},\widehat{{\bm{{\theta}}}})$. If the augmented log-likelihood is smooth and the data ${\bm{{D}}}=\left\{{\bm{{D}}}_{1},\ldots,{\bm{{D}}}_{N}\right\}$ consists of i.i.d. samples, one can expect this condition to be satisfied asymptotically in $N$ since by the law of large numbers $\tilde{\ell}({\bm{{D}}},\tilde{{\bm{{\theta}}}})$ converges almost surely to its mean $E_{\theta_{0}}[\tilde{\ell}({\bm{{D}}},\tilde{{\bm{{\theta}}}})]$. Under such conditions, similarly to Theorem 2 for the gap function $G$ (25), Theorem 3 can be used to justify the use of the mean function $m_{g}$ to explore candidate reparameterized embeddings using the augmented validation function $g$ (27).

In analogy to (16), for gap functions $G$ or $g$, we propose a one sided reparameterized embedding test of $H_{0}:\widehat{{\bm{{\theta}}}}=\widehat{{\bm{{\theta}}}}_{Global}$ vs. $H_{1}:\widehat{{\bm{{\theta}}}}\not=\widehat{{\bm{{\theta}}}}_{Global}$. For the augmented validation function $g$ the proposed test is of the form:

[TABLE]

and similarly for the max gap function $G$. As above, for ${\bm{{\theta}}}_{0},{\bm{{\theta}}}_{1}\in{\bm{{\Theta}}}$, $m_{g}\!\left({{\bm{{\theta}}}_{0},{\bm{{\theta}}}_{1}}\right)={\mathrm{E}}_{{\bm{{\theta}}}_{0}}\!\left[g({\bm{{D}}},{\bm{{\theta}}}_{1})\right]$ and we have defined the variance $v_{g}({\bm{{\theta}}}_{0})={\mathrm{Var}}_{{\bm{{\theta}}}_{0}}\!\left[g\!\left({{\bm{{D}}},{\bm{{\theta}}}_{0}}\right)\right]$.

In the following sections we show examples of how a well chosen reparameterized embedding space can result in significant improvement of global maximum testing. Given a statistical model and a candidate embedding $\tilde{{\bm{{\theta}}}}\in\tilde{{\bm{{\Theta}}}}$. Below we specify a computational procedure for selecting a reparameterized embedding space $\tilde{{\bm{{\Theta}}}}$. The procedure is inspired by Rao's locally optimal test cox:1974 ; rao:1948 of the hypotheses $H_{0}:{\bm{{\theta}}}={\bm{{\theta}}}_{0}$ vs. $H_{1}:{\bm{{\theta}}}={\bm{{\theta}}}_{0}+{\bm{{\delta}}}$, where ${\bm{{\delta}}}$ is a small local perturbation of magnitude $\|{\bm{{\delta}}}\|=\alpha$. Rao's score test solves a generalized eigenvalue problem of the form $\max_{\|{\bm{{u}}}\in{\bm{{\Theta}}}\|}{\left|{\bm{{u}}}^{T}{\bm{{s}}}({\bm{{d}}},{\bm{{\theta}}}_{0})\right|^{2}}/{{\bm{{u}}}^{T}{\bm{{I}}}({\bm{{\theta}}}_{0}){\bm{{u}}}}$, where ${\bm{{s}}}({\bm{{d}}},{\bm{{\theta}}})=\nabla_{{\bm{{\theta}}}}\ell({\bm{{d}}};{\bm{{\theta}}})$ is the score function and ${\bm{{I}}}({\bm{{\theta}}})=-E_{{\bm{{\theta}}}}[\nabla_{{\bm{{\theta}}}}^{2}\ell({\bm{{d}}},{\bm{{\theta}}})]$ is the Fisher information matrix. The solution of the generalized eigenvalue problem gives the optimal direction vector ${\bm{{u}}}={\bm{{I}}}^{-1/2}({\bm{{\theta}}}_{0}){\bm{{s}}}({\bm{{d}}},{\bm{{\theta}}}_{0})$ specifying the locally optimal perturbation as ${\bm{{\delta}}}=\alpha{\bm{{u}}}$.

Spectral embedding procedure: Motivated by Theorems 2 and 3 and Rao's test, we propose a heuristic procedure to improve on a randomly initialized embedding space using a singular value decomposition (SVD) on a set of Rao locally optimal direction vectors that are computed offline. The effectiveness of this procedure is demonstrated in the numerical results sections below.

First, an $\epsilon$-net of sampled parameters ${\bm{{\theta}}}_{0}^{(i)}$, $i=1,\ldots,p$, is constructed on ${\bm{{\Theta}}}$. For each ${\bm{{\theta}}}_{0}^{(i)}$ the non-global local maxima $\left\{\widehat{{\bm{{\theta}}}}^{(i,j)}\right\}_{i,j=1}^{p,q_{i}}$ of $m({\bm{{\theta}}}_{0}^{(i)},{\bm{{\theta}}})$ are found. The likelihood function is reparameterized into an initial $\tilde{{\bm{{\Theta}}}}$ space of higher dimensional and the mean $m({\bm{{\theta}}}_{0}^{(i)},\tilde{{\bm{{\theta}}}})=E_{{\bm{{\theta}}}_{0}^{(i)}}[\tilde{\ell}({\bm{{D}}};\tilde{{\bm{{\theta}}}})]$ of the augmented log-likelihood function is computed. Let $\tilde{{\bm{{\theta}}}}^{(i,j)}\in{\bm{{\Theta}}}\times{\bm{{\Theta}}}^{{}^{\prime}}$ be the point in the initial space associated with $\widehat{{\bm{{\theta}}}}^{(i,j)}$ for each $i=1,\ldots,p,\ j=1,\ldots,q_{i}$. Then the solutions $\left\{{\bm{{u}}}(\tilde{{\bm{{\theta}}}}^{(i,j)})\right\}_{i,j}$ of the generalized eigenvalue problem are arranged as the columns of a matrix ${\bm{{G}}}$, whose left principal singular-vector ${\bm{{r}}}$ is taken as a basis for the final spectral embedding subspace in ${\bm{{\Theta}}}^{\prime}$. If more than one additional embedding dimension is desired, this procedure can be repeated with previously identified embedded subspaces successively removed.

3.4 Example: Sinusoidal frequency estimation

We return to the sinusoid in additive Gaussian noise example presented in Section 2.2 to illustrate the theory presented in the previous section and to demonstrate the advantages of the proposed one-sided version of the Biernacki test (16) and the one-sided reparameterized embedding test (32).

We embed the one dimensional frequency parameter $\theta_{0}\in\mathbb{R}$ into the expanded parameter $\tilde{{\bm{{\theta}}}}=(\theta_{0},\theta_{1},\ldots,\theta_{k})\in\mathbb{R}^{k+1}$ and link $\tilde{{\bm{{\theta}}}}$ to the augmented log-likelihood through the augmented mean function $\tilde{{\bm{{\mu}}}}\in\mathbb{R}^{N}$

[TABLE]

This reparameterization embedding introduces variation into the instantaneous frequency of the sinusoid, which enhances the gap between the native and augmented log-likelihood functions at the local maxima, leading to improved detection performance using the test (32). With this embedding, the max gap function $G$ defined in (25) takes the form

[TABLE]

The max mean gap function $m_{G}(\theta_{0},\widehat{\theta})$ (26), which, by Theorem 2, lower bounds the mean shift in $G({\bm{{d}}};\widehat{\theta})$, has the representation:

[TABLE]

where, for large $N$,

[TABLE]

is the maximum of a Fresnel integral. The forms $g$ and $m_{g}$ for the augmented validation function (27) and the ambiguity function (28) are similar except that the minimizations include $\theta_{0}$ and a local minimization is performed.

To simplify the discussion we will respectively reparameterize the augmented and native log-likelihood functions by the associated mean functions: $\tilde{{\bm{{\mu}}}}(\tilde{{\bm{{\theta}}}})$, defined in (33), and ${\bm{{\mu}}}(\theta)=\tilde{{\bm{{\mu}}}}(\theta,{{\bm{{0}}}})$, which lie in different subspaces of $\mathbb{R}^{N}$. With this reparameterization the native log-likelihood $\ell({\bm{{d}}};\theta)$ becomes $\ell({\bm{{d}}};{\bm{{\mu}}})$, with ${\bm{{\mu}}}=\mathrm{sin}(\theta{\bm{{x}}})$ lying in a one dimensional subspace while the augmented log-likelihood $\tilde{\ell}({\bm{{d}}};\tilde{{\bm{{\theta}}}})$ becomes $\tilde{\ell}({\bm{{d}}};\tilde{{\bm{{\mu}}}})$, with $\tilde{{\bm{{\mu}}}}$ lying in a higher dimensional subspace. In particular, the spectral embedding procedure, described at the end of the previous subsection, yields a mean vector $\tilde{{\bm{{\mu}}}}$ lying in the two dimensional subspace $\left\{\tilde{{\bm{{\mu}}}}:\tilde{{\bm{{\mu}}}}=\mathrm{sin}(\theta_{0}{\bm{{x}}})+\theta_{1}{\bm{{r}}},\theta_{0},\theta_{1}\in\mathbb{R}\right\}$, where ${\bm{{r}}}$ is the left principal singular vector of a $N\times\left(\sum_{i=1}^{p}q_{i}\right)$. The obtained vector ${\bm{{r}}}$ is shown in Figure 2 for the case that the number of time samples is $N=100$, the measurement time interval is $T=1$, the Gaussian noise variance is $\sigma^{2}=1$ and the local maxima search region $\theta\in[0,4\pi]$. This particular reparameterized embedding vector ${\bm{{r}}}$ defines the embedded space used in all simulations described in this subsection.

Figure 3 gives a graphical depiction of the proposed reparameterized embedding method in the context of this example. The left panel of the figure represents the case where the hypothesis $H_{0}:\widehat{\theta}=\theta_{Global}$ is rejected, while in the right panel $H_{0}$ is not rejected. The figure shows the likelihood trajectories in the native parameterization ${\bm{{\mu}}}(\theta)$ as a green curve labeled $\mathcal{U}$. The corresponding likelihood surface in the embedded parameterization $\tilde{{\bm{{\mu}}}}(\tilde{{\bm{{\theta}}}})$ is shown as the disk labeled $\mathcal{\tilde{U}}$.

The left and right panels in Figure 3 each show two local maxima, a local maximum proximal to the true parameter $\theta_{0}$, denoted $\widehat{\theta}_{0}$, and a local maximum distant from the true parameter, denoted $\widehat{\theta}$. Associated with these local maxima are the expanded parameters $\tilde{{\bm{{\theta}}}}_{0}$ and $\widehat{\tilde{{\bm{{\theta}}}}}$, respectively, that are depicted in basins of attraction $S\!\left({\widehat{\theta}}\right)$ of the augmented log-likelihood $\max_{\tilde{{\bm{{\theta}}}}}\tilde{\ell}\!\left({{\bm{{d}}},{\bm{{\mu}}}\!\left({\widehat{\theta},\tilde{{\bm{{\theta}}}}}\right)}\right)$. On the right panel the proximal local maximum $\widehat{\theta}_{0}$ is close to the true parameter $\theta_{0}$ so that the means $\mu\!\left({\widehat{\theta}_{0}}\right)$ and $\tilde{\mu}\!\left({\widehat{\tilde{{\bm{{\theta}}}}}_{0}}\right)$ are close to each other. By contrast, on the left panel these two means are not close to each other when the local maximum $\widehat{\theta}$ is distant from $\theta_{0}$. The reparameterized embedding enhances this contrast by increasing the gap between the augmented log-likelihood and the native log-likelihood.

Figure 4 shows the surface corresponding to the (negative) mean augmented log-likelihood, $\left\{-E_{\theta_{0}}[{\ell}({\bm{{d}}};{\bm{{\mu}}}(\tilde{{\bm{{\theta}}}}))]\right\}_{\tilde{{\bm{{\theta}}}}\in\tilde{\Theta}}$ as a surface over $\tilde{\theta}=(\theta,\theta^{{}^{\prime}})$ along with its cross-section (blue curve) along the line $\left\{(\theta,0):\theta\in\Theta\right\}$, which is the (negative) native mean log-likelihood trajectory $\left\{-E_{\theta_{0}}[\ell({\bm{{d}}},{\bm{{\mu}}}(\theta))]\right\}_{\theta\in\Theta}$. The mean augmented log-likelihood has two local maxima, a global maximum $\widehat{\tilde{{\bm{{\theta}}}}}_{0}=(\widehat{\theta}_{0},\widehat{\theta}^{\prime}_{0})$ near $(\theta_{0},0)$ and a local maximum $\widehat{\tilde{{\bm{{\theta}}}}}=(\widehat{\theta},\widehat{\theta}^{\prime})$. The green curve traces out $-\max_{\theta^{\prime}\in\Theta^{\prime}}E_{\theta_{0}}[\ell({\bm{{d}}},\tilde{{\bm{{\mu}}}}(\theta,\theta^{\prime}))]$ and the gap between the blue and green curves is the augmented ambiguity function $m_{g}(\theta_{0},\theta)$ (28).

Figure 5 demonstrates the improvements achieved using the proposed reparameterized embedding test (32) with respect to testing accuracy, as quantified by receiver operating characteristic (ROC) curves. The ROC curve sweeps out the probability of detection (PD) against the probability of false alarm (PFA) achieved by a test that a local maximum is the global maximum. Tests with higher ROC curves are more accurate. 10,000 simulations were performed with the following parameters: a true sinusoidal frequency $\theta_{0}=3\pi$ and, $N=100$, $T=1$, $\sigma^{2}=1$, and search region $\theta\in[0,4\pi]$. In order of increasing accuracy these tests are: the two-sided Biernacki test biernacki:2005 , the one-sided version of the Biernacki test (16), equivalent to the proposed test with $k=0$, followed by the proposed one-sided with the embedding given by (33) for $k=1$ and $k=3$ respectively. Using $k=3$ or the one-sided test based on a single embedding basis chosen according to the proposed spectral approach, achieves significantly better performance than Biernacki tests. For all of the reparameterized embedding tests the maximization in the augmented validation function (34) was computed numerically using a limited-memory quasi-Newton solver nocedal:1999 .

4 Application to Wavefront Sensing

We apply the proposed global maximum testing framework to the problem of jointly estimating camera blur and pose from a known calibration target in the presence of aliasing, studied in leblanc:2018 . Wavefront phase-aberrations characterize the camera's point spread function (PSF) through deviations of an otherwise ideal optical system, and a Zernike polynomial basis Zernike:1934 is used to parameterize these phase aberrations resulting in a parametric blur model. The solution to the inverse problem is the global maximum of a highly nonconvex log-likelihood function, that admits many local maxima barakat:1992 ; isernia:1995 ; gerwe:2008 ; moretta:2019 . The inverse problem is typically solved using iterative optimization algorithms that converge to local maxima. Thus, in this context, a test for global maximum is a test of global convergence of the iterative algorithm. We demonstrate how the proposed test of global convergence can be used to reduce the computation burden.

The spectral embedding procedure described in Section 3.3 is used in conjunction with the one-sided test (32) to substantially improve the detectability of convergence of the estimated blur to a suboptimal maximum of the log-likelihood function. Under the generalized imaging model goodman:1996 , the PSF $h$ is described in terms of a non-negative aperture function $A$ and a real-valued phase function $\Psi$, expressed as a linear combination of Zernike polynomials, where the parameter vector ${\bm{{\theta}}}$ contain the basis coordinates. Specifically, $h$ is given as

[TABLE]

where $c_{0}$ is a normalizing constant that ensures the PSF integrates to 1, $g\!\left({w_{x},w_{y};{\bm{{\theta}}}}\right)$ is the coherent transfer function (CTF), and $\mathcal{F}^{-1}$ is the inverse Fourier transform. A test for global convergence is constructed using the proposed spectral method to embed this parametric blur model into the space of non-negative PSFs on a Nyquist-sampled grid. Monte Carlo simulations performed at a moderate signal-to-noise ratio (20 dB) and blur strength (0.025 waves RMS) demonstrate the substantial power of the proposed test. Over 100 such trials, a limited-memory quasi-Newton search started from the point representing no phase aberrations (an ideal imaging system) led to non-global local maxima 96% of the time. Table 1 compares the power of the Biernacki's test biernacki:2005 with the proposed approach when both tests are operated at a false alarm rate of 0.01. The observed improvement is consistent with the simulation example described in Section 3.4.

Given a local maximum of the log-likelihood suspected of being a suboptimal solution, one would like to exploit knowledge of this local maximum to identify alternative regions of the parameter space likely to contain a better solution. The PSF $\tilde{h}$ corresponding to a perturbed phase-screen $\Psi+\beta$ can be expressed as the modulus squared of a convolution of CTFs associated with $\Psi$ an $\beta$ respectively

[TABLE]

where $A_{B}$ is the binary aperture corresponding to the support of $A$. Letting $\left[h\right]_{m,n}$ be the $(m,n)^{\textrm{th}}$ element of a Nyquist sampled representation of the PSF, then a point-wise bound on the magnitude of the change induced by $\beta$ is given by

[TABLE]

where $\delta$ is the Kronecker delta function, and $a$ an arbitrary complex constant such that ${\left|a\right|}=1$. This point-wise bound on the PSF perturbation $\epsilon$, associated with the wavefront perturbation $\beta$, is minimized when $\angle a=\angle\left[\tilde{g}\right]_{0,0}$. Under this condition, the right-hand side of (35) is monotonic in the Strehl ratio martial:1991 associated with $\beta$, which we will denote as $c_{0}{\left|\left[\tilde{g}\right]_{0,0}\!\left({\beta}\right)\right|}^{2}$. Thus, the set of wavefronts that maximize the Strehl ratio for a fixed root mean square (RMS) perturbation strength also minimizes the worst-case, point-wise error in the perturbed PSF. These wavefronts are given by

[TABLE]

A PSF $h$ perturbed by a wavefront in (36) will result in a new PSF $\tilde{h}$ that is point-wise close to $h$ despite its wavefront $\Psi+\beta$ being $\tau$ waves RMS from $\Psi$. The proposed test for global convergence can be used in conjunction with this restarting strategy to search for globally optimal solutions.

A Monte Carlo simulation was performed to assess the efficacy of the proposed approach for identifying the global maximum of the likelihood function. A simulated imaging system leblanc:2018 was configured to provide moderate SNR images (20 dB), and the number of blur aberration parameters was varied to alter the difficulty of the resulting inverse problem. As the number of Zernike modes in the model increases, so does the probability of encountering local maxima. A limited-memory quasi-Newton search nocedal:1999 was used to identify stationary points of the log-likelihood starting from a diffraction-limited model. If the reparameterized embedding approach described in Section 3.3 failed to reject the null hypothesis at a false alarm level $\alpha=0.01$ the search was terminated, otherwise a new starting point 0.2 waves RMS away from the current maximum was chosen according to (36), and the search was continued. Figure 6 shows the mean and standard errors of runtimes corresponding to 10 independent realizations of the same camera model. For comparison, the simulated annealing algorithm provided in $\textrm{MATLAB}^{\circledR}$ Optimization Toolbox version 8.0 was used as a point of reference. The simulated annealing algorithm was provided the objective function gradients and was terminated according to an oracle criterion: terminate the first time that any local maximum fell within 0.01 waves RMS of the true solution. In the astronomical imaging community, wavefront descriptions of optical systems typically include Zernike models up to at least radial-order 3 (7 Zernike modes). For models of such high complexity, the proposed reparameterized embedding strategy resulted in a five times reduction in total runtime. Figure 7 illustrates a typical sequence of PSFs associated with the global search procedure when 12 Zernike modes parameterize the blur. Despite the relatively small differences between the PSFs, the two non-global local maxima (local #1 and local #2) of the log-likelihood are associated with relatively large wavefront perturbations errors of 0.117 and 0.115 waves RMS, respectively.

5 Concluding Remarks

This paper addresses a principal computational bottleneck in non-convex imaging and vision problems: It determines if a local maximum found by a non-global optimization algorithm is a global maximum. Specifically, we introduced a powerful new method for validating that a local maximum is a global maximum when the objective function is specified as the likelihood function associated with a parametric statistical model. The proposed method implements a one-sided threshold test on a novel validation function defined as the normalized difference between the log-likelihood function and an augmented log-likelihood function, each evaluated at a local maximum point. The augmented log-likelihood is constructed by embedding the original parameter vector into a higher dimensional parameter space, a procedure we call reparameterized embedding, and the validation function is evaluated at the local maximum before thresholding. We proposed a computational spectral embedding procedure for identifying good reparameterized embeddings, and numerical results are presented exhibiting an extraordinarily high level of detection accuracy, e.g., achieving significantly better accuracy than the two-sided test proposed by Biernacki. Finally, to demonstrate how our results can dramatically impact non-convex imaging applications, we applied the proposed test to a set of local maxima generated from multiple restarts of an iterative maximum likelihood algorithm for reconstructing camera blur from images of a calibration target. When the test is used as a stopping rule, i.e., the restarts are stopped when the test declares a global maximum has been found, it reduced total runtime by a factor of five.

Code for reproducing the key figures from this document is available at https://github.com/jwleblan/localMinima.

Acknowledgements.

This work was partially supported by ARO grant W911NF-15-1-0479 and a DOE NNSA grant to the University of Michigan Consortium on Verification Technology.