Robust Fusion Methods for Big Data

TL;DR

This paper explores robust fusion techniques to combine estimators from subsamples in Big Data, aiming to improve estimation accuracy when handling large datasets without processing the entire sample.

Contribution

It introduces novel robust fusion methods specifically designed for Big Data scenarios where full data processing is infeasible.

Findings

Effective fusion procedures that enhance estimator robustness.

Improved estimation accuracy in large-scale data settings.

Method demonstrates scalability and robustness in experiments.

Abstract

We address one of the important problems in Big Data, namely how to combine estimators from different subsamples by robust fusion procedures, when we are unable to deal with the whole sample.

Tables

Table 1. Table 1: Percentage of estimator breakdowns for 5000 5000 5000 replications and different values of m 𝑚 m for n = 30000 𝑛 30000 n=30000 ; p 𝑝 p is the proportion of outliers

$m$	$p=0.45$	$p=0.49$	$p=0.495$	$p=0.499$
5	0	0.0020	0.0820	0.3892
10	0	0.0088	0.1564	0.5352
30	0	0.0052	0.1426	0.5186
50	0	0.0080	0.1598	0.5412
100	0	0.0192	0.2162	0.6084
150	0	0.0278	0.2728	0.6780

Table 2. Table 2: Covariance operator estimator in presence of outliers. Using the classical and robust estimators over the entire sample, and aggregating by average or fusion of M subsamples estimates. p=0.15, T=20

n	M	time0	time1	Cov	AvCov	Cov.RFM	CovRob	AvRob	RFM
$0.05$	$20$	$553$	$18.20$	$24.3$	$24.3$	$24.7$	$5.16$	$5.21$	$5.52$
$0.05$	$50$	$543$	$7.81$	$24.3$	$24.3$	$24.9$	$5.20$	$5.24$	$5.60$
$0.05$	$100$	$528$	$4.79$	$24.3$	$24.3$	$25.2$	$5.20$	$5.17$	$5.58$
$0.05$	$1000$	$459$	$19.40$	$24.3$	$24.3$	$27.0$	$5.13$	$5.54$	$6.58$
$0.10$	$20$	$2300$	$69.00$	$24.2$	$24.2$	$24.4$	$5.14$	$5.22$	$5.43$
$0.10$	$50$	$2300$	$28.10$	$24.2$	$24.2$	$24.6$	$5.04$	$5.09$	$5.13$
$0.10$	$100$	$2290$	$15.20$	$24.2$	$24.2$	$25.0$	$5.06$	$5.15$	$5.43$
$0.10$	$1000$	$1850$	$21.60$	$24.3$	$24.2$	$26.1$	$5.21$	$5.35$	$6.13$

Table 3. Table 3: Covariance operator estimator in the presence of outliers. Using classical and robust estimators over the entire sample, and aggregating by average or fusion of m 𝑚 m subsamples estimates. p = 0.2 𝑝 0.2 p=0.2 , T=20

n	$m$	time0	time1	Cov	AvCov	Cov.RFM	CovRob	AvRob	RFM
$0.05$	$20$	$572$	$17.90$	$30.5$	$30.5$	$30.9$	$0.879$	$3.96$	$1.45$
$0.05$	$50$	$649$	$7.88$	$30.5$	$30.5$	$31.3$	$0.876$	$7.34$	$2.10$
$0.05$	$100$	$633$	$4.61$	$30.5$	$30.5$	$31.6$	$0.839$	$8.86$	$2.43$
$0.05$	$1000$	$478$	$19.50$	$30.5$	$30.5$	$32.3$	$0.864$	$13.10$	$7.08$
$0.10$	$20$	$1970$	$69.10$	$30.4$	$30.4$	$30.6$	$0.914$	$3.83$	$1.36$
$0.10$	$50$	$2030$	$28.10$	$30.4$	$30.4$	$31.1$	$0.921$	$4.32$	$1.55$
$0.10$	$100$	$2020$	$15.10$	$30.4$	$30.4$	$31.3$	$0.840$	$8.44$	$2.35$
$0.10$	$1000$	$1840$	$21.60$	$30.4$	$30.4$	$32.9$	$0.961$	$12.10$	$5.20$

Equations

$$g_Y(y)=\frac{(2k+1)!}{(k!{)}^2}{F}_X(t{)}^k(1-F_X(t){)}^k{f}_X(t).$$

$$h_Y(y)=\frac{(2k+1)!}{(k!{)}^2}{t}^k(1-t){)}^k{1}_{[0,1]}(t),$$

$$D(x,P)=1-\Vert {\mathbb{E}}_P\left(\frac{X-x}{\Vert X-x\Vert }\right)\Vert ,$$

$${\mathcal{X}}_p=\{\mathbf{y}:card(\mathbf{y})=n,card(\mathbf{x}\cap \mathbf{y})=n-p\},$$

$$\Vert S{\Vert }_{HS}^2=\sum _{k=1}^{\infty }\Vert S(e_k){\Vert }^2<\infty ,⟨S_1,S_2{⟩}_{\mathcal{F}}=\sum _{k=1}^{\infty }⟨S_1(e_k),S_2(e_k)⟩,$$

$$W_i:E\to E,W_i(f)=⟨X_i,f⟩X_i(.),i=1,\dots n.$$

$$d_{ij}^2:=\Vert W_i-W_j{\Vert }_{HS}^2=\Vert X_i{\Vert }^4+\Vert X_j{\Vert }^4-2⟨X_i,X_j{⟩}^2\text{ for }1\le i\le j\le n.$$

$$X(t)=\mu (t)+\sqrt{2}\sum _{k=1}^{10}{\lambda }_k{a}_k\sin (2\pi kt)+\sqrt{2}\sum _{k=1}^{10}{\nu }_k{b}_k\cos (2\pi kt)$$

Full text

1 Introduction

Big Data covers a large list of different problems, see for instance Wang et al. (2016), Yu (2014), and references therein. We address one of them, namely how to combine, using robust techniques, estimators obtained from different subsamples in the case where we are unable to deal with the whole sample. In what follows we will refer to this as Robust Fusion Methods (RFM).

To fix ideas, we start by describing one of the simplest problems in this area as a toy example. Suppose we are interested in the median of a huge set of iid random variables $\{X_{1},\ldots,X_{n}\}$ with common density $f$, and we split the sample into $m$ subsamples of length $l$, $n=ml$. We calculate the median of each subsample and obtain $m$ random variables $Y_{1},\ldots,Y_{m}$. Then we take the median of the set $Y_{1},\ldots,Y_{m}$, i.e. we consider the well known median of medians. It is clear that it does not coincide with the median of the whole original sample $\{X_{1},\ldots,X_{n}\}$, but it will be close. What else could we say about this estimator regarding efficiency and robustness?

The estimator is nothing but the median of $m$ iid random variables but now with a different distribution given by the distribution of the median of $l$ random variables with density $f_{X}$. Suppose for simplicity that $l=2k+1$. Then, the density of the random variables $Y_{i}$ is given by

[TABLE]

In particular, if $f_{X}$ is uniform on $(0,1)$, (1) is given by

[TABLE]

a $Beta(k+1,k+1)$ distribution.

On the other hand, we have that asymptotically, for a sample of size $n$ the empirical median $\theta=med(X_{1},\ldots,X_{n})$ behaves as a normal distribution centred at the true median $\theta$ with variance $\frac{1}{4nf_{X}(\theta)^{2}}$, while the median of medians behaves asymptotically as a normal distribution centred at $\theta$, the median of the median distribution, and with variance $\frac{1}{4mg_{Y}(\theta)^{2}}$. For the uniform case, both are centred at 1/2, $f_{X}(0.5)=1$ and $g_{Y}(0.5)=(1/2)^{2k}(2k+1)!/(k!)^{2}$, so we can explicitly calculate the relative efficiency.

In section 2 we generalize this idea and study the breakdown point, efficiency, and computational time of the robust fusion method. In section 3 we tackle, as a particular case, the robust estimation of the covariance operator and show, in section 3.3, the performance throughout a simulation study.

2 A general setup for RFM.

In this short note we present briefly a general framework for RFM methods for several multivariate and functional data problems. We illustrate our procedure considering only the problem of robust covariance operator estimation, based on a new simple robust estimator. Our approach is to consider RFM methods based on data depth functions. The idea is quite simple: given a statistical problem, (such as multivariate location, covariance operators, linear regression, principal components, among many others), we first split the sample into subsamples. For each subsample we calculate a robust estimator for the statistical problem considered. We will use them all to obtain an RFM estimator that is more accurate. More precisely, the RFM estimator is defined as the deepest point (with respect to the appropriate norm associated to the problem) among all the estimators obtained from the subsamples. Since we need to be able to calculate depths for large sample sizes and eventually high dimensional and infinite dimensional data, we will consider the spatial median corresponding to maximizing the spatial depth function

[TABLE]

where $P$ is a probability in some Banach space $(E,\|\cdot\|)$ and $x\in E$, introduced by Chaudhuri (1996), formulated (in a different way) by Vardi and Zhand (2000), and extended to a very general setup by Chakraborty and Chaudhuri (2014). We want to address the consistency, efficiency, robustness and computational time properties of the RFM proposals.

To be more precise, the general algorithm is as follows. a) We observe $X_{1},\ldots X_{n}$ iid random elements in a metric space $E$ (for instance $E=\mathbb{R}^{d}$), b) we split the sample into subsamples $\{X_{1},\ldots X_{l}\},\{X_{l+1},\ldots X_{2l}\},\ldots,\{X_{(m-1)l+1},\ldots X_{lm}\}$ with $n=ml$, c) we solve our statistical problem on each subsample with a robust procedure (for example, estimate a parameter $\theta$ on each subsample, obtaining $\hat{\theta}_{1},\ldots,\hat{\theta}_{m}$), d) we take the fusion of the results at each subsample, (for instance $\hat{\hat{\theta}}$ can be the deepest point among $\hat{\theta}_{1},\ldots,\hat{\theta}_{m}$.

2.1 Breakdown point

Breakdown point for the RFM. Following Donoho (1982) we consider the finite-sample breakdown point.

Definition 1

Let $\hat{\theta}_{n}=\hat{\theta}_{n}(x)$ be an estimate of $\theta$ defined for samples $\mathbf{x}=\{x_{1},\dots,x_{n}\}$. Let us assume that $\theta$ takes values in $\Theta\subset\mathbb{R}^{d}$ (it can be $\Theta=\mathbb{R}^{d}$). Let $\mathcal{X}_{p}$ be the set of all data sets $\mathbf{y}$ of size $n$ having $n-p$ elements in common with $x$:

[TABLE]

then $\epsilon_{n}^{*}(\hat{\theta}_{n},\mathbf{x})=\frac{p^{*}}{n},$ where $p^{*}=\max\{p\geq 0:\hat{\theta}_{n}(\mathbf{y})$ is bounded and also bounded away from $\partial\Theta\ \forall\mathbf{y}\in\mathcal{X}_{p}\}.$

Let us consider, for $n=ml$, the random walk $S_{n}$ with $S_{0}=0$, and $S_{j}=B_{1}+\ldots+B_{j}$ for $j=1,\dots,n$, $B_{i}$ being iid $\text{Bernoulli}(p)$ for $i=1,\dots,n$, where a one represents the presence of an outlier, while a zero represents no presence of an outlier. Then to compute the breakdown point for the median of medians, we need to count how many times the sequence $\{S_{l},S_{2l}-S_{l},\ldots S_{n}-S_{n-l}\}$ is larger than $k$ (recall that $l=2k+1$). Let us define, $U_{m,n}:=card\{1\leq j\leq m:S_{jl}-S_{(j-1)l}\leq k\}/m,$ since the median has breakdown point $0.5$ the fusion will break down if $U_{m,n}$ is greater than $0.5$.

This will also be true if we take the median of any robust estimate with breakdown point equal to 0.5 calculated at each subsample.

To have a glance at the breakdown point, we performed 5000 replicates of the vector $S_{30000}$ and calculated the percentage of times the estimator breaks down for $p=0.45,0.49,0.495$ and $0.499$. The results are in the following table.

Since the number $Y$ of outliers in the subsamples of length $l$ follows a Binomial distribution, $\text{Binom}(l,p)$, as a direct application of Theorem 1 in Short (2013) we can bound the probability, $q=\mathbb{P}(Y>l/2)$, of breakdown.

2.2 Efficiency of Fusion of $M$-estimators

In this section we will obtain the asymptotic variance of the RFM method, for the special case of $M$-estimators. Recall that $M$-estimators can be defined (see Section 3.2 in Huber and Ronchetti (2009)) by the implicit functional equation $\int\psi(x,T(F))F(dx)=0,$ where $\psi(x;\theta)=(\partial/\partial\theta)\rho(x;\theta)$, for some function $\rho$. The estimator $T_{n}$ is given by the empirical version of $T$, based on a sample $\mathcal{X}_{n}=\{X_{1},\dots,X_{n}\}$. It is well known that $\sqrt{n}(T_{n}-T(F))$ is asymptotically normal with mean 0, variance $\sigma^{2}$, and can be calculated in general as the integral of the square of the influence function. The asymptotic efficiency of $T_{n}$ is defined as ${\rm Eff}(T_{n})=\frac{\sigma_{ML}^{2}}{\sigma^{2}}$, where $\sigma_{ML}^{2}$ is the asymptotic variance of the maximum likelihood estimator. Then, the asymptotic variance of a $M$-estimator built from a sample $T_{n}^{1},\dots,T_{n}^{p}$ of $p$ $M$-estimators of $T$ can be calculated easily.

2.3 Computational time

We want to calculate the computational time of our robust fusion method for a sample $\mathcal{X}_{n}=\{X_{1},\dots,X_{n}\}$ iid of $X$, where we have split $\mathcal{X}_{n}$ into $m$ subsamples of length $l$, then apply a robust estimator to every subsample of length $l$, and fuse them by taking the deepest point among the $m$ subsamples. If we denote by compRE($l$) the computational time required to calculate the robust estimator based on every subsample of length $l$, and compDeph($m$) the computational time required to compute the deepest robust estimator based on the $m$ estimators, then the computational time of our robust fusion method is $m\times$compRE($l$)+compDeph($m$).

3 Robust Fusion for covariance operator

The estimation of the covariance operator of a stochastic process is a very important topic that helps to understand the fluctuations of the random element, as well as to derive the principal functional components from its spectrum. Several robust and non-robust estimators have been proposed, see for instance Chakraborty and Chaudhuri (2014), and the references therein. In order to perform RFM, we introduce a computationally simple robust estimator to apply to each of the $m$ subsamples, that can be performed using parallel computing. It is based on the notion of impartial trimming applied on the Hilbert–Schmidt space, where covariance operators are defined. The RFM estimator is the deepest point among the $m$ estimators corresponding to each subsample, where the norm in (3) is given by (4) below.

3.1 A resistant estimate of the covariance operator

Let $E=L^{2}(T)$, where $T$ is a finite interval in $\mathbb{R}$, and $X,X_{1},\ldots X_{n},\ldots$ iid random elements taking values in $(E,\mathcal{B}(E))$, where $\mathcal{B}(E)$ stands for the Borel $\sigma$-algebra on $E$. Assume that $\mathbb{E}(X(t)^{2})<\infty$ for all $t\in T$, and $\int_{T}\int_{T}\rho^{2}(s,t)dsdt<\infty$, so the covariance function is well defined and given by $\rho(s,t)=\mathbb{E}((X(t)-\mu(t))(X(s)-\mu(s))),\ \ \mbox{where}\ \mathbb{E}(X(t))=\mu(t)$.

For notational simplicity we assume that $\mu(t)=0,\forall t\in T$. Under these conditions, the covariance operator, given by $S_{0}(f)=\mathbb{E}(\langle X,f\rangle X)$, is diagonalizable, with nonnegative eigenvalues $\lambda_{i}$ such that $\sum_{i}\lambda^{2}_{i}<\infty$. Moreover $S_{0}$ belongs to the Hilbert–Schmidt space $HS(E,E)$ of linear operators with square norm and inner product given by

[TABLE]

respectively, where $\{e_{k}:k\geq 1\}$ is any orthonormal basis of $E$, and $S,S_{1},S_{2}\in HS(E,E)$. In particular, $\|S_{0}\|^{2}=\sum_{i=1}^{\infty}\lambda_{i}^{2}$, where $\lambda_{i}$ are the eigenvalues of $S_{0}$. Given an iid sample $X_{1},\ldots,X_{n}$, define the Hilbert–Schmidt operators of rank one,

[TABLE]

Let $\phi_{i}=X_{i}/\|X_{i}\|$. Then, $W_{i}(\phi_{i})=\|X_{i}\|^{2}\phi_{i}=:\eta_{i}\phi_{i}.$

The standard estimator of $S_{0}$ is just the average of these operators, i.e. $\hat{S}_{n}=\frac{1}{n}\sum_{i=1}^{n}W_{i}$, which is a consistent estimator of $S_{0}$ by the Law of Large Numbers in the space $HS(E,E)$. Our proposal is to consider an impartial trimmed estimate as a resistant estimator. The notion of impartial trimming was introduced in Gordaliza (1991), while the functional data setting was considered in Cuesta-Albertos and Fraiman (2006), from were it can be obtained the asymptotic theory for our setting. In order to perform the algorithm we will derive an exact formula for the matrix distances $\|W_{i}-W_{j}\|$, $1\leq i\leq j\leq n$.

**Lemma 1 **We have that

[TABLE]

3.2 The impartial trimmed mean estimator

Following Cuesta-Albertos and Fraiman (2006), we define the impartial trimmed covariance operator estimator, which is calculated by the following algorithm.

Given the sample $X_{1}(t),\ldots,X_{n}(t)$ (which we have assumed with mean zero for notational simplicity) and $0<\alpha<1$, we provide a simple algorithm to calculate an approximate impartial trimmed mean estimator of the covariance operator of the process $S_{0}:E\rightarrow E,\hskip 5.69054ptS_{0}(f)(t)=\mathbb{E}(\langle X,f\rangle X(t))$, that will be strongly consistent.

STEP 1: Calculate $d_{ij}=\|W_{i}-W_{j}\|_{HS}$, $1\leq i\leq j\leq n$, using Lemma 1.

STEP 2: Let $r=\lfloor(1-\alpha)n\rfloor+1$. For each $i=1,\ldots n$, consider the set of indices $I_{i}\subset\{1,\ldots,n\}$ corresponding to the $r$ nearest neighbours of $W_{i}$ among $\{W_{1},\ldots W_{n}\}$, and the order statistic of the vector $(d_{i1},\ldots,d_{in})$, $d_{i}^{(1)}<\ldots<d_{i}^{(n)}$.

STEP 3: Let $\gamma=\text{argmin}\{d_{1}^{(r)},\ldots,d_{n}^{(r)}\}$.

STEP 4: The impartial trimmed mean estimator of $S_{0}$ is given by $\hat{S}=$ the average of the $m$ nearest neighbours of $W_{\gamma}$ among $\{W_{1},\ldots,W_{n}\}$, i.e the average of the observations in $I_{\gamma}$.

This estimator corresponds to estimating $\rho(s,t)$ by $\hat{\rho}(s,t)=\frac{1}{r}\sum_{j\in I_{\gamma}}X_{j}(s)X_{j}(t)$. Observe that steps 1 and 2 of the algorithm can be performed using parallel computing.

3.3 Simulation results for the covariance operator

Simulations were done using a PC Intel core i7-3770 CPU, 8GO of RAM using 64 bit version of Win10, and R software ver. 3.3.0.

We vary the sample size $n$ within the set $\{0.1e6,1e6,5e6,10e6\}$ and the number of subsamples $m\in\{100,500,1000,10000\}$. The proportion of outliers was fixed to $p=13\%$ and $p=15\%$. We replicate each simulation case $K=5$ times and report a mean average of the results over these replicates.

We report the average time in seconds necessary for both a global estimate (time0, over the whole sample), and time1 the estimate obtained by fusion (including computing the estimates over subsamples and aggregating them by fusion).

We compare the classical estimator (Cov), the mean of the classical estimators obtained from the subsamples (AvCov), the Fusion estimate of the classical estimator (Cov.RFM), the global robust estimate (CovRob), the robust fusion estimate RFM, and the average of the robust estimates from the subsamples AvRob.

To generate the data, we have used a simplified version of the simulation model used in Kraus and Panaretos (2012):

[TABLE]

where $\nu_{k}=\left(\frac{1}{3}\right)^{k},\lambda_{k}=k^{-3}$, and $a_{k}$ and $b_{k}$ are random standard Gaussian independent observations. The central observations were generated using $\mu(t)=0$ whereas for the outliers we took $\mu(t)=2-8\sin(\pi t)$. For $t$ we used an equally spaced grid of $T=20$ points in $[0,1]$.

The covariance operator of this process was computed for the comparisons:

$Cov(s,t)=\sum_{k=1}^{10}A_{k}(s)A_{k}(t)+B_{k}(s)B_{k}(t)$, where $A_{k}(t)=\sqrt{2}\lambda_{k}\sin(2\pi kt)$ and $A_{k}(t)=\sqrt{2}\nu_{k}\cos(2\pi kt)$.

The results are shown in the following two tables for two proportions of outliers, $p=0.15$ and $p=0.2$.

If the proportion of outliers is moderate $p=15\%$, the average of the robust estimators still behaves well, better than RFM, but if we increase the proportion of outliers to $p=0.2$, RFM clearly outperforms all the other estimators.

Acknowledgement

We thanks an anonymous referee for helpful suggestions on a first version.