Maximum Likelihood Estimation based on Random Subspace EDA: Application to Extrasolar Planet Detection

TL;DR

This paper introduces a novel subspace-based estimation of distribution algorithm to improve maximum likelihood estimation for complex, high-dimensional models in extrasolar planet detection, addressing the curse of dimensionality.

Contribution

It proposes a new interactive random subspace mechanism for EDAs, enhancing scalability in high-dimensional, nonlinear likelihood surface exploration.

Findings

Effective in numerical benchmarks

Successful application to real data

Reduces dimensionality challenges in ML estimation

Abstract

This paper addresses maximum likelihood (ML) estimation based model fitting in the context of extrasolar planet detection. This problem is featured by the following properties: 1) the candidate models under consideration are highly nonlinear; 2) the likelihood surface has a huge number of peaks; 3) the parameter space ranges in size from a few to dozens of dimensions. These properties make the ML search a very challenging problem, as it lacks any analytical or gradient based searching solution to explore the parameter space. A population based searching method, called estimation of distribution algorithm (EDA), is adopted to explore the model parameter space starting from a batch of random locations. EDA is featured by its ability to reveal and utilize problem structures. This property is desirable for characterizing the detections. However, it is well recognized that EDAs can not scale well to large scale problems, as it consists of iterative random sampling and model fitting procedures, which results in the well-known dilemma curse of dimensionality. A novel mechanism to perform EDAs in interactive random subspaces spanned by correlated variables is proposed and the hope is to alleviate the curse of dimensionality for EDAs by performing the operations of sampling and model fitting in lower dimensional subspaces. The effectiveness of the proposed algorithm is verified via both benchmark numerical studies and real data analysis.

Tables

Table 1. Table 1: Parameter value ranges of RV models

$P_{\min }$	$1$ day	$P_{\max }$	$1,000$ years
$K_{\min }$	$1$ m/s	$K_{\max }$	$2128$ m/s
$C_{min}$	$-2128$ m/s	$C_{max}$	$2128$ m/s
$s_{\min }$	$1$ m/s	$s_{\max }$	$2128$ m/s

Table 2. Table 2: The calculated BIC metrics of ℳ 1 subscript ℳ 1 {\cal{M}}_{1} and ℳ 2 subscript ℳ 2 {\cal{M}}_{2}

	$\ln \left(\text{ML}\right)$	BIC
$\text{One Planet Model}:{ℳ}_1$	-148.4024	320.6132
$\text{Two Planet Model}:{ℳ}_2$	-111.5458	263.9060

Equations

$$v_i\mid {\mathcal{M}}_0\sim \mathcal{N}\left(C,{\sigma }_i^2+s^2\right),$$

$$v_i\mid {\mathcal{M}}_1\sim \mathcal{N}\left(C+△V(t_i|{\varphi }_1),{\sigma }_i^2+s^2\right),$$

$$△V(t|\varphi )=K[\cos (\omega +T(t))+e\cos (\omega )],$$

$$T(t)=2\arctan \left[\tan \left(\frac{E(t)}{2}\right)\sqrt{\frac{1+e}{1-e}}\right],$$

$$E(t)-e\sin (E(t))=\text{mbox}mod\left(\frac{2\pi }{P}t+{\mu }_0,2\pi \right).$$

$$△V(t_i|{\varphi }_1,\dots ,{\varphi }_j)=\sum _{a=1}^j△V(t_i|{\varphi }_a).$$

$$\text{mbox}{BIC}_i=-2\cdot \ln {\hat{L}}_i+k_i\cdot \ln (n)$$

$${\hat{\theta }}_i=\underset{\theta \in {\Theta }_i}{\mathrm{argmax}}p(\mathbf{v}|\theta ,{\mathcal{M}}_i).$$

$$f(\theta )={\left(\frac{1}{d-1}\sum _{i=1}^{d-1}\left(\sqrt{z_i}+\sqrt{z_i}{\sin }^2\left(50z_i^{0.2}\right)\right)\right)}^2-800$$

Full text

1 Introduction

This paper presents an evolutionary computation based maximum likelihood (ML) estimation method for multivariate highly nonlinear time series models. This work is motivated by a challenging signal detection task, which aims to detect extrasolar planets (exoplanets) based on observations collected by astronomical instruments such as NASA’s Kepler space telescope [1]. The terminology “exoplanets” denotes planets outside our solar system. The goal of exoplanet science is to answer the scientific quest whether we are alone or whether there are other planets that might support life in the universe [2, 1, 3]. Exoplanet science has become a booming field in astrophysics since 1992 when the first detection of an exoplanet was confirmed [4]. In this paper, we focus on the radial velocity (RV) method of exoplanet detection [5, 6, 3]. This method has become one of the most productive techniques for detecting exoplanets so far.

Signal processing (SP) plays an important part in exoplanet detection, which contributes to improve the signal-to-noise ratio of the observations, detect signals of potential planets and so forth. Current SP techniques fall short of meeting the fundamental requirement for the future of this field. For example, the ML based periodogram method was developed to deal with correlated noise in RV time series [7], while, it is just limited to detect one signal. Many planetary systems are found to contain more than one planet, which means the RV time series should exhibit more than one signal. The Bayesian simulation techniques have been applied to explore the parameter space of a global RV model using Markov Chain Monte Carlo (MCMC) or adaptive importance sampling methods [3, 8, 6], while such methods require very large computational overhead to guarantee a satisfactory performance in parameter estimation and signal detection.

The objective of this paper is to propose a computationally efficient method to address the problem of exoplanet detection based on ML fitting of complex RV models. This problem is featured by the following properties: 1) the candidate models under consideration are highly nonlinear; 2) the likelihood surface has a huge number of peaks; 3) the parameter space ranges in size from a few to dozens of dimensions. This problem lacks any analytical or gradient based searching solution to explore the parameter space. The proposed approach is based on an evolutionary computation method called estimation of distribution algorithm (EDA) [9, 10, 11]. A novel mechanism is proposed to perform EDAs in a series of random subspaces spanned by correlated variables. The basic idea is that, since the dimension of each subspace will be smaller than that of the full parameter space, the number of occurrences of sample-starved model fitting will decrease and then the curse of dimensionality is hoped to be alleviated. A benchmark numerical study and a real data analysis are used to demonstrate the effectiveness of the new algorithm.

2 RV Models and the ML based Model Fitting

In this section, we present the RV time series models and then introduce the ML parameter fitting problem in the context of exoplanet detection.

2.1 RV Models

A succinct introduction of the RV models is presented here. For more details, readers are referred to [6]. We use ${\cal{M}}_{j}$, $j=0,1,\ldots,J$, to denote the $j$-planet model, corresponding to the hypothesis that there is (are) $j$ planet(s) in the extrasolar system under consideration. In the 0-planet model ${\cal{M}}_{0}$, the $i$th element of the RV data $v_{i}$ is modeled to be Gaussian distributed as follows

[TABLE]

where $C$ and $\sigma^{2}_{i}+s^{2}$ are its mean and variance, respectively. Here $C$ denotes constant center-of-mass velocity of the star relative to earth and $s$ denotes the square root of the “stellar jitter”, which represents the random fluctuations in a star’s luminosity or fluctuations stemming from other systematic sources, e.g., starspots. The additional variance component $\sigma^{2}_{i}$ is a calculated error of $v_{i}$ due to the observation procedure.

In the 1-planet model ${\cal{M}}_{1}$, $v_{i}$ is modeled as follows

[TABLE]

where $\bigtriangleup V(t_{i}|\phi)$ is the velocity shift caused by the presence of the planet. Such velocity shift is a family of curves parameterized by a 5-dimensional vector $\phi\triangleq(K,P,e,\omega,\mu_{0})$ defined as follows

[TABLE]

where $T(t)$ is the “true anomaly at time $t$” given by

[TABLE]

and $E(t)$ is called the “eccentric anomaly at time $t$”, which is the solution to a transcendental equation

[TABLE]

In the above expressions, $K$ denotes the velocity semi-amplitude, $P$ the orbital period, $e$ the eccentricity $(0\leq e\leq 1)$, $\omega$ the argument of periastron $(0\leq\omega\leq 2\pi)$, and $\mu_{0}$ the mean anomaly at time $t=0$, $(0\leq\mu_{0}\leq 2\pi)$. The parameters $C$, $K$ and $s$ have the same unit as velocity; the velocity semi-amplitude $K$ is usually restricted to be non-negative to avoid identification problems; $C$ may be positive or negative. The eccentricity parameter $e$, $0\leq e<1$, is unitless, with $e=0$ corresponding to a circular orbit, and larger $e$ means more eccentric orbits. Periastron is the point at which the planet is closest to the star and the argument of periastron $\omega$ measures the angle at which we observe the elliptical orbit. The mean anomaly $\mu_{0}$ is an angular distance of a planet from periastron.

More generally, a $j$-planet model ($j\geq 1$) represents the expected velocity by $C+\bigtriangleup V(t_{i}|\phi_{1},\ldots,\phi_{j})$, in which the overall velocity shift $\bigtriangleup V$ takes the form of the summation of velocity shifts of each individual planet, i.e.,

[TABLE]

Therefore the parameter dimension of a $j$-planet model is $2+5j$. So the more planets covered by the model, the higher dimensional it is. Moreover, the RV models ${\cal{M}}_{1},{\cal{M}}_{2},\ldots$ are highly nonlinear due to the velocity shift item $\bigtriangleup V$ included in these models.

2.2 ML Parameter Estimation of RV Models

Here we treat exoplanet detection as a problem of model selection, which means, given a set of RV observations $\mathbf{v}\triangleq(v_{1},\ldots,v_{n})$, how to select from the candidate models $\{{\cal{M}}_{0},{\cal{M}}_{1},\ldots,{\cal{M}}_{j}\}$ the one that fits the data best in terms of Bayesian criterion. A full Bayesian solution needs to calculate the marginal likelihood of each candidate model, which involves large scale stochastic integrations over the whole parameter space of the candidate models [3, 6]. It is computationally expensive to solve such stochastic integrations. Here we resort to the Bayesian information criterion (BIC) to evaluate the fitness of the candidate models to the data. A Bayesian argument for adopting BIC was presented in [12]. We use $\theta_{i}$ and ${\Theta_{i}}$ to denote the parameters of ${\cal{M}}_{i}$ and the corresponding value space, respectively. The BIC metric of ${\cal{M}}_{i}$ is defined as

[TABLE]

where $k_{i}$ is the number of free parameters to be estimated for ${\cal{M}}_{i}$, $\hat{L}_{i}$ is the maximal likelihood function value associated with ${\cal{M}}_{i}$, i.e., $\hat{L}_{i}=p(\mathbf{v}|\hat{\theta}_{i},{\cal{M}}_{i})$, where $\hat{\theta}_{i}$ are the parameter values that maximize the likelihood function, namely

[TABLE]

Given a finite set of models, the model with the lowest BIC value is preferred, according to the BIC criterion. Since the second item on the right-hand side of Eqn. (7) is a constant given the model, calculation of $\mbox{BIC}_{i}$ then translates to how to solve the maximization problem defined in Eqn. (8).

3 General EDA Procedure

To address an optimization problem such as that defined in Eqn. (8), a general EDA procedure typically works with a population of candidate solutions defined over the full parameter space. The initial population is generated according to the uniform distribution over all admissible solutions. The fitness function gives a numerical ranking for each solution. Here the likelihood function $p(\mathbf{v}|\theta,{\cal{M}}_{i})$ plays the role of a fitness function. A subset of the most promising solutions is selected by the selection operator. A commonly used selection operator selects a certain proportion, e.g., the best 50% of solutions. A probabilistic model is then constructed to estimate the probability distribution of the selected solutions. Given the above model, the algorithm generates new solutions by sampling the distribution defined by the model. The new population of solutions replaces the old population and then the modeling and sampling procedure is repeated until some termination criteria are met. The main scheme for an iteration of the EDA method is summarized as follows: starting from a population of solutions $P$,

• •

Select a population of promising solutions $S$ from $P$;

• •

Build a probabilistic model $M$ from $S$;

• •

Sample new candidate solutions $Q$ based on $M$;

• •

Replace the old population with the new population, namely set $P$ to be $Q$.

For more details about the EDA algorithm, readers are referred to [11]. As an iterative sampling and modeling procedure, the general EDA method suffers from the well-known dilemma curse of dimensionality. Specifically speaking, with the increase in the dimension of the solution space, the volume of the space increases so fast that the available data points used for constructing the model become sparse, making the resulting model not qualified for guiding the searching process to find better solutions.

4 Random Subspace EDA (RS-EDA)

In this section, we propose a new EDA algorithm, namely RS-EDA. The idea is to partition the original multivariate parameter space into a series of random subspaces, and then perform EDAs in these subspaces, rather than the full parameter space, with the hope to alleviate the curse of dimensionality for EDA type methods. Fig.1 shows one iteration of the RS-EDA method, wherein $d$ denotes the dimension of $\theta$. The iteration ends when the estimate of the global optimal solution keeps unchanged for a fixed number of iterations. Taking the maximization problem as an instance, we describe in what follows the details of the four operators included in Fig.1.

4.1 The 1st Step: Estimation of Variable Correlations

This step corresponds to the leftmost box in Fig.1. The purpose is to provide a rough estimate on correlations among the variables (or dimensions) of $\theta$. Such estimation is made based on a pool of promising solutions (PPS), given by the previous one iteration of the RS-EDA. If it is currently in the first iteration, we initialize the PPS based on available prior knowledge on the variable correlations. If there is no such prior knowledge, we just draw a number of random samples uniformly from the solution space and then initialize the PPS via an operation of truncation selection. Specifically, we preserve the fittest 20% of the sample for use and throw away the others. The data in PPS is formulated by an $a\times d$ matrix, where $a$ denotes the number of samples in PPS. The operation of correlation estimation returns the sample linear partial correlation coefficients between pairs of variables in PPS, controlling for the remaining variables in PPS. We use MATLAB’s built-in function “partialcorr” to perform the above operation. The output of this function is a symmetric $d\times d$ matrix $C$, whose $(i,j)$-th entry is the sample linear partial correlation corresponding to the $i$th and $j$th columns of the PPS matrix.

4.2 The 2nd Step: Random Subspace Partition

This step corresponds to the second box in Fig.1. The purpose is to partition the whole parameter space into subspaces based on the correlation matrix $C$ given by the 1st step. This 2nd step consists of a series of procedures performed on each row of $C$. Take the procedures corresponding to $C_{i}$, the $i$th row of $C$, for example. We first sort the elements of $C_{i}$, namely $\{C_{i,1},\ldots,C_{i,d}\}$, according to their values from large to small. This procedure outputs $\{C_{i,j1},\ldots,C_{i,jd}\}$, where $\{j1,\ldots,jd\}$ is a rearrangement of $\{1,\ldots,d\}$ with $C_{i,jm}\geq C_{i,jn}$ as long as $m<n$. Let $S=\sum_{m=1}^{d}C_{i,jm}$ and then set $C_{i,jm}=C_{i,jm}/S$ for each $m$ in $\{1,\ldots,d\}$. Now we have $\sum_{m=1}^{d}C_{i,jm}=1$. Then we set $C^{\prime}_{i,m}=\sum_{k=1}^{m}C_{i,jk}$ for $m=1,\ldots,d$, draw a random number $r$ from a uniform distribution between 0 and 1, find the minimum index $m$ from $\{1,\ldots,d\}$ that satisfies $C^{\prime}_{i,m}\geq r$, and finally return the indexes $\{j1,\ldots,jm\}$, which we use to constitute the coordinates of the $i$th subspace. After traversing each row of $C$, we build up a series of overlapped random subspaces.

The basic idea underlying the above operations is that the more correlation between a pair of variables, the more likely they will be incorporated into the same lower dimensional subspace. Employing the random mechanism presented above, we do not need to introduce any free parameter, such as a threshold, to cluster the variables into different subspaces.

4.3 The 3rd Step: Performing EDAs in Subspaces

This step corresponds to the elliptical box in Fig.1. Now we focus on a specific subspace and present the EDA operations endowed with it. Assume that the current estimate of the global optimal solution takes value at $\hat{\theta}=(\hat{\theta}_{1},\ldots,\hat{\theta}_{d})$ and the subspace under consideration is associated with variables $(\theta_{j1},\ldots,\theta_{jm})$. We use a Gaussian model to fit the PPS data mapped to this subspace. Specifically, we calculate the empirical mean and sample covariance of this Gaussian model based on the PPS data mapped to the dimensions $\{j1,\ldots,jm\}$. Then we draw $R\times m$ new samples from this Gaussian distribution, where $R$ is a constant prescribed beforehand. As the dimensions of all the exoplanet models of our concern are less than 100, we select $R=100$ in our experiments, which is big enough to prevent from sample-starved model fitting in the follow-up EDA operations. We assign these newly generated samples’ values to the corresponding variables $\{j1,\ldots,jm\}$ of $\hat{\theta}$ and keep the other variables unchanged. Then we get a set of full dimensional samples. We calculate the fitness values of these samples, based on which we select 20% fittest samples for use in updating the Gaussian model. During the above process, once better solutions are found, we shall update the estimate of the global optimal solution and the PPS accordingly to guarantee that the PPS maintains the fittest samples that have been found so far. Then we iterate the model fitting, sampling and selection procedure until the estimate of the global optimal solution keeps unchanged in the most recent continuous five iterations. We regard this phenomenon as an indication of algorithm convergence.

4.4 The 4th Step: Updating the estimate of the global optimum and the PPS

This step occupies the rightmost box in Fig.1 for ease of presentation, while, in practice, it is totally interactive with the 3rd step presented above. Once better solutions have been found in an EDA procedure included in the 3rd step, the estimate of the global optimal solution and the PPS will be updated accordingly. In this way, the PPS always keeps a set of the fittest solutions. The EDA procedure of the next subspace will be performed based on the most recent estimate of the global optimal solution and the most recently updated PPS. After carrying out the EDA procedures of all subspaces, the finally outputted PPS will then act as the input for the 1st step in the next iteration.

5 Numerical Study with Benchmark Function

To test the potential of our idea and the ability of our algorithm to improve the scalability of EDAs to search a near-optimal solution in large-scale problem settings, we tested it based on a benchmark function listed in the suite of benchmark test functions released by a special session on real-parameter optimization of the 2013 Congress on Evolutionary Computation [13]. We selected the Rotated Schaffers F7 (RSF7) function, which is multimodal, nonseparable, asymmetrical, and has a huge number of the local optimum. It is defined as follows [13]

[TABLE]

where $z_{i}=\sqrt{y_{i}^{2}+y_{i+1}^{2}}$ for $i=1,\ldots,d$, $y=\Lambda_{d}^{10}\mbox{M}_{2}T_{asy}^{0.5}(\mbox{M}_{1}(\theta-o))$. In the above expressions, $o$ is the shifted global optimum, $\Lambda_{d}^{\alpha}$ denotes a $d$ dimensional diagonal matrix, the $i$th diagonal element of which is $\alpha^{\frac{i-1}{2(d-1)}}$, $i=1,2,\ldots,d$. $T_{asy}^{\beta}(x)$ is an operator that transforms $x_{i}$ to be $x_{i}^{1+\beta\frac{i-1}{d-1}\sqrt{x_{i}}}$ for $x_{i}>0$, $i=1,2,\ldots,d$. $\mbox{M}_{1}$ and $\mbox{M}_{2}$ are orthogonal matrices whose entries are standard normally distributed. More details about this function can be found in [13].

The RSF7 function is designed for minimization tasks. We consider maximizing the function $g(\theta)=-f(\theta)$ in order to mimic the ML estimation problem. The maximum function value of $g(\theta)$ is 800, which is the target for the algorithm to search. The first instance we considered was a low dimensional case, in which $d$ was set at 5. We initialized the PPS of RS-EDA using $N=20000$ random samples and we allowed a fixed budget of $10^{6}$ function evaluations. We compared RS-EDA with other three leading evolutionary computation methods, namely the Trelea type vectorized Particle Swarm Optimization (PSO) algorithm [14], the Covariance Matrix Adaptation Evolution Strategy (CMAES) method [15] and a classical EDA method termed EMNA${}_{\mbox{global}}$ [16]. The difference between EMNA${}_{\mbox{global}}$ and our method is that the former always employs a full-dimensional multivariate Gaussian distribution model, while the latter performs EDA operations in subspaces. The population size of EMNA${}_{\mbox{global}}$ is set at $N=20000$. The population size of the PSO is set at $1000$. The CMAES method of [15] was included in our comparison because it is developed to handle high dimensional problems and it currently represents the gold-standard for comparisons in new EDA research. We used the Matlab implementation available from the authors with the diagonal option, default parameter settings and random initialization [17]. We ran all these baseline methods with a maximum number of, i.e., $10^{6}$, function evaluations, the same as for the RS-EDA method. We ran each method for 50 times and calculated the mean of its estimate on the maximal function value. Fig. 2 gives a visual summary of the results obtained in comparison. It is shown that the proposed RS-EDA method finds the global optimum with a slower convergence speed than PSO and CMAES and a faster convergence speed than EDA.

We then focused on a higher dimensional instance with $d$ set at 20. Compared with the first instance, the population size of each method involved increased by 10 times. The maximum function evaluations allowed for each method is $5\times 10^{6}$. Fig. 3 shows the summary of the results corresponding to 50 independent runs of each method. The result of the EDA method totally diverged, so it is not included in Fig. 3. We see that the RS-EDA beats all the other methods. Combining the two instances for a joint analysis, we see a potential of the proposed idea of random subspace in improving the scalability of EDAs in dealing with higher dimensional problems. However, it is worthy of further investigation.

6 Analysis of Real Exoplanet Data

We applied the proposed RS-EDA method to analyze a real data set released in [18]. This data set was claimed by the astronomers to have two planets [18]. An adaptive annealed importance sampling (AAIS) method was developed in [6], which calculates out the probabilities of the hypothetical models based on the given RV measurements.

We fit the data based on ${\cal{M}}_{1}$ and ${\cal{M}}_{2}$, respectively. For each model, we apply RS-EDA to estimate the ML model parameters. The search space is constrained by the value ranges of the model parameters as listed in Table 1. These ranges are provided by the astronomers based on their experiences [6]. Fig. 4 shows the fitting results based on ${\cal{M}}_{1}$ and ${\cal{M}}_{2}$, respectively. The calculated logarithm of MLs and BIC metrics associated with ${\cal{M}}_{1}$ and ${\cal{M}}_{2}$ are listed in Table 2. Both the visual fitting result and the quantitative comparison of the BIC metrics suggest that ${\cal{M}}_{2}$ holds. This result is consistent with that reported in [18] and [6].

7 Concluding remarks

In this paper, we proposed an RS-EDA algorithm in the context of exoplanet detection based on ML estimation. The most important feature of the RS-EDA method lies in the proposed operation of constructing random subspaces and then endowing the routine sampling and model fitting operations of EDAs into the lower dimensional parameter spaces, with the hope to alleviate the curse of dimensionality for EDA type methods. The effectiveness of the proposed method was demonstrated by numerical studies and real data analysis. The results show a potential of the proposed technique in dealing with complex nonlinear models with multimodal likelihood functions and in solving high dimensional optimization tasks. The scalability of the proposed method and new approaches for constructing random subspaces are both worthy of further investigations.

Acknowledgement

This work was partly supported by the National Natural Science Foundation (NSF) of China under grant No. 61571238, China Postdoctoral Science Foundation under grant Nos. 2015M580455 and 2016T90483, the Six Talents Peak Foundation of Jiangsu Province under grant No. XYDXXJS-CXTD-006 and the Scientific and Technological Support Project (Society) of Jiangsu Province under grant No. BE2016776.