Robust and Scalable Bayes via a Median of Subset Posterior Measures

Introduction

Contemporary data analysis problems pose several general challenges. One is resource limitations: massive data require computer clusters for storage and processing. Another problem occurs when data are severely contaminated by “outliers” that are not easily identified and removed. Following Box and Tiao (1968), an outlier can be defined as “being an observation which is suspected to be partially or wholly irrelevant because it is not generated by the stochastic model assumed.” While the topic of robust estimation has occupied an important place in the statistical literature for several decades and significant progress has been made in the theory of point estimation, robust Bayesian methods are not sufficiently well-understood.

Our main goal is to make a step towards solving these problems, proposing a general Bayesian approach that is

provably robust to the presence of outliers in the data without any specific assumptions on their distribution or reliance on preprocessing;

scalable to big data sets through allowing computational algorithms to be implemented in parallel for different data subsets prior to an efficient aggregation step.

The proposed approach consists in splitting the sample into disjoint parts, implementing Markov chain Monte Carlo (MCMC) or another posterior sampling method to obtain draws from each “subset posterior” in parallel, and then using these draws to obtain weighted samples from the median posterior (or M-Posterior), a new probability measure which is a (properly defined) median of a collection of subset posterior distributions. We show that, despite the loss of “interactions” among the data in different groups, the final result still admits strong guarantees; moreover, splitting the data gives certain advantages in terms of robustness to outliers.

In particular, we demonstrate that the M-posterior is a probability measure centered at the “robust” estimator of the unknown parameter, the associated credible sets are often of the same “width” as the credible sets obtained from the usual posterior distribution and admit strong “frequentist” coverage guarantees (see section 3.3 for exact statements).

The paper is organized as follows: section 1.1 contains an overview of the existing literature and explains the goals that we aim to achieve in this work. Section 2 introduces the mathematical background and key facts used throughout the paper. Section 3 describes the main theoretical results for the median posterior. Section 4 presents details of algorithms, implementation, and numerical performance of the median posterior for several models. The simulation study and analysis of data examples convincingly show the robustness properties of the median posterior. We have also implemented the matrix completion example based on the MovieLens data set (GroupLens Research, 2013) which illustrates the scalability of the model. Proofs that are omitted in the main text are contained in the appendix.

A. Dasgupta remarks that (see the discussion following Berger (1994)): “Exactly what constitutes a study of Bayesian robustness is of course impossible to define.” The popular definition (which also indicates the main directions of research in this area) is due to J. Berger (Berger, 1994): “Robust Bayesian analysis is the study of the sensitivity of Bayesian answers to uncertain inputs. These uncertain inputs are typically the model, prior distribution, or utility function, or some combination thereof.” Outliers are typically accommodated by either employing heavy-tailed likelihoods (e.g., Svensen and Bishop (2005)) or by attempting to identify and remove them as a first step (as in Box and Tiao (1968) or Bayarri and Berger (1994)). The usual assumption in the Bayesian literature is that the distribution of the outliers can be modeled (e.g., using a $t$-distribution, contamination by a larger variance parametric distribution, etc). In this paper, we instead bypass the need to place a model on the outliers and do not require their removal prior to analysis. We base inference on the median posterior, whose robustness can be formally and precisely quantified in terms of concentration properties around the true delta measure under the potential influence of outliers and contaminations of arbitrary nature.

Also relevant is the recent progress in scalable Bayesian algorithms. Most methods designed for distributed computing share a common feature: they efficiently use the data subset available to a single machine and combine the “local” results for “global” learning, while minimizing communication among cluster machines (Smola and Narayanamurthy (2010)). A wide variety of optimization-based approaches are available for distributed learning (Boyd et al., 2011); however, the number of similar Bayesian methods is limited. One of the reasons for this limitation is related to Markov chain Monte Carlo (MCMC), the dominating approach for approximating the posterior distribution of parameters in Bayesian models. While there are many efficient MCMC techniques for sampling from posterior distributions based on small subsets of the data (called “subset posteriors” in the sequel), to the best of our knowledge, there is no general rigorously justified approach for combining the subset posteriors into a single distribution for improved performance.

Three major approaches exist for scalable Bayesian learning in a distributed setting. The first approach independently evaluates the likelihood for each data subset across multiple machines and returns the likelihoods to a “master” machine, where they are appropriately combined with the prior using conditional independence assumptions of the probabilistic model. These two steps are repeated at every MCMC iteration (see Smola and Narayanamurthy (2010); Agarwal and Duchi (2012)). This approach is problem-specific and involves extensive communication among machines. The second approach uses a so-called stochastic approximation (SA) and successively learns “noisy” approximations to the full posterior distribution using data in small mini-batches. The accuracy of SA increases as it uses more data. A group of methods based on this approach uses sampling-based techniques to explore the posterior distribution through modified Hamiltonian or Langevin dynamics (e.g., Welling and Teh (2011); Ahn, Korattikara and Welling (2012); Korattikara, Chen and Welling (2013)). Unfortunately, these methods fail to accommodate discrete-valued parameters and multimodality. Another subgroup of methods uses deterministic variational approximations and learns the variational parameters of the approximated posterior through an optimization-based approach (see Wang, Paisley and Blei (2011); Hoffman et al. (2013); Broderick et al. (2013)). Although these techniques often have excellent predictive performance, it is well known (Bishop, 2006) that variational methods tend to substantially underestimate posterior uncertainty and provide a poor characterization of posterior dependence, while lacking theoretical guarantees.

Our approach instead falls in a third class of methods which avoid extensive communication among machines by running independent MCMC chains for each data subset and obtaining draws from subset posteriors. These subset posteriors can be combined in a variety of ways. Some of these methods simply average draws from each subset (Scott et al., 2013). Other alternatives use an approximation to the full posterior distribution based on kernel density estimates (Neiswanger, Wang and Xing, 2013) or the so-called Weierstrass transform (Wang and Dunson, 2013). These methods have limitations related to the dimension of the parameter, moreover, their applicability and theoretical justification are restricted to parametric models. Unlike the method proposed below, none of the aforementioned algorithms are provably robust.

Our work was inspired by recent multivariate median-based techniques for robust estimation developed in Minsker (2013) (see also Hsu and Sabato (2013); Alon, Matias and Szegedy (1996); Lerasle and Oliveira (2011); Nemirovski and Yudin (1983) where similar ideas were applied in different frameworks).

Preliminaries

We proceed by recalling key definitions and facts which will be used throughout the paper.

Finally, given two nonnegative sequences $\{a_{n}\}$ and $\{b_{n}\}$, we write $a_{n}\lesssim b_{n}$ if $a_{n}\leq Cb_{n}$ for some $C>0$ and all $n$. Other objects and definitions are introduced in the course of exposition when necessity arises.

2 Generalizations of the univariate median

$j_{\ast}:=j(\varepsilon_{\ast}),\text{ where ties are broken arbitrarily, }$ and set

We will say that $x_{\ast}$ is the metric median of $x_{1},\ldots,x_{m}$. Note that $x_{\ast}$ always belongs to $\{x_{1},\ldots,x_{m}\}$ by definition. Advantages of this definition are its generality (only metric space structure is assumed) and simplicity of numerical evaluation since only the pairwise distances $d(x_{i},x_{j}),\ i,j=1,\ldots,m$ are required to compute the median. This construction was previously employed in Nemirovski and Yudin (1983) in the context of stochastic optimization and is further studied in Hsu and Sabato (2013). A closely related notion of the median was used in Lopuhaa and Rousseeuw (1991) under the name of the “minimal volume ellipsoid” estimator.

Finally, we recall an important property of the median (shared both by $\mbox{{\rm med}}_{g}$ and $\mbox{{\rm med}}_{0}$) which states that it transforms a collection of independent, “weakly concentrated” estimators into a single estimator with significantly stronger concentration properties. Given $q,\alpha$ such that $0<q<\alpha<1/2$, define a nonnegative function $\psi(\alpha,q)$ via

The following result is an adaptation of Theorem 3.1 in Minsker (2013):

Let $\hat{\theta}_{\ast}=\mbox{{\rm med}}_{g}(\hat{\theta}_{1},\ldots,\hat{\theta}_{m})$ be the geometric median of $\{\hat{\theta}_{1},\ldots,\hat{\theta}_{m}\}$. Then

Let $\hat{\theta}_{\ast}=\mbox{{\rm med}}_{0}(\hat{\theta}_{1},\ldots,\hat{\theta}_{m})$. Then

While we require $\kappa<1/3$ above for clarity and to keep the constants small, we prove a slightly more general result that holds for any $\kappa<1/2$.

Theorem 2.1 implies that the concentration of the geometric median of independent estimators around the “true” parameter value improves geometrically fast with respect to the number of such estimators, while the estimation rate is preserved, up to a constant. In our case, the role of $\hat{\theta}_{j}$’s will be played by posterior distributions based on disjoint subsets of observations, viewed as elements of the space of signed measures equipped with a suitable distance.

Parameter $\kappa$ allows taking corrupted observations into account: if the initial sample contains not more than $\lfloor\kappa m\rfloor$ outliers (of arbitrary nature), then at most $\lfloor\kappa m\rfloor$ estimators amongst $\{\theta_{1},\ldots,\theta_{m}\}$ can be affected but their median remains stable, still being close to the unknown $\theta_{0}$ with high probability. To clarify the notion of “robustness” that such a statement provides, assume that $\hat{\theta}_{1},\ldots,\hat{\theta}_{m}$ are consistent estimators of $\theta_{0}$ based on disjoint samples of size $n/m$ each. If $\frac{n}{m}\to\infty$, then $\frac{\kappa m}{n}\to 0$, hence the breakdown point of the estimator $\hat{\theta}_{\ast}$ is is general. However, it is able to handle a number of outliers that grows like $o(n)$ while preserving consistency, which is the best one can hope for without imposing any additional assumptions on the underlying distribution, parameter of interest or nature of the outliers.

Let us also mention that the the geometric median of a collection of points in a Hilbert space belongs to the convex hull of these points. Thus, one can think about “downweighing” some observations (potential outliers) and increasing the weight of others, and geometric median gives a way to formalize this approach. The median $\mbox{{\rm med}}_{0}$ defined in (2.3) corresponds to the extreme case when all but one weight are equal to . Its potential advantage lies in the fact that its evaluation requires only the knowledge of pairwise distances $d(\hat{\theta}_{i},\hat{\theta}_{j}),\ i,j=1,\ldots,m$, see (2.2).

3 Distances between probability measures

Next, we discuss the special family of distances between probability measures that will be used throughout the paper. These distances provide the necessary structure to define and evaluate medians in the space of measures, as discussed above. Since one of our goals was to develop computationally efficient techniques, we focus on distances that admit accurate numerical approximation.

Important special cases include the situation when

where $\|f\|_{L}:=\sup\limits_{x_{1}\neq x_{2}}\frac{|f(x_{1})-f(x_{2})|}{\rho(x_{1},x_{2})}$ is the Lipschitz constant of $f$.

It is well-known (Dudley (2002), Theorem 11.8.2) that in this case $\|P-Q\|_{\mathcal{F}_{L}}$ is equal to the Wasserstein distance (also known as the Kantorovich-Rubinstein distance)

where $\mathcal{L}(\bm{Z})$ denotes the law of a random variable $\bm{Z}$ and the infimum on the right is taken over the set of all joint distributions of $(\bm{X},\bm{Y})$ with marginals $P$ and $Q$.

Note that when $P$ and $Q$ are discrete measures (e.g., $P=\sum\limits_{j=1}^{N_{1}}\beta_{j}\delta_{z_{j}}$ and $Q=\sum\limits_{j=1}^{N_{2}}\gamma_{j}\delta_{y_{j}}$), then

In this paper, we will only consider characteristic kernels, which means that $\|P-Q\|_{\mathcal{F}_{k}}=0$ if and only if $P=Q$. It follows from Theorem 7 in Sriperumbudur et al. (2010) that a sufficient condition for $k$ to be characteristic is its strict positive definiteness: we say that $k$ is strictly positive definite if it is bounded, measurable, and such that for all non-zero signed Borel measures $\nu$

If $\phi$ is compactly supported, then $k$ is characteristic.

For i.i.d samples, a useful and favorable fact is that $e_{m,n}$ often does not depend on $D$: under weak assumptions on kernel $k$, $e_{m,n}$ has an upper bound of order $m^{-1/2}+n^{-1/2}$ (that is, $\lim_{m,n\to\infty}\Pr\left(e_{m,n}\geq C(m^{-1/2}+n^{-1/2})\right)$ can be made arbitrarily small by choosing $C$ big enough, see Corollary 12 in Sriperumbudur et al. (2009)). On the other hand, the bound for the (stronger) Wasserstein distance is not dimension-free and is of order $m^{-1/(D+1)}+n^{-1/(D+1)}$. Similar error rates hold for empirical measures based on samples from Markov Chains used to approximate invariant distributions, including MCMC samples (see Boissard and Le Gouic (2014) and Fournier and Guillin (2013)).

it will be useful to assume that the class $\mathcal{F}$ is chosen such that the distance between the measures is lower bounded by the distance between their means, namely

Then $k$ is characteristic and satisfies (2.13) with $C=1$.

The total variation distance between two probability measures defined on a $\sigma$-algebra $\mathfrak{B}$ is

Contributions and main results

This section explains the construction of “median posterior” (or M-Posterior) distribution, along with the theoretical guarantees for its performance.

and assume that the metric space $(\Theta,\rho)$ is separable.

Let $k$ be a characteristic kernel defined on $\Theta\times\Theta$. Kernel $k$ defines a metric on $\Theta$ via

All subsequent results apply to this special case. While this is a “natural” metric for the problem, the disadvantage of $k_{H}(\cdot,\cdot)$ is that it is often difficult to evaluate numerically. Instead, we will consider metrics $\rho_{k}$ that are “dominated” by $\rho$ (this is formalized in assumption 3.4).

for all Borel measurable sets $B\subseteq\Theta$. It is known (see Ghosal, Ghosh and Van Der Vaart (2000)) that under rather general assumptions the posterior distribution $\Pi_{n}$ “contracts” towards $\theta_{0}$, meaning that

almost surely or in probability as $n\to\infty$ for a suitable sequence $\varepsilon_{n}\to 0$.

One of the questions that we address can be formulated as follows: what happens if some observations in $\mathcal{X}_{n}$ are corrupted, e.g., if $\mathcal{X}_{n}$ contains outliers of arbitrary nature and magnitude? Even if there is only one “outlier”, the usual posterior distribution might concentrate most of its mass “far” from the true value $\theta_{0}$.

We proceed with a general description of our proposed algorithm for constructing a robust version of the posterior distribution. Let $1\leq m\leq n/2$ be an integer. Divide the sample $\mathcal{X}_{n}$ into $m$ disjoint groups $G_{1},\ldots,G_{m}$ of size $|G_{j}|\geq\lfloor n/m\rfloor$ each:

A good choice of $m$ efficiently exploits the available computational resource while ensuring that the groups $G_{j}$s are sufficiently large.

Let $\Pi$ be a prior distribution over $\Theta$, and let

be the family of subset posterior distributions depending on disjoint subgroups $G_{j},\ j=1,\ldots,m$:

where the medians $\mbox{{\rm med}}_{g}(\cdot)$ and $\mbox{{\rm med}}_{0}(\cdot)$ are evaluated with respect to $\|\cdot\|_{\mathcal{F}_{L}}$ or $\|\cdot\|_{\mathcal{F}_{k}}$ introduced in section 2.2 above. Note that $\hat{\Pi}_{n,g}$ and $\hat{\Pi}_{n,0}$ are always probability measures: indeed, due to the aforementioned properties of a geometric median, there exists $\alpha_{1}\geq 0,\ldots,\alpha_{m}\geq 0,\ \sum\limits_{j=1}^{m}\alpha_{j}=1$ such that $\hat{\Pi}_{n,g}=\sum\limits_{j=1}^{m}\alpha_{j}\Pi^{(j)}$, and $\hat{\Pi}_{n,0}\in\{\Pi^{(1)}(\cdot),\ldots,\Pi^{(m)}(\cdot)\}$ by definition.

To overcome this difficulty, we propose a modification of our approach where the random measures $\Pi_{n}^{(j)}$ are replaced by the stochastic approximations $\Pi_{|G_{j}|,m}(\cdot|G_{j})$, $j=1,\ldots,m$ of the full posterior distribution. To this end, define the “stochastic approximation” based on the subsample $G_{j}$ as

where we assume that $p_{\theta}^{m}(\cdot)$ is an integrable function for all $\theta$. In other words, $\Pi_{|G_{j}|,m}(\cdot|G_{j})$ is obtained as a posterior distribution given that each data point from $G_{j}$ is observed $m$ times. While each of $\Pi_{|G_{j}|,k}(\cdot|G_{j})$ might underestimate uncertainly, the median $\hat{\Pi}^{\mbox{\rm st}}_{n,g}$ (or $\hat{\Pi}^{\mbox{\rm st}}_{n,0}$) of these random measures yields credible sets with much better coverage. This approach shows good performance in numerical experiments. One of our main results (see section 3.3) provides a justification for this observation (albeit, under rather strong assumptions and for the parametric case).

2 Convergence of posterior distribution and robust Bayesian inference

In this subsection, we study the contraction and robustness properties of the median posterior.

Our first result establishes the “weak concentration” property of the posterior distribution around the true parameter. Let $\delta_{0}:=\delta_{\theta_{0}}$ be the Dirac measure supported on $\theta_{0}\in\Theta$. Recall the following version of Theorem 2.1 in Ghosal, Ghosh and Van Der Vaart (2000) (we state the result for the Wasserstein distance $d_{W_{1,\rho}}(\Pi_{n}(\cdot|\mathcal{X}_{l}),\delta_{0})$ rather than the (closely related) contraction rate of the posterior distribution). Here, the Wasserstein distance is evaluated with respect to the “Hellinger metric” $\rho(\cdot,\cdot)$ defined in (3.1).

Let $\mathcal{X}_{l}=\{X_{1},\ldots,X_{l}\}$ be an i.i.d. sample from $P_{0}$. Assume that $\varepsilon_{l}>0$ and $\Theta_{l}\subset\Theta$ are such that for some constant $C>0$

The proof closely mimics the argument behind Theorem 2.1 in Ghosal, Ghosh and Van Der Vaart (2000). Details are outlined in section A.2. ∎

Conditions of Theorem 3.1 are standard assumptions guaranteeing that the resulting posterior distribution contracts to the true parameter $\theta_{0}$ at the rate $\varepsilon_{n}$. Note that the bounds for the distance $d_{W_{1,\rho}}(\delta_{0},\Pi_{l}(\cdot|\mathcal{X}_{l})$ slightly differ from the contraction rate itself: indeed, we have

hence to obtain the inequality $d_{W_{1,\rho}}(\delta_{0},\Pi_{l}(\cdot|\mathcal{X}_{l}))\lesssim\varepsilon_{l}$, we usually require

which adds an extra logarithmic factor in the parametric case.

Therefore, $\mathcal{F}_{k}\subseteq\mathcal{F}_{L}$ and $\|P-Q\|_{\mathcal{F}_{k}}\leq\|P-Q\|_{\mathcal{F}_{L}}$, where $\mathcal{F}_{k}$ and $\mathcal{F}_{L}$ were defined in (2.10) and (2.8) respectively, and the underlying metric structure is given by $\rho$. In particular, convergence with respect to $\|\cdot\|_{\mathcal{F}_{L}}$ implies convergence with respect to $\|\cdot\|_{\mathcal{F}_{k}}$.

Let $X_{1},\ldots,X_{n}$ be an i.i.d. sample from $P_{0}$, and assume that $\hat{\Pi}_{n,g}$ is defined with respect to the norm $\|\cdot\|_{\mathcal{F}_{L}}$ as in (3.4) above. Set $l:=\lfloor n/m\rfloor$, assume that conditions of Theorem 3.1 hold, and, moreover, that $\varepsilon_{l}$ satisfies

It is enough to apply part (a) of Theorem 2.1 with $\kappa=0$ to the independent random measures $\Pi_{n}(\cdot|G_{j}),\ j=1,\ldots,m$. Note that the “weak concentration” assumption (A.2) is implied by (3.6). ∎

Once again, note the exponential improvement of concentration as compared to Theorem 3.1. It is easy to see that a similar statement holds for the median $\hat{\Pi}_{n,0}(\cdot)$ defined in (3.4) (even for the stronger Wasserstein distance $d_{W_{1,\rho}}(\delta_{0},\hat{\Pi}_{n,0})$), modulo changes in constants.

While the result of the previous statement is promising, numerical approximation and sampling from the “robust posterior” $\hat{\Pi}_{n,g}$ is often problematic due to the underlying geometry defined by the Hellinger metric, and the associated distance $\|\cdot\|_{\mathcal{F}_{k_{H}}}$ is hard to estimate in practice. Our next goal is to derive similar guarantees for the M-posterior evaluated with respect to the computationally tractable family of distances discussed in section 2.3 above.

To transfer the conclusions of Theorem 3.1 and Corollary 3.2 to the case of other kernels $k(\cdot,\cdot)$ and associated metrics $\rho_{k}(\cdot,\cdot)$, we need to guarantee the existence of tests versus the complements of the balls in these distances. Such tests can be obtained from comparison inequalities between distances.

where $d$ is the Hellinger distance or the Euclidean distance (in the parametric case).

When $d$ is the Euclidean distance, we will impose an additional mild assumption guaranteeing existence of test versus the complements of the balls (for the Hellinger distance, this is always true, see Ghosal, Ghosh and Van Der Vaart (2000)). Namely, we will assume that for every $n$ and every pair $\theta_{1},\theta_{2}\in\Theta$, there exists a test $\phi_{n}:=\phi_{n}(X_{1},\ldots,X_{n})$ such that for some $\gamma>0$ and a universal constant $K>0$

Below, we provide several examples of kernels satisfying the stated assumption.

For finite-dimensional models, we will be especially interested in kernels $k(\cdot,\cdot)$ such that the associated metric $\rho_{k}(\cdot,\cdot)$ is bounded by the Euclidean distance. The following proposition gives a sufficient condition for this to hold.

Moreover, the result of the previous proposition clearly remains valid for kernels of the form

We are ready to state our main result for convergence with respect to the RKHS-induced distance $\|\cdot\|_{\mathcal{F}_{k}}$.

Assume that conditions of Theorem 3.1 hold with $\rho$ being the Hellinger or the Euclidean distance, and that assumption 3.4 is satisfied. In addition, let prior $\Pi$ be such that

The result essentially follows from the combination of Theorem 3.1 and assumption 3.4, see section A.3 in the appendix for details. ∎

Theorem 3.8 yields the “weak” estimate that is needed to obtain the stronger bound for the M-Posterior distribution $\hat{\Pi}_{n,g}$. This is summarized in the following corollary:

Let $X_{1},\ldots,X_{n}$ be an i.i.d. sample from $P_{0}$, and assume that $\hat{\Pi}_{n,g}$ is defined with respect to the distance $\|\cdot\|_{\mathcal{F}_{k}}$ as in (3.4) above. Let $l:=\lfloor n/m\rfloor$. Assume that conditions of Theorem 3.8 hold, and, moreover, $\varepsilon_{l}$ is such that

It is enough to apply parts (a) and (b) of Theorem 2.1 with $\kappa=0$ to the independent random measures $\Pi_{|G_{j}|}(\cdot|G_{j}),\ j=1,\ldots,m$. Note that the “weak concentration” assumption (A.1) is implied by (3.10). ∎

where $\theta_{\ast}=\int_{\Theta}\theta d\hat{\Pi}_{n,g}(\theta)$ is the mean of $\hat{\Pi}_{n,g}$. In other words, this shows that the M-posterior mean is the “robust” estimator of $\theta_{0}$.

3 Bayesian inference based on stochastic approximation of the posterior distribution

As we have already mentioned in section 3.1, when the number of disjoint subgroups $m$ is large, the resulting M-Posterior distribution is “too flat”, which results in large credible sets and overestimation of uncertainty. Clearly, the source of the problem is the fact that each individual random measure $\Pi_{|G_{j}|}(\cdot|G_{j}),\ j=1,\ldots,m$ is based on sample of size $l\simeq\frac{n}{m}$ which can be much smaller than $n$.

where we have assumed that $p_{\theta}(\cdot)$ is integrable. Here, $\left(\prod_{i\in G_{j}}p_{\theta}(X_{i})\right)^{m}$ can be viewed as an approximation of the full data likelihood. We call the random measure $\Pi_{|G_{j}|,m}(\cdot|G_{j})$ the $j$-th stochastic approximation to the full posterior distribution.

Of course, such a “correction” negatively affects coverage properties of the credible sets associated with each measure $\Pi_{|G_{j}|}(\cdot|G_{j})$. However, taking the median of stochastic approximations yields improved coverage of the resulting M-posterior distribution. The main goal of this section is to establish an asymptotic statement in spirit of a Bernstein-von Mises theorem for the M-posterior based on stochastic approximations $\Pi_{|G_{j}|,m}(B|G_{j}),\ j=1,\ldots,m$.

We will start by showing that under certain assumptions the upper bounds for the convergence rates of $\Pi_{|G_{j}|,m}(\cdot|G_{j})$ towards $\delta_{0}$ are the same as for $\Pi_{|G_{j}|}(\cdot|G_{j})$, the “standard” posterior distribution given $G_{j}$.

be the bracketing entropy. In what follows, $B(\theta_{0},r):=\{\theta\in\Theta:\ h(P_{\theta},P_{\theta_{0}})\leq r\}$ denotes the “Hellinger ball” of radius $r$ centered at $\theta_{0}$.

There exist constants $c_{j},\ j=1,\ldots,4$ and $\zeta>0$ such that if

In particular, one can choose $c_{1}=1/24,\ c_{2}=(4/27)(1/1926),\ c_{3}=10$ and $c_{4}=(2/3)^{5/2}/512$.

whenever $h\left(P_{\theta},P_{\theta_{0}}\right)\leq r_{0}$, and

there exists $\alpha>0$ such that for $\theta_{1},\theta_{2}\in B(\theta_{0},r_{0})$,

Application of theorem 3.10 to the analysis of “stochastic approximations” yields the following result.

Let $\varepsilon_{l}>0$ be such that conditions of Theorem 3.10 hold with $\zeta:=\varepsilon_{l}$, and

Note that for the kernel $k(\cdot,\cdot)$ of the form (3.9), assumption 3.4 reduces to the inequality between the Hellinger and Euclidean distances.

As before, Theorem 2.1 combined with the “weak concentration” inequality of Theorem 3.11 gives stronger guarantees for the median $\hat{\Pi}^{\mbox{\rm st}}_{n,g}$ (or its alternative $\hat{\Pi}^{\mbox{\rm st}}_{n,0}$) of $\Pi_{|G_{1}|,m}(\cdot|G_{1}),\ldots,\Pi_{|G_{m}|,m}(\cdot|G_{m})$. Exact statement is very similar in spirit to Corollary 3.9.

We will first state a preliminary result for each individual “subset posterior” distribution:

Let $X_{1},\ldots,X_{l}$ be an i.i.d. sample from $P_{\theta_{0}}$ for some $\theta_{0}$ in the interior of $\Theta$. Assume that

the family $\{P_{\theta},\ \theta\in\Theta\}$ is differentiable in quadratic mean;

the prior $\Pi$ has a density (with respect to the Lebesgue measure) that is continuous and positive in the neighborhood of $\theta_{0}$;

conditions of Theorem 3.10 hold with $\zeta=\frac{C}{\sqrt{l}}$ for some $C>0$ and $l$ large enough.

in $P_{\theta_{0}}$-probability as $l\to\infty$.

The proof follows standard steps (e.g., Theorem 10.1 in Van der Vaart (2000)), where the existence of tests is substituted by the inequality of Theorem 3.10. See section A.5 in the appendix for more details. ∎

The implication of this result for the M-posterior is the following: if $k$ is the kernel of type (3.9), for sufficiently regular parametric families (differentiable in quadratic mean, with “well-behaved” bracketing numbers, satisfying assumption 3.4 for the Euclidean distance with $\gamma=1$) and regular priors, then

the M-posterior is well approximated by a normal distribution centered at the “robust” estimator $\theta^{\ast}$ of unknown $\theta_{0}$;

the estimator $\theta^{\ast}$ is a center of the confidence set of level $1.15^{-m}$ and diameter of order $\sqrt{\frac{m}{n}}$ (same as we would expect for this level for the usual posterior distribution - however, the bound for the M-posterior holds for finite sample sizes).

in $P_{\theta_{0}}$-probability when $n\to\infty$, where $\theta^{\ast}$ is the mean of $\hat{\Pi}^{\mbox{\rm st}}_{n,0}$.

Assume that conditions (a), (b) of Theorem 3.11 hold with

and $\gamma=1$. Then for all $n\geq n_{0}$ and $\bar{R}$ large enough,

(a) It is easy to see that convergence in total variation norm, together with an assumption that the prior distribution satisfies

implies that the expectations converge in $P_{\theta_{0}}$-probability as well:

Together with an observation that the total variation distance between $N(\mu_{1},\Sigma)$ and $N(\mu_{2},\Sigma)$ is bounded by the multiple of $\left\|\mu_{1}-\mu_{2}\right\|_{2}$, it implies that we can replace $\theta_{0}+\frac{\Delta_{l,\theta_{0}}}{\sqrt{l}}$ by the mean

in other words, the conclusion of Proposition 3.13 can be stated as

in $P_{\theta_{0}}$-probability. Now assume that $m=\lfloor\frac{n}{l}\rfloor$ is fixed, and let $n,l\to\infty$. As before, let $G_{1},\ldots,G_{m}$ be disjoint groups of i.i.d. observations from $P_{\theta_{0}}$ of cardinality $l$ each. Recall that, by the definition (2.3) of $\mbox{{\rm med}}_{0}(\cdot)$, $\Pi^{\mbox{\rm st}}_{n,0}=\Pi_{l,m}(\cdot|\mathcal{X}_{l_{\ast}})$ for some $l_{\ast}\leq m$, and $\theta^{\ast}:=\theta_{l_{\ast},m}$ is the mean of $\Pi^{\mbox{\rm st}}_{n,0}$. Clearly, we have

(b) Let $\varepsilon_{l}\geq C\sqrt{\frac{1}{l}}$ where $C$ large enough so that

In particular, for $m=A\log(n)$ and $\varepsilon_{l}\simeq\sqrt{\frac{m}{n}}$, we obtain the bound

for some constant $\bar{R}$ independent of $m$. Note that $\theta^{\ast}$ itself depends on $m$, hence this bound is not uniform, and holds only for a given confidence level $1-n^{-A}$.

It is convenient to interpret this (informally) in terms of the credible sets: to obtain the credible set with “frequentist” coverage level $\geq 1-n^{-A}$, pick $m=A\log n$ and use the $(1-n^{-A})$ - credible set of the M-posterior $\hat{\Pi}^{\mbox{\rm st}}_{n,0}$.

Numerical algorithms and examples

In this section, we consider examples and applications in which comparisons are made for the inference based on the usual posterior distribution and on the M-Posterior. One of the well-known and computationally efficient ways to find the geometric median in Hilbert spaces is the famous Weiszfeld’s algorithm (introduced in Weiszfeld (1936)). Details of implementation are described in Algorithms 1 and 2. Algorithm 1 is a particular case of Weiszfeld’s algorithm applied to subset posterior distributions and distance $\|\cdot\|_{\mathcal{F}_{k}}$, while Algorithm 2 shows how to obtain an approximation to M-Posterior given the samples from $\Pi_{n,m}(\cdot|G_{j}),\ j=1\ldots m$. Note that the subset posteriors $\Pi_{n,m}(\cdot|G_{j})$ whose “weights” $w_{\ast,j}$ in the expression of the M-Posterior are small (in our case, smaller than $1/(2m)$) are excluded from the analysis. Our extensive simulations show the empirical evidence in favor of this additional thresholding step.

Detailed discussion of convergence rates and acceleration techniques for Weiszfeld’s method from the viewpoint of modern optimization can be found in Beck and Sabach (2013). For alternative approaches and extensions of Weiszfeld’s algorithm, see Bose, Maheshwari and Morin (2003), Ostresh (1978), Overton (1983), Chandrasekaran and Tamir (1990), Cardot, Cénac and Zitt (2012), Cardot, Cénac and Zitt (2013), among other works.

In all numerical simulations below, we use “stochastic approximations” and the corresponding median measure $\hat{\Pi}^{\mbox{\rm st}}_{n,g}$, unless noted otherwise.

Before presenting the results of numerical analysis, let us remark on two important computational aspects.

The number of subsets $m$ appears as a “free parameter” entering the theoretical guarantees for M-Posterior. One interpretation of $m$ (in terms of the credible sets) is given in the end of section 3.3. Our results also imply that partitioning the data into $m=2k+1$ subsets guarantees robustness to the presence of $k$ outliers of arbitrary nature.

In many applications, $m$ is dictated by the sample size and computational resources (e.g., the number of available machines). In section B.3 of the appendix, we describe a heuristic approach to selection of $m$ that shows good practical performance. As a rule of a thumb, we recommend choosing $m\lesssim\sqrt{n}$ as larger values of $m$ lead to an M-posterior that overestimates uncertainty. This heuristic is supported by the numerical results presented below.

It is easy to get a general idea regarding the potential improvement in computational time complexity achieved by the M-Posterior. Given the data set $\mathcal{X}_{n}=\{X_{1},\ldots,X_{n}\}$ of size $n$, let $t(n)$ be the running time of the algorithm (e.g., MCMC) that outputs a single observation from the posterior distribution $\Pi_{n}(\cdot|\mathcal{X}_{n})$. If the goal is to obtain $S$ samples from the posterior, then the total running time is $O\left(S\cdot t(n)\right)$. Let us compare this time with the running time needed to obtain $S$ samples from the $M$-posterior given that the algorithm is running on $m$ machines in parallel. In this case, we need to generate $O\left(S\right)$ samples from each of $m$ subset posteriors, which is done in time $O\left(S\cdot t\left(\frac{n}{m}\right)\right)$, where $S$ is typically large and $m\ll n$. According to Theorem 7.1 in Beck and Sabach (2013), Weiszfeld’s algorithm approximates the M-Posterior to degree of accuracy $\varepsilon$ in at most $O(1/\varepsilon)$ steps, and each of these steps has complexity $O(S^{2})$ (which follows from (2.12)), so that the total running time is

The term $\frac{S^{2}}{\varepsilon}$ can be refined in several ways via application of more advanced optimization techniques (see the aforementioned references). If, for example, $t(n)\simeq n^{r}$ for some $r\geq 1$, then $\frac{S}{m}\cdot t\left(\frac{n}{m}\right)\simeq\frac{1}{m^{1+r}}Sn^{r}$ which should be compared to $S\cdot n^{r}$ required by the standard approach.

To give a specific example, consider an application of (4.1) in the context of Gaussian process (GP) regression. If $n$ is the number of training samples, then GP regression has $O(n^{3})+O(Sn^{2})$ asymptotic time complexity to obtain $S$ samples from the posterior distribution of GP (Rasmussen and Williams, 2006, Algorithm 2.1). Assuming we have access to $m$ machines, the time complexity to obtain $S$ samples from M-Posterior in GP regression is $O\left(\left(\tfrac{n}{m}\right)^{3}+S\left(\tfrac{n}{m}\right)^{2}+\tfrac{S^{2}}{\varepsilon}\right)$. If for example $S=cn$ for some $c>0$ and $m^{2}<n\varepsilon$, we get $O(m^{2})$ improvement in running time.

In many cases, replacing the “subset posterior” by the stochastic approximation does not result in increased sampling complexity: indeed, the log-likelihood in the sampling algorithm for the subset posterior is simply multiplied by $m$ to obtain the sampler for the stochastic approximation. We have included the description of a modified Dirichlet mixture model in section B.2 of the appendix as an illustration.

M-Posterior was more robust than its competitors and its performance improved with increasing magnitude of the outlier. We compared the performance of “consensus posterior”, the overall posterior, and the M-Posterior using the empirical coverage of (1-$\alpha$)100% credible intervals (CIs) calculated across 50 replications for $\alpha=0.2,0.15,0.10$, and 0.05. The empirical coverages of M-Posterior’s CIs showed robustness to magnitude of the outlier. On the contrary, performance of the consensus and overall posteriors deteriorated fairly quickly across all $\alpha$’s leading to 0% empirical coverage as magnitude of the outlier increased from $i=1$ to $i=25$ (Figure 1). We compared uncertainty quantification of the M-Posterior with that of the overall posterior using relative lengths of their CIs, with zero value corresponding to identical lengths and a positive value to wider CIs of the M-Posterior. We found that widths of CIs for both posteriors were fairly similar for $i=1,\ldots,25$, with M-Posterior’s CIs being slightly wider in absence of large outliers (Figure 2).

Stochastic approximation was important for proper calibration of uncertainty quantification. The empirical coverages of (1-$\alpha$)100% CIs of the M-Posterior without stochastic approximation overcompensated for uncertainty at all levels of $\alpha$ (Figure 3). Similarly, lengths of the CIs of M-Posterior without stochastic approximation are wider than those with stochastic approximation (Figure 4). Both these observations showed that stochastic approximation led to shorter CIs for M-Posterior that had empirical coverages close to their theoretical values.

The number of subsets ($m$) had an effect on credible interval length of the M-Posterior. We modified the simulation above and generated 1000 observations $\operatorname{\mathbf{x}}$, with the last 10 observations in $\operatorname{\mathbf{x}}$ being outliers with value $x_{j}=25\max(|x_{1}|,\ldots,|x_{990}|)$, $j=991,\ldots,1000$. The simulation setup was replicated 50 times. M-Posteriors were obtained for $m=16,18,\ldots,40,50,60$. Across all values of $m$, M-Posterior’s CI was compared to the CI of the overall posterior after removing the outliers; the relative difference of M-Posterior and the overall posterior CI lengths decreases for $m\geq 22>2k$, where $k=10$ is the number of outliers, remains stable as $m$ increases to $m=38$, and grows for larger values of $m$ (Figure 5a). This demonstrates that inference based on M-Posterior was not too sensitive to the choice of $m$ for a wide range of values.

2 Real data analysis: General social survey

The General Social Survey (GSS; gss.norc.org) has collected responses to questions about evolution of American society since 1972. We selected data for 9 questions from different social topics: “happy” (happy), “Bible is a word of God” (bible), “support capital punishment” (cap), “support legalization of marijuana” (grass), “support premarital sex” (pre), “approve bible prayer in public schools” (prayer), “expect US to be in world war in 10 years” (uswar), “approve homosexual sex relations” (homo), “support abortion” (abort). These questions were in the survey since 1988 and their answers were converted to two levels: yes or no. Missing data were imputed based on the average, resulting in a data set with approximately 28,000 respondents.

M-Posterior had similar uncertainty quantification as the overall posterior while being more efficient. M-Posterior was at least 10 ($m=20$) and 8 times ($m=10$) faster than the overall posterior and it used less than 25% of the memory resources required by the overall posterior (Figure 5b). The overall posterior was more concentrated than the M-Posterior for $m=10$ and 20; however, its coverage of maximum likelihood estimators for $\operatorname{{\bm{\pi}}}_{c_{i},c_{j}}$ obtained from the test data was worse than that of the two M-Posteriors (Table 1).

References

Appendix A Remaining proofs.

We will prove a slightly more general result:

Let $\hat{\theta}_{\ast}=\mbox{{\rm med}}_{g}(\hat{\theta}_{1},\ldots,\hat{\theta}_{m})$ be the geometric median of $\{\hat{\theta}_{1},\ldots,\hat{\theta}_{m}\}$. Then

where $C_{\alpha}=(1-\alpha)\sqrt{\frac{1}{1-2\alpha}}$.

Let $\hat{\theta}_{\ast}=\mbox{{\rm med}}_{0}(\hat{\theta}_{1},\ldots,\hat{\theta}_{m})$. Then

To get the bound stated in the paper, take $q=\frac{1}{7}$ and $\alpha=\frac{3}{7}$ in part (a) and $q=\frac{1}{4}$ in part b.

We start by proving part a. To this end, we will need the following lemma (see lemma 2.1 in Minsker (2013)):

and $r>0$. Then there exists a subset $J\subseteq\left\{1,\ldots,m\right\}$ of cardinality $|J|>\alpha m$ such that for all $j\in J$, $\|x_{j}-z\|>r$.

Assume that event $\mathcal{E}:=\left\{\|\hat{\theta}_{\ast}-\theta_{0}\|>C_{\alpha}\varepsilon\right\}$ occurs. Lemma A.1 implies that there exists a subset $J\subseteq\{1,\ldots,m\}$ of cardinality $|J|\geq\alpha k$ such that $\|\hat{\theta}_{j}-\theta_{0}\|>\varepsilon$ for all $j\in J$, hence

If $W$ has Binomial distribution $W\sim B(\lfloor(1-\gamma)m\rfloor+1,q)$, then

(see Lemma 23 in Lerasle and Oliveira (2011) for a rigorous proof of this fact). Chernoff bound (e.g., Proposition A.6.1 in van der Vaart and Wellner (1996)), together with an obvious bound $\lfloor(1-\gamma)m\rfloor+1>(1-\gamma)m$, implies that

To establish part b, we proceed as follows: let $\mathcal{E}_{1}$ be the event

The rest of the proof repeats the argument of part a since

where $\mathcal{E}_{1}^{c}$ is the complement of $\mathcal{E}_{1}$.

A.2 Proof of Theorem 3.1

By the definition of Wasserstein distance $d_{W_{1}}$,

(recall that $\rho$ is the Hellinger distance). Let $R$ be a large enough constant to be determined later. Note that the Hellinger distance is uniformly bounded by 1. Using (A.3), it is easy to see that

To this end, it remains to estimate the second term in the sum above. We will follow the proof of Theorem 2.1 in Ghosal, Ghosh and Van Der Vaart (2000). Bayes formula implies that

For any $C_{1}>0$, Lemma 8.1 Ghosal, Ghosh and Van Der Vaart (2000) yields

for every probability measure $Q$ on the set $A_{l}$. Moreover, by the assumption on the prior $\Pi$,

Consequently, with probability at least $1-\frac{1}{C_{1}^{2}l\varepsilon_{l}^{2}}$,

Define the event B_{l}=\left\{\int\limits_{\Theta}\prod\limits_{i=1}^{l}\frac{p_{\theta}}{p_{0}}(X_{i})d\Pi(\theta)\leq\exp\Big{(}-(1+C_{1}+C)l\varepsilon_{l}^{2}\Big{)}\right\}.

Let $\Theta_{l}$ be the set satisfying conditions of Theorem 3.1. Then by Theorem 7.1 in Ghosal, Ghosh and Van Der Vaart (2000), there exist test functions $\phi_{l}:=\phi_{l}(X_{1},\ldots,X_{l})$ and a universal constant $K$ such that

Next, by the definition of $B_{l}$, we have

To estimate the second term of last equation, note that

for $R\geq\sqrt{(C+4)/K}$. Set $C_{1}=1$ and note that $I\{B_{l}\}=1$ with probability $P(B_{l})\leq 1/l\varepsilon_{l}^{2}$. It follows from (A.6), (A.2) and (A.2) and Chebyshev’s inequality that for any $t>0$

A.3 Proof of Theorem 3.8

From equation (3.7) in the paper and proceeding as in proof of Theorem 3.2, we get that

Following in the proof of Theorem 3.1, we note that $\Pr(B_{l})\geq 1-\frac{1}{C_{1}^{2}l\varepsilon_{l}^{2}}$ (where constants $C,C_{1}$ are the same as in the proof of Theorem 3.1). Also, note that

By Theorem 7.1 in Ghosal, Ghosh and Van Der Vaart (2000), there exist test functions $\phi_{l}:=\phi_{l}(X_{1},\ldots,X_{l})$ and a universal constant $K$ such that

Repeating the steps of the proof of Theorem 3.1, we see that if $R$ is chosen large enough, then

A.4 Proof of Theorem 3.11

The proof strategy is similar to Theorem 3.1. Note that

where $h(\cdot,\cdot)$ is the Hellinger distance.

Let $\mathcal{E}_{l}:=\{\theta:h(P_{\theta},P_{0})\geq R\varepsilon_{l}\}$. By the definition of $\Pi_{n,m}$, we have

To bound the denominator from below, we proceed as before. Let

where $Q$ is a probability measure supported on $\Theta_{l}$. Lemma 8.1 in Ghosal, Ghosh and Van Der Vaart (2000) yields that $\Pr(B_{l})\leq\frac{1}{l\varepsilon_{l}^{2}}$ for any $Q$, in particular, for the conditional distribution $\Pi(\cdot|\Theta_{l})$. We conclude that

To estimate the numerator in (A.11), note that if Theorem 3.10 holds for $\gamma=\varepsilon_{l}$, then it also holds for $\gamma=L\varepsilon_{l}$ for any $L\geq 1$. This observation implies that

with probability $\geq 1-4e^{-c_{2}R^{2}l\varepsilon_{l}^{2}}$, hence

with the same probability. Choose $R=R(C)$ large enough so that $c_{1}mR^{2}\geq 3m+C$. Putting the bounds for the numerator and denominator of (A.11) together, we get that with probability $\geq 1-\frac{1}{l\varepsilon_{l}^{2}}-4e^{-c_{2}R^{2}l\varepsilon_{l}^{2}}$,

A.5 Proof of Proposition 3.13

The proof follows a standard pattern: on the first step, we show that it is enough to consider the posterior obtained from the prior restricted to a large compact set, and then proving the theorem for the prior with compact support. The second part mimics the classical argument exactly (e.g., see Van der Vaart (2000)).

To show that one can restrict the prior to the compact set, it is enough to establish that for $R$ large enough and $\mathcal{E}_{l}:=\{\theta:\|\theta-\theta_{0}\|\geq\frac{R}{\sqrt{l}}\}$,

can be made arbitrarily small. This follows from the inclusion

(due to the assumed inequality between Hellinger and Euclidean distances) and the bounds for the numerator and the denominator of (A.12) established in the proof of Theorem 3.11.

Appendix B Numerical simulation: additional examples and details.

The generative model of p-parafac has two levels. First, prior probabilities of latent classes and parameters are sampled. Discrete random measure $\nu(\cdot)=\sum_{h=1}^{\infty}\nu_{h}\delta_{h}(\cdot)$ is generated using the stick-breaking construction of DP

where $\nu_{h}$ is the prior probability of responders responses belonging to the latent class $h$. The prior probability of a response to $k$th question depending on the latent class $h$

Second, latent variables and parameters are generated, which are specific to responders in the training data. The latent class of $n$th responder

Finally, the response of $n$th responder for question $k$

This generative model in turn implies that

All sampling algorithms were implemented in Matlab. The samplers ran for 10,000 iterations and every fifth sample was collected after a burn-in of 5000. All experiments were performed on Oracle Grid Engine cluster with 2.6GHz 16 core compute nodes. Memory resources were capped at 16GB and 64GB for sampling from subset and overall posteriors, respectively.

B.2 Dirichlet process mixture model sampling for the stochastic approximation

As we have mentioned, using the stochastic approximation in place of the subset posterior does not lead to increased sample complexity in many cases. Many Bayesian models involve hierarchical exponential family specifications, in which case conditional distributions in Gibbs sampling or acceptance probabilities in Metropolis-Hastings algorithms can be trivially modified to account for the weighted likelihood. Here, we present one particular example of the Dirichlet process mixture model. We augment the original sampling model with latent variables and raise the complete data likelihood to an appropriate power. We have observed excellent performance with this approach in other contexts, including the matrix completion example.

and $i=1,\ldots,n$; $\sigma^{2}$ and $\alpha$ are respectively assigned Inverse-Gamma($a_{0},b_{0}$) and Gamma($a_{\alpha}$, $b_{\alpha}$) priors. Assume that data are partitioned into $m$ subsets of size $l$ such that data on subset $j$ are $\operatorname{\mathcal{X}}_{j}=\{X_{j1},\ldots,X_{jl}\}$. The complete data likelihood for subset $j$ after stochastic approximation is

where $1(\cdot)$ is an indicator function, $K^{*}$ is the maximum number of atoms in the stick breaking representation for $G$, and $l_{j}^{m}(\theta_{1},\ldots,\theta_{K^{*}})$ also depends on latent variables and $\sigma$. Full conditionals of latent variables and unknown parameters are tractable in terms of standard distributions. The Gibbs sampler iterates between the following five steps:

Sample $\theta_{h}\mid\text{rest}$ from Normal($\mu_{h}$, $\sigma_{h}^{2}$) for $h=1,\ldots,K^{*}$, where

Sample $Z_{ji}\mid\text{rest}$ for $i=1,\ldots,l$ from the categorical distribution

Sample $\sigma^{2}\mid\text{rest}$ from Inverse-Gamma($a_{\sigma},b_{\sigma}$), where

Sample $V_{h}\mid\text{rest}$ from Beta($1+\sum_{i=1}^{l}1(z_{ji}=h)$, $\alpha+\sum_{i=1}^{l}1(z_{ji}>h)$) for $h=1,\ldots,K^{*}$.

Sample $\alpha\mid\text{rest}$ from Gamma($a_{\alpha}+K^{*}$, $b_{\alpha}-\sum_{h=1}^{K^{*}}\log(1-V_{h})$).

We fix $\sigma_{\theta}=100,a_{\alpha}=b_{\alpha}=1,$ and $a_{0}=b_{0}=0.01$ following standard conventions.

B.3 Selection of the optimal number of subsets m𝑚m

The following heuristic approach picks the median among the candidate $M$-posteriors. Namely, start by evaluating the M-Posterior for each $m$ in the range of candidate values $[m_{1},m_{2}]$:

and choose $m_{\ast}\in[m_{1},m_{2}]$ such that

where $\mbox{{\rm med}}_{0}$ is the metric median defined in (2.3) in the section 2 of the paper.