Shibin Parameswaran

Charles-Alban Deledalle
email: cdeledalle@eng.ucsd.edu

Loïc Denis
email: loic.denis@univ-st-etienne.fr.

Truong Q Nguyen

C.-A Deledalle

C Deledalle

Accelerating GMM-based patch priors for image restoration: Three ingredients for a 100× speed-up

Keywords: Image restoration, image patch, Gaussian mixture model, efficient algorithms I

Image restoration methods aim to recover the underlying clean image from corrupted observations. The Expected Patch Log-likelihood (EPLL) algorithm is a powerful image restoration method that uses a Gaussian mixture model (GMM) prior on the patches of natural images. Although it is very effective for restoring images, its high runtime complexity makes EPLL ill-suited for most practical applications. In this paper, we propose three approximations to the original EPLL algorithm. The resulting algorithm, which we call the fast-EPLL (FEPLL), attains a dramatic speed-up of two orders of magnitude over EPLL while incurring a negligible drop in the restored image quality (less than 0.5 dB). We demonstrate the efficacy and versatility of our algorithm on a number of inverse problems such as denoising, deblurring, super-resolution, inpainting and devignetting. To the best of our knowledge, FEPLL is the first algorithm that can competitively restore a 512×512 pixel image in under 0.5s for all the degradations mentioned above without specialized code optimizations such as CPU parallelization or GPU implementation.

Mixture Models (GMMs) to capture the prior distribution of patches in natural images. Note that a similar idea was introduced concurrently in [START_REF] Yu | Solving inverse problems with piecewise linear estimators: From gaussian mixture models to structured sparsity[END_REF]. The success of this method is evident from the large number of recent works that extend the original EPLL formulation [START_REF] Guillemot | Covariance Trees for 2D and 3D Processing[END_REF], [START_REF] Sulam | Expected patch log likelihood with a sparse prior[END_REF], [START_REF] Luo | Adaptive image denoising by mixture adaptation[END_REF], [START_REF] Cai | Image denoising via patch-based adaptive gaussian mixture prior method[END_REF], [START_REF] Papyan | Multi-scale patch-based image restoration[END_REF], [START_REF] Ren | Example-based image synthesis via randomized patch-matching[END_REF], [START_REF] Houdard | High-Dimensional Mixture Models For Unsupervised Image Denoising (HDMI)[END_REF]. However, a persistent problem of EPLL-based algorithms is their high runtime complexity. For instance, it is orders of magnitude slower than the well-engineered BM3D image denoising algorithm [8]. However, extensions of BM3D that perform super-resolution [START_REF] Danielyan | Spatially adaptive filtering as regularization in inverse imaging: Compressive sensing superresolution and upsampling[END_REF] and other inverse problems [START_REF] Katkovnik | Nonlocal image deblurring: Variational formulation with nonlocal collaborative L0-norm prior[END_REF] require fundamental algorithmic changes, making BM3D far less adaptable than EPLL. Other approaches that are as versatile as EPLL [START_REF] Schmidt | Shrinkage fields for effective image restoration[END_REF], [START_REF] Chen | Insights into analysis operator learning: From patch-based sparse models to higher order mrfs[END_REF], [START_REF] Jancsary | Regression tree fields -an efficient, non-parametric approach to image labeling problems[END_REF] either lack the algorithmic efficiency of BM3D or the restoration efficacy of EPLL.

Another class of techniques that arguably offers better runtime performance than EPLL-based methods (but not BM3D) are those based on deep learning. With the advancements in computational resources, researchers have attempted to solve some classical inverse problems using multi-layer perceptrons [3] and deep networks [START_REF] Chen | On learning optimized reaction diffusion processes for effective image restoration[END_REF], [START_REF] Dong | Image super-resolution using deep convolutional networks[END_REF], [START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF]. These methods achieve very good restoration performance, but are heavily dependent on the amount of training data available for each degradation scenario. Most of these methods learn filters that are suited to restore a specific noise level (denoising), blur (deblurring) or upsampling factor (super-resolution), which makes them less attractive to serve as generic image restoration solutions. More recently, Zhang et al. [START_REF] Zhang | Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising[END_REF] demonstrated the use of deep residual networks for general denoising problems, single-image superresolution and compression artifact removal. Unlike earlier deep learning efforts, their approach can restore images with different noise levels using a single model which is learned by training on image patches containing a range of degradations. Even in this case, the underlying deep learning model requires retraining whenever a new degradation scenario different from those considered during the learning stage is encountered.

More recently, [START_REF] Chang | One network to solve them all-solving linear inverse problems using deep projection models[END_REF] proposed training a single deep network to solve many inverse problems. This work uses an iterative scheme that alternatively enforces a good fit to the learned prior model of natural images (via a projection operator performed by a deep neural network) and a satisfying fidelity to the data (via the direct model). The projection operator takes the form of an auto-encoder, trained by an adversarial strategy, that processes 64 × 64 patches with 1024-dimensional latent space (i.e., 4 times smaller). This framework shows promise but, in its current form, requires a very large image database for training. In addition, almost all of the training and testing in the original publication are conducted on small 64 × 64 images. Another limitation, as noted by the authors, is that the Because of ArXiv file size limitation, images have been compressed in this document. Assessing restoration quality by visiual inspection is then impossible due to compression artifacts. Please refer to https://www.math.u-bordeaux.fr/ ~cdeledal/fepll to access a pdf file with optimal quality. 1 regularization parameter that controls the contribution of prior is fixed during training stage. In other words, regularization parameter cannot be changed during test time which can be an issue in certain situations. In contrast, our approach needs much less data during training stage and, during test time, the regularization parameter can be tuned according to the signal-to-noise ratio of the image being addressed.

Moreover, it is much harder to gain insight into the actual model learned by a deep architecture compared to a GMM. For this reason, even with the advent of deep learning methods, flexible algorithms like EPLL that have a transparent formulation remain relevant for image restoration.

Recently, researchers have tried to improve the speed of EPLL by replacing the most time-consuming operation in the EPLL algorithm with a machine learning-based technique of their choice [START_REF] Wang | Discriminative indexing for probabilistic image patch priors[END_REF], [START_REF] Rosenbaum | The return of the gating network: combining generative models and discriminative training in natural image priors[END_REF]. These methods were successful in accelerating EPLL to an extent but did not consider tackling all of its bottlenecks. In contrast, this paper focuses on accelerating EPLL by proposing algorithmic approximations to all the prospective bottlenecks present in the original algorithm proposed by Zoran et al. [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]. To this end, we first provide a complete computational and runtime analysis of EPLL, present a new and efficient implementation of original EPLL algorithm and then finally propose innovative approximations that lead to a novel algorithm that is more than 100× faster compared to the efficiently implemented EPLL (and 350× faster than the runtime obtained by using the original implementation [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]).

Contributions: The main contributions of this work are the following. We introduce three strategies to accelerate patchbased image restoration algorithms that use a GMM prior. We show that, when used jointly, they lead to a speed-up of the EPLL algorithm by two orders of magnitude. Compared to the popular BM3D algorithm, which represents the current state-ofthe-art in terms of speed among CPU-based implementations, the proposed algorithm is almost an order of magnitude faster. The three strategies introduced in this work are general enough to be applied individually or in any combination to accelerate other related algorithms. For example, the random subsampling strategy is a general technique that could be reused in any algorithm that considers overlapping patches to process images; the flat tail spectrum approximation can accelerate any method that needs Gaussian log-likelihood or multiple Mahalanobis metric calculations; finally, the binary search tree for Gaussian matching can be included in any algorithm based on a GMM prior model and can be easily adapted for vector quantization techniques that use a dictionary.

For reproducibility purposes, we release our software on GitHub along with a few usage demonstrations (available at https://goo.gl/xjqKUA).

II. EXPECTED PATCH LOG-LIKELIHOOD (EPLL)

We consider the problem of estimating an image x ∈ R N (N is the number of pixels) from noisy linear observations y = Ax + w, where A : R N → R M is a linear operator and w ∈ R M is a noise component assumed to be white and Gaussian with variance σ 2 . In a standard denoising problem A is the identity matrix, but in more general settings, it can account for loss of information or blurring. Typical examples for operator A are: a low pass filter (for deconvolution), a masking operator (for inpainting), or a projection on a random subspace (for compressive sensing). To reduce noise and stabilize the inversion of A, some prior information is used for the estimation of x. The EPLL introduced by Zoran and Weiss [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] includes this prior information as a model for the distribution of patches found in natural images. Specifically, the EPLL defines the restored image as the maximum a posteriori estimate, corresponding to the following minimization problem:

argmin x P 2σ 2 ||Ax -y|| 2 - i∈I log p (P i x) (1)
where I = {1, . . . , N } is the set of pixel indices, P i : R N → R P is the linear operator extracting a patch with P pixels centered at the pixel with location i (typically, P = 8 × 8), and p(.) is the a priori probability density function (i.e., the statistical model of noiseless patches in natural images).

While the first term in eq. ( 1) ensures that Ax is close to the observations y (this term is the negative log-likelihood under the white Gaussian noise assumption), the second term regularizes the solution x by favoring an image such that all its patches fit the prior model of patches in natural images. The authors of [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] showed that this prior can be well approximated (upon removal of the DC component of each patch) using a zero-mean Gaussian Mixture Model (GMM) with K = 200 components, that reads for any patch z ∈ R P , as

p (z) = K k=1 w k 1 (2π) P/2 |Σ k | 1/2 exp - 1 2 z t Σ -1 k z , (2)
where the weights w k (such that w k > 0 and k w k = 1) and the covariance matrices Σ k ∈ R P ×P are estimated using the Expectation-Maximization algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] on a dataset consisting of 2 million "clean" patches extracted from the training set of the Berkeley Segmentation (BSDS) dataset [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF].

Half-quadratic splitting: Problem (1) is a large non-convex problem where A couples all unknown pixel values x and the patch prior is highly non-convex. A classical workaround, known as half-quadratic splitting [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF], [START_REF] Krishnan | Fast image deconvolution using hyperlaplacian priors[END_REF], is to introduce N auxiliary unknown vectors z i ∈ R P , and consider instead the penalized optimization problem that reads, for β > 0, as

argmin x,z1,...,z N P 2σ 2 ||Ax -y|| 2 + β 2 i∈I ||P i x -z i || 2 - i∈I log p (z i ) . (3)
When β → ∞, the problem (3) becomes equivalent to the original problem (1). In practice, an increasing sequence of β is considered, and an alternating optimization scheme is used:

ẑi ← argmin zi β 2 ||P i x -z i || 2 -log p (z i ) i=1..N (4)
x ← argmin Algorithm 1 The five steps of an EPLL iteration [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] for all i ∈ I zi ← Pi x (Patch extraction)

x P 2σ 2 ||Ax -y|| 2 + β 2 i∈I ||P i x -ẑi || 2 . ( 5
)
k i ← argmin 1 k i K log w -2 k i + log Σ k i + 1 β IdP + zt i Σ k i + 1 β IdP -1 zi (Gaussian selection) ẑi ← Σ k i + 1 β IdP -1 Σ k i zi (Patch estimation) x ← i∈I P t i Pi -1 i∈I P t i ẑi (Patch reprojection) x ← A t A + βσ 2 IdN -1 A t y + βσ 2 x (Image estimation)
Algorithm: Subproblem (5) corresponds to solving a linear inverse problem with a Tikhonov regularization, and has an explicit solution often referred to as Wiener filtering:

x = A t A + βσ 2 P i∈I P t i P i -1 A t y + βσ 2 P i∈I P t i ẑi , (6)
where P t i P i is a diagonal matrix whose i-th diagonal element corresponds to the number of patches overlapping the pixel of index i. This number is a constant equal to P (assuming proper boundary conditions are used), which allows to split the computation into two steps Patch reprojection and Image estimation as shown in Alg. 1. Note that the step Patch reprojection is simply the average of all overlapping patches. In contrast, subproblem (4) cannot be obtained in closed form as it involves a term with the logarithm of a sum of exponentials. A practical solution proposed in [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] is to keep only the component k i maximizing the likelihood for the given patch assuming it is a zero-mean Gaussian random vector with covariance matrix Σ ki + 1 β Id P . With this approximation, the solution of (4) is also given by Wiener filtering, and the resulting algorithm iterates the steps described in Alg. 1. The authors of [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] found that using T = 5 iterations, with the sequence β = 1 σ 2 {1, 4, 8, 16, 32}, for the initialization x = y, provides relevant solutions in denoising contexts for a wide range of noise level values σ 2 .

A. Complexity via eigenspace implementation

The algorithm summarized in the previous section may reveal cumbersome computations as it requires performing numerous matrix multiplications and inversions. Nevertheless, as the matrices Σ k are known prior to any calculation, their eigendecomposition can be computed offline to improve the runtime. If we denote the eigendecomposition (obtained offline) of

Σ k = U k S k U t k , such that U k ∈ R P ×P is
unitary and S k is diagonal with positive diagonal elements ordered in decreasing order, steps Gaussian selection and Patch estimation can be expressed in the space of coefficients c as

ck i ← U t k zi k=1..K i=1..N O(N KP 2 ) ( 7
)
k i ← argmin 1 k K ι k + P j=1 log ν k j + [c k i ] 2 j ν k j i=1..N O(N KP ) (8) [ĉ i ] j ← γ k i j [c k i i ] j j=1..P i=1..N O(N P ) (9) ẑi ← U k i ĉi i=1..N O(N P 2 ) (10)
where [c] j denotes the j-th entry of vector c,

ι k = -2 log w k , ν k j = [S k ] jj + 1 β , and γ k j = [S k ] jj /ν k j with [S k
] jj the jth entry on the diagonal of matrix S k . The complexity of each operation is indicated on its right and corresponds to the number of operations per iteration of the alternate optimization scheme. The steps Patch extraction and Patch reprojection share a complexity of O(N P ). Finally, the complexity of step Image estimation depends on the transform A. In many scenarios of interest, A t A can be diagonalized using a fast transform and the inversion of A t A+βσ 2 Id N can be performed efficiently in the transformed domain (since Id N is diagonal in any orthonormal basis). For instance, it leads to O(N ) operations for denoising or inpainting, and O(N log N ) for periodical deconvolutions or super-resolution problems, thanks to the fast Fourier transform (these are the settings we have adopted in this paper). If A cannot be easily diagonalized, this step can be performed using conjugate gradient (CG) method, as done in [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF], at a computational cost that depends on the number of CG iterations (i.e., on the conditioning of A t A + βσ 2 Id N ). In any case, as shown in the next section, this step has a complexity independent of P and K and is one of the faster operations in the image restoration problems considered in this paper.

B. Computation time analysis

In order to uncover the practical computational bottlenecks of EPLL, we have performed the following computational analysis. To identify clearly which part is time consuming, it is important to make the algorithm implementation as optimal as possible. Therefore, we refrain from using the MATLAB code provided by the original authors [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] for speed comparisons. Instead, we use a MATLAB/C version of EPLL based on the eigenspace implementation described above, where some steps are written in C language and interfaced using mex functions. This version, which we refer to as EPLLc, provides results identical to the original implementation while being 2-3 times faster. The execution time of each step for a single run of EPLLc is reported in the second column of Table I. Reported times fit our complexity analysis and clearly indicate that the

Step Gaussian selection causes significant bottleneck due to O(N P 2 K) complexity.

In the next section, we propose three independent modifications leading to an algorithm with a complexity of O(N P r log K/s 2 ) with two constants 1 s 2 P and 1 r P that control the accuracy of the approximations introduced. The algorithm, in practice, is more than 100 times faster as shown by its runtimes reported in the third column.

III. FAST EPLL: THE THREE KEY INGREDIENTS

We propose three accelerations based on (i) scanning only a (random) subset of the N patches, (ii) reducing the number of mixture components matched, and (iii) projecting on a smaller subspace of the covariance eigenspace. We begin by describing this latter acceleration strategy in the following paragraph.

A. Speed-up via flat tail spectrum approximation

To avoid computing the P coefficients of the vector ck i in eq. ( 7), we rely on a flat-tail approximation. The k-th Gaussian model is said to have a flat tail if there exists a rank r k such that for any j > r k , the eigenvalues are constant:

[S k ] j,j = λ k . Denoting by Ū k ∈ R P ×r k (resp. Ū c k ∈ R P ×r c k ) the matrix formed by the r k first (resp. r c k = P -r k last) columns of U k , we have Ū c k ( Ū c k ) t = Id P -Ū k Ū t k . It follows (Σ k + 1 β Id P ) -1 = Ū k ( Sk + 1 β Id r k ) -1 Ū t k + (λ k + 1 β ) -1 (Id P -Ū k Ū t k ), (11)
(Σ k + 1 β Id P ) -1 Σ k = Ū k ( Sk + 1 β Id r k ) -1 Sk Ū t k + λ k (λ k + 1 β ) -1 (Id P -Ū k Ū t k ), (12)
where Sk ∈ R r k ×r k is the diagonal matrix formed by the r k first rows and columns of S k . Steps Gaussian selection and Patch estimation can thus be rewritten as

ck i ← Ū t k zi k=1..K i=1..N O(N KP r) ( 13
)
k i ← argmin 1 k K ι k + r c k log ν k P + ||z i || 2 ν k P + r k j=1 log ν k j + [c k i ] 2 j ν k j - [c k i ] 2 j ν k P i=1..N O(N K r) (14) [ĉ i ] j ← (γ k i j -γ k i P ) ck i i j j=1..r k i i=1..N O(N P r) (15) ẑi ← Ū k i ĉi +γ k i P zi i=1..N O(N P r) (16)
where

ν k P = λ k + 1 β , γ k P = λ k /ν k P .
As ||z i || 2 can be computed once for all k, the complexity of each step is divided by P/r, where r = 1 K K k=1 r k is the average rank after which eigenvalues are considered constant.

In practice, covariance matrices Σ k are not flat-tail but can be approximated by a flat-tail matrix by replacing the lowest eigenvalues by a constant λ k . To obtain a small value of r (hence a large speed-up), we preserve a fixed proportion ρ ∈ (0, 1] of the total variability and replace the smallest eigenvalues accounting for the remaining 1-ρ fraction of the variability by their average (see Fig. 1): r k is the smallest integer such that Tr( Sk ) ρ Tr(S k ). Choosing ρ = 0.95 means that 5% of the variability, in the eigendirections associated to the smallest eigenvalues, is assumed to be evenly spread in these directions. In practice, the choice of ρ = 0.95 leads to an average rank of r = 19.6 (for P = 8 × 8 ) for a small drop of PSNR as shown in Fig. 4. Among several other covariance approximations that we tested, for instance, the one consisting in keeping only the r k first directions, the flat tail approximation provided the best trade-off in terms of acceleration and restoration quality. The analyses showing the superiority of the proposed approximation over the more common approach of keeping only the first r k directions, and the effect of ρ on image quality are included in the supplementary document.

B. Speed-up via a balanced search tree

As shown in Table I, the step Gaussian selection has a complexity of O(N P 2 K), reduced to O(N P K r) using the flat tail spectrum approximation. This step remains the biggest bottleneck since each query patch has to be compared to all the K components of the GMM. To make this step even more efficient, we reduce its complexity using a balanced search tree. As described below, such a tree can be built offline by adopting a bottom-up strategy that repeatedly collapses the original GMM to models with fewer components, until the entire model is reduced to a single Gaussian model.

We progressively combine the GMM components from one level to form the level above, by clustering the K components into L < K clusters of similar ones, until the entire model is reduced to a single component. The similarity between two zero-mean Gaussian models with covariance Σ 1 and Σ 2 is measured by the symmetric Kullback-Leibler (KL) divergence

KL(Σ 1 , Σ 2 ) = 1 2 Tr(Σ -1 2 Σ 1 + Σ -1 1 Σ 2 -2Id P ). (17)
Based on this divergence, at each level n, we look for a partition Ω n of the K Gaussian components into L clusters (with about equal sizes) minimizing the following optimization problem argmin

Ω n L l=1 k1,k2∈Ω n l KL(Σ k1 , Σ k1 ), (18)
such that

L l=1 Ω n l = [K] and Ω n l1 ∩ l1 =l2
Ω n l2 = ∅, where Ω n l is the l-th set of Gaussian components for the GMM at level n. This clustering problem can be approximately solved using the genetic algorithm of [START_REF] Kirk | Source code for "multiple traveling salesmen problemgenetic algorithm[END_REF] for the Multiple Traveling Salesmen Problem (MTSP). MTSP is a variation of the classical Traveling Salesman Problem where several salesmen visit a unique set of cities and return to their origins, and each city is visited by exactly one salesman. This attempts to minimize the total distance traveled by all salesmen. Hence, it is similar to our original problem given in eq. ( 18) where the Gaussian components and the clusters correspond to K cities and L salesmen, respectively. Given the clustering at level n, the new GMM at level n -1 is obtained by combining the zero-mean Gaussian components such that, for all 1 l L:

w n-1 l = k∈Ω n l w n k and Σ n-1 l = 1 w n-1 l k∈Ω n l w n k Σ n k , ( 19
)
where Σ n k and w n k are the corresponding covariance matrix and weight of the k-th Gaussian component at level n. Following this scheme, the original GMM of K = 200 components is collapsed into increasingly more compact GMMs with K = 64, 32, 16, 8, 4, 2 and 1 components. The main advantage of using MTSP compared to classical clustering approaches, is that this procedure can be adapted easily to enforce approximately equal sized clusters, simply by enforcing that each salesman visits at least 3 cities for the last level and 2 for the other ones.

We also experimented with other clustering strategies such as the hierarchical kmeans-like clustering in [START_REF] Goldberger | Hierarchical clustering of a mixture model[END_REF] and hierarchical agglomerative clustering. With no principled way to enforce even-sized clusters, these approaches, in general, lead to unbalanced trees (with comb structured branches) which result in large variations in computation times from one image to another. Although they all lead to similar denoising performances, we opted for MTSP based clustering to build our Gaussian tree in favor of obtaining a stable speed-up profile for our resulting algorithm. Please refer to the supplementary document for timing comparisons of MTSP vs. other tree building strategies.

In Fig. 2 we show that the tree obtained using MTSP-based clustering is almost a binary tree (left) and also display the types of patches it encodes along a given path (right). Such a balanced tree structure lets one avoid testing each patch against all K components. Instead, a patch is first compared to the two first nodes in level 1 of the tree, then the branch providing the smallest cost is followed and the operation is repeated at higher levels until a leaf has been reached. Using this balanced search tree reduces the complexity of step Gaussian selection to O(N P r log K).

C. Speed-up via the restriction to a random subset of patches

The simplest and most effective proposed acceleration consists of subsampling the set I of N patches to improve the complexity of the four most time-consuming steps, see Table I. One approach, followed by BM3D [8], consists of restricting the set I to locations on a regular grid with spacing s ∈ [1,

√ P ] pixels in both directions, leading to a reduction of complexity by a factor s 2 . We refer to this approach as the regular patch subsampling. A direct consequence is that |I| = N/s 2 and the complexity is divided by s 2 . However, we observed that this strategy consistently creates blocky artifacts revealing the regularity of the extraction pattern. A random sampling approach, called "jittering", used in the computer graphics community [START_REF] Cook | Stochastic sampling in computer graphics[END_REF] is preferable to limit this effect. This procedure ensures that each pixel is covered by at least one patch. The location (i 0 , j 0 ) of a point of the grid undergoes a random perturbation, giving a new location (i, j) such that

i 0 - √ P -s 2 i i 0 + √ P -s 2
and j 0 -

√ P -s 2 j j 0 + √ P -s 2 , ( 20
)
where • denotes the flooring operation. We found experimentally that independent and uniform perturbations offered the best performance in terms of PSNR and visual quality against all other tested strategies. In addition, we also resample these positions at each of the T iterations and add a (random) global shift to ensure that all pixels have the same expected number of patches covering them. Figure 3 illustrates the difference between a regular grid and a jittered grid of period s = 6 for patches of size P = 8×8. In both cases, all pixels are covered by at least one patch, but the stochastic version reveals an irregular pattern.

Nevertheless, when using random subsampling, a major bottleneck occurs when A t A is not diagonal because the inversion involved in eq. ( 6) cannot be simplified as in Alg. 1. Using a conjugate gradient is a practical solution but will negate the reduction of complexity gained by using subsampling. To the best of our knowledge, this is the main reason why patch subsampling has not been utilized to speed up EPLL. Here, we follow a different path. We opt for approximating the solution of the original problem (involving all patches) rather than evaluating the solution of an approximate problem (involving random subsample of patches). More precisely, we speed up Alg. 1 by replacing the complete set of indices by the random subset of patches. In this case, step Patch reprojection consists of averaging only this subset of overlapping restored patches. This novel and nuanced idea avoids additional overhead and attains dramatic complexity improvements compared to the standard approach. Note that even in the case of some inverse problems, such as deblurring, super-resolution, inpainting and devignetting, this strategy can still be used in order to avoid conjugate gradient and maintain a large speed-up.

Experiments conducted on our validation dataset show that this strategy used with s = 6 leads to an acceleration of about 36× with less than a 0.2dB drop in PSNR. In comparison, for a similar drop of PSNR, the regular patch subsampling can only achieve a 9× acceleration with s = 3 (plots included in the supplementary).

D. Performance analysis

Figure 4 shows the image restoration performance and speedup obtained when the three ingredients are applied separately or in combination. The results are averaged over 40 images from the test set of BSDS dataset [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF] that is set aside for validation purposes. The speed-up is calculated with respect to the EPLLc implementation which is labeled "original" in Fig. 4. Among the three ingredients, random subsampling or jittering (labeled "subsampling") leads to the largest speed-up (32×), while the usage of the search tree provides more than 7× faster processing. The average speed-up obtained when combining all three ingredients is around 179× on our validation set consisting of images of size 481×321, for an average drop of PSNR less than 0.5dB. Note that most of the runtime gain is achieved by combining binary search tree and random subsampling (166×), while the flat tail acceleration leads to a modest acceleration of less than 10%. Nevertheless, we include it in our approach because it introduces another parameter (ρ) that addresses the bottleneck presented by the dimensionality (P 2 ) in a unique way. Adjusting ρ lets a practitioner choose an appropriate operating point suitable for their need in the speed-vs-quality trade-off space without losing too much in terms of quality.

IV. RELATED METHODS

To the best of our knowledge, there are only two other approaches [START_REF] Wang | Discriminative indexing for probabilistic image patch priors[END_REF], [START_REF] Rosenbaum | The return of the gating network: combining generative models and discriminative training in natural image priors[END_REF] that have attempted to accelerate EPLL. Unlike our approach, these methods focus on accelerating only one of the steps of EPLL namely the Gaussian selection step. Both use machine learning techniques to reduce its runtime.

In [START_REF] Wang | Discriminative indexing for probabilistic image patch priors[END_REF], the authors use a binary decision tree to approximate the mapping zi → k * i performed in step Gaussian selection. At each node k and level n of the tree, the patch zi is confronted with a linear separator in order to decide if the recursion should continue on the left or right child given by

a n k , zi + b n k 0 (21)
where (a n k , b n k ) are the parameters of the hyperplane for the k-th node at level n. These separators are trained offline on all pairs of (z i , k * i ) obtained after the first iteration of EPLL for a given β and noise level σ. Once a leaf has been reached, its index provides a first estimate for the index k * i . To reduce errors due to large variations among the neighboring pixels, this method further employs a Markov random fields on the resulting map of Gaussian components which runs in O(N K) complexity.

Hence, their overall approach reduces the complexity of step Gaussian selection from O(N KP 2 ) to O(N (P D+K)), where D = 12 is the depth of the learned decision tree.

In [START_REF] Rosenbaum | The return of the gating network: combining generative models and discriminative training in natural image priors[END_REF], the authors approximate the Gaussian selection step, by using a gating (feed-forward) network with one hidden layer

zi → ι k + Q j=1 log ν k j +ω k j ( di ) 2 j with di = V t zi (22)
where Q is the size of the hidden layer. The matrix V ∈ R P ×Q encodes the weights of the first layer, ω k corresponds to the weights of the hidden layer and they are learned discriminatively to approximate the exact posterior probability:

zi → ι k + P j=1 log ν k j + (c k i ) 2 j ν k j with ck i = U t k zi (23)
that we encounter in eq. ( 7) and ( 8). Theoretically, a new network will need to be trained for each type of degradations, noise levels and choices of β (recall that ν k j = (S k ) jj + 1 β ). However, the findings of [START_REF] Rosenbaum | The return of the gating network: combining generative models and discriminative training in natural image priors[END_REF] indicate that applying a network learned on clean patches and with 1 β = 0 is effective regardless of the type of degradation or the value of β. Their main advantage can be highlighted by comparing eq. ( 22) and [START_REF] Kim | Accurate image super-resolution using very deep convolutional networks[END_REF] where complexity is reduced from O(N KP 2 ) to O(N Q(K + P )). The authors utilize this benefit by choosing Q = 100.

Unlike these two approaches, our method does not try to learn the Gaussian selection rule directly (which depends on both the noise level through 1/β and the prior model through the GMM). Instead, we simply define a hierarchical organization of the covariance matrices Σ k . In other words, while the two other approaches try to infer the posterior probabilities (or directly the maximum a posteriori), our approach provides an approximation to the prior model. During runtime, this approximation of the prior is used in the posterior for the Gaussian selection task. Please note that the value of β does not play a role in determining the prior. This allows us to use the same search tree independently of the noise level, degradations, etc. Given that the main advantage of EPLL is that the same model can be used for any type of degradations, it is important that this property remains true for the accelerated version. Last but not least, the training of our search tree takes a few minutes while the training steps for the above mentioned approach take from several hours to a few days [START_REF] Wang | Discriminative indexing for probabilistic image patch priors[END_REF].

Apart from methods that accelerate EPLL, [START_REF] Houdard | High-Dimensional Mixture Models For Unsupervised Image Denoising (HDMI)[END_REF] and [START_REF] Guillemot | Covariance Trees for 2D and 3D Processing[END_REF] are two works that use strategies sharing similarities to our proposed flat-tail and binary search tree approximations, respectively. In [START_REF] Houdard | High-Dimensional Mixture Models For Unsupervised Image Denoising (HDMI)[END_REF], the authors use a flat-tail GMM to model distribution of noisy patches. Since the GMM is learned directly on noisy patches, the constant value of the tail corresponds to the noise variance. This allows to improve inference when learning on noisy datasets, estimating the noise level, and retrieving the intrinsic dimension of each cluster. This is different from our approach which uses GMM priors learned on clean patch and provides flat-tail approximation using mean of least significant coefficients. In [START_REF] Guillemot | Covariance Trees for 2D and 3D Processing[END_REF], the authors introduce a general data structure called covariance tree (CovTree). A CovTree is constructed by first building a binary space partition tree on patches (data points) and patches in each node are then modeled using a Gaussian distribution. While the resulting tree of Gaussians do share some similarities with our binary search tree, the learning process is very different. The covariance tree is built directly on data points, our approach is applied on an already learned GMM. In addition, unlike CovTree, our proposed strategy can handle any number of components (not just powers of two) and also encodes all parameters of mixture components including the mixing weights.

In the next section, we show that our proposed accelerations produce restoration results with comparable quality to competing methods while requiring a smaller amount of time.

V. NUMERICAL EXPERIMENTS

In this section, we present the results obtained on various image restoration tasks. Our experiments were conducted on standard images of size 512 × 512 such as Barbara, Boat, Couple, Fingerprint, Lena, Mandrill and on 60 test images of size 481×321 from the Berkeley Segmentation Dataset (BSDS) [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF] (the original BSDS test set contains 100 images, the other 40 was used for validation purposes while setting parameters ρ and s). For denoising, we compare the performance of our fast EPLL (FEPLL) to the original EPLL algorithm [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] and BM3D [8]. For the original EPLL, we have included timing results given by our own MATLAB/C implementation (EPLLc) and the MATLAB implementation provided by the authors (EPLLm). We also compare our restoration performance and runtime against other fast restoration methods introduced to achieve competitive trade-off between runtime efficiency and image quality. These methods include RoG [START_REF] Rosenbaum | The return of the gating network: combining generative models and discriminative training in natural image priors[END_REF] (a method accelerating EPLL based on feedforward networks described in Sec. IV), and CSF [START_REF] Schmidt | Shrinkage fields for effective image restoration[END_REF] (a fast restoration technique using random field-based architecture).

For deblurring experiments, we additionally compare with field-of-experts (FoE)-based non-blind deconvolution [START_REF] Chen | Insights into analysis operator learning: From patch-based sparse models to higher order mrfs[END_REF] denoted as iPiano. We contacted the corresponding author of [START_REF] Wang | Discriminative indexing for probabilistic image patch priors[END_REF] and got confirmation that the implementation of their algorithm (briefly described in IV) is not publicly available. Due to certain missing technical details, we were unable to reimplement it faithfully. However, the results reported in [START_REF] Wang | Discriminative indexing for probabilistic image patch priors[END_REF] indicate that their algorithm performs in par with BM3D in terms of both PSNR and time. Hence, BM3D results can be used as a faithful proxy for the expected performance of Wang et al.'s algorithm [START_REF] Wang | Discriminative indexing for probabilistic image patch priors[END_REF].

To explicitly illustrate the quality vs. runtime tradeoff of FEPLL, we include results obtained using a slightly slower version of FEPLL referred to as FEPLL , that does not use the balanced search tree and uses a flat tail spectrum approximation with ρ = 0.98. Please note that FEPLL is not meant to be better or worse than FEPLL, it is just another version running at a different PSNR/time tradeoff which allows us to compare our algorithm to others operating in different playing fields.

Finally, to illustrate the versatility of FEPLL, we also include results for other inverse problems such as devignetting, superresolution, and inpainting. Additional results can be found in the supplementary.

Parameter settings: In our experiments, we use patches of size P = 8 × 8, and the GMM provided by Zoran et al. [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] with K = 200 components. The 200-components GMM is progressively collapsed into smaller GMMs with K = 64, 32, 16, 8, 4, 2 and 1, and then all Gaussians of the tree are modified offline by flat-tail approximations with ρ = 0.95. The final estimate for the restored image is obtained after 5 iterations of our algorithm with β set to λσ -2 {1, 4, 8, 16, 32} where λ = min{N -1 ||A t A|| 2 F ||A|| -2 2 , 250σ 2 }. For denoising, where A is identity, λ = 1 which boils down to the setting used by Zoran et al. [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]. For inverse problems, we found that the initialization x = (A t A + 0.2σ 2 /λ∇) -1 A t y, with ∇ the image Laplacian, provides relevant solutions whatever the linear operator A and the noise level σ 2 . While the authors of [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] do not provide any further direction for setting β and the initialization in general inverse problems, our proposed setting leads to competitive solutions irrespective of A and σ 2 . For BM3D [8], EPLLm [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF], RoG [START_REF] Rosenbaum | The return of the gating network: combining generative models and discriminative training in natural image priors[END_REF], CSF [START_REF] Schmidt | Shrinkage fields for effective image restoration[END_REF] and iPiano [START_REF] Chen | Insights into analysis operator learning: From patch-based sparse models to higher order mrfs[END_REF] we use the implementations provided by the original authors and use the default parameters prescribed by them. EPLLc (our own MATLAB/C implementation), RoG [START_REF] Rosenbaum | The return of the gating network: combining generative models and discriminative training in natural image priors[END_REF] BM3D [8], CSF [START_REF] Schmidt | Shrinkage fields for effective image restoration[END_REF] and DnCNN [START_REF] Zhang | Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising[END_REF]. We evaluate the algorithms under low-, mid-and high-noise settings by using Gaussian noise of variance 5 2 , 20 2 and 60 2 , respectively. The result labeled "Berkeley" is an average over 60 images from the BSDS testing set [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF]. All numbers are averages obtained on 10 independent noise realizations. of 0.5dB may not affect the visual quality of the image. To illustrate this, we show a sample image denoised by BM3D, EPLL and FEPLL in Fig. 5.

Denoising:

Deblurring: Table III shows the performance of FE-PLL when used for deblurring as compared to RoG [START_REF] Rosenbaum | The return of the gating network: combining generative models and discriminative training in natural image priors[END_REF], iPiano [START_REF] Chen | Insights into analysis operator learning: From patch-based sparse models to higher order mrfs[END_REF] and CSF [START_REF] Schmidt | Shrinkage fields for effective image restoration[END_REF]. Once again, FEPLL uses Alg. 1 with operator A in the Image estimation step defined by

Ax = F -1 [F[h] F[x]],
where F is the fast Fourier transform and F -1 its inverse, indicates element-wise product and h is the blur kernel. For these experiments, we use the blur kernel provided by Chen et al. [START_REF] Chen | Insights into analysis operator learning: From patch-based sparse models to higher order mrfs[END_REF] along with their algorithm implementation. The results under the label "Berkeley" are averaged over 60 images from the BSDS test dataset [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF]. The results labeled "Classic" is averaged over the six standard images (Barbara, Boat, Couple, Fingerprint, Lena and Mandrill). FEPLL consistently outperforms its efficient competitors both in terms of quality and runtime. Although the GPU version of CSF is faster, the restoration quality obtained by CSF is 2-3dB lower than FEPLL. The proposed algorithm outperforms RoG by 1-1.8dB while running 3 and 5 times faster on "Berkeley" and "Classic" datasets, respectively.

The superior qualitative performance of FEPLL is demonstrated in Fig. 7. For brevity, we only include the deblurring results obtained from the top competitors of FEPLL algorithm in terms of both quality and runtime. As observed, FEPLL provides the best quality vs. runtime efficiency trade-off. In contrast, a deblurring procedure using the regular EPLL is around 350 times slower than FEPLL with the original implementation [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]. Specifically, on the sample image shown in Fig. 7, EPLL provides a qualitatively similar result (not shown in the figure) with a PSNR of 32.7 dB and SSIM of 0.922 in 142 seconds.

Other inverse problems: Unlike BM3D, EPLL and FEPLL are more versatile and handle a wide range of inverse problems without any change in formulation. In Fig. 8, we show the results obtained by FEPLL on problems such as (a) devignetting, which involves a progressive loss of intensity, (b) superresolution and (c) inpainting. To show the robustness of our method, the input images of size 481×321 were degraded with zero-mean Gaussian noise with σ = 2. All of the restoration results were obtained within/under 0.4 seconds and with the same set of parameters explained above (cf. Parameter settings).

VI. CONCLUSION

In this paper, we accelerate EPLL by a factor greater than 100 with negligible loss of image quality (less than 0.5dB). This is achieved by combining three independent strategies: a flat tail approximation, matching via a balanced search tree, and stochastic patch sampling. We show that the proposed accelerations are effective in denoising and deblurring problems, as well as in other inverse problems such as super-resolution and devignetting. An important distinction of the proposed accelerations is their genericity: the accelerated EPLL prior can be applied to many restoration tasks and various signalto-noise ratios, in contrast to existing accelerations based on learning techniques applied to specific conditions (such as image size, noise level, blur kernel, etc.) and that require an expensive re-training to address a different problem.

Since the speed-up is achieved solely by reducing the algorithmic complexity, we believe that further inclusion of accelerations based on parallelization and/or GPU implementations will allow for real-time video processing. Moreover, the acceleration techniques introduced in this work are general strategies that can be used to speed up other image restoration and/or related machine learning algorithms. For reproducibility purposes, the code of FEPLL is made available on GitHub 1 . artifacts start appearing especially when using a regular patch sub-sampling. In comparison, the stochastic sub-sampling leads to competing results even for s = √ P , which corresponds to a regular grid that is globally shifted by a random shift at each iteration of FEPLL. These experiments show that the choice of s = 6 leads to good trade-off in terms of speed and visual quality for a PSNR drop of less than 0.2dB.

Additional

results As the first set of additional results, we provide the visual results obtained on five of the standard images: Barbara, Boat, Couple, Fingerprint, Lena and Mandrill, in the setting σ = 20. Results are displayed in Supplemental Figure 6.

Deconvolution results

As the second set additional results, we show deconvolution results for five of testing images from the BSDS dataset [3] with a Gaussian blur setting of 3 pixels (standard deviation 3) and noise level σ = 2. Results are displayed in Supplemental Figure 7.

Super-resolution results

As the last set of additional results, super-resolution results are given for five of the testing images from the BSDS dataset [3] for super-resolution by a factor of 3 with a noise level of σ = 2. Note that the sub-sampling operator also includes a small Gaussian blur of 0.5 pixels as well as a Kaiser window apodization for anti-aliasing. See Supplemental Figure 8.

Fig. 1 .

 1 Fig. 1. Flat tail approximation: (a) with eigenvalues display on linear and (b) logarithmic scale. The r k = 15 first eigen components explains ρ = 95% of the variability.

Fig. 2 .

 2 Fig. 2. (left) Our obtained search tree (with numbers of nodes for each level). (right) Four patches well represented by each of the eight nodes along the branch highlighted in red in the tree. These patches are randomly sampled from the generative model encoded at each node of the branch (patches on the same column are generated from the same random seed).

Fig. 3 .

 3 Fig. 3. Illustration of patch subsampling. Instead of extracting all patches, only a subset of patches is extracted either (a) regularly or (b) with some randomizations. Patches are represented by 8×8 squares and the red intensity represents the number of patches overlapping the corresponding pixel.

Fig. 4 .

 4 Fig. 4. (top) Average speed-up and (bottom) average PSNR for our three accelerations, and all possible combinations of them, on the 40 images of the BSDS validation set.

Fig. 5 .

 5 Fig. 5. Illustration of a denoising problem with noise standard deviation σ = 20. Part of: (a) the original image (b) its noisy version (c-e) denoised results of competitive methods with PSNR/SSIM and time (inset).

Fig. 6 .

 6 Fig. 6. (a) PSNR and (b) SSIM versus time for different restoration methods in a denoising problem with noise standard deviation σ = 20. PSNR, SSIM and time are averaged on the 60 BSDS test images, each of size 481 × 321. Optimal methods tend to be in the top-left corner.

Fig. 7 .

 7 Fig. 7. Illustration of a deblurring problem with noise standard deviation σ = 0.5. Part of: (a) the original image and the blur kernel (inset), (b) blurry version (c-e) deblurred results of competitive methods with PSNR/SSIM and time (inset). The '*' indicates runtime on GPU while others are CPU times.

  with the other accelerations (c) Reference (d) Stochastic s = 2 (e) Stochastic s = 4 (f) Stochastic s = 6 (g) Stochastic s = 8 (h) Noisy (i) Regular s = 2 Regular s = 4 (k) Regular s = 6 (l) Regular s = 8 Supplemental Figure PSNR as a function of speed-up for varying values of (top) the patch extraction expected period s (a) no other accelerations are used and (b) the other two accelerations are also used. Stochastic and regular sub-sampling are compared. The parameter offering the best speed-up within a drop of at most -0.2dB is indicated. (c) A reference image, (h) its noisy version with σ = 20 and results obtained for increasing value of s (including the other two accelerations as well) for (d-g) the stochastic version and (i-l) the regular version.

TABLE I COMPARISON

 I OF THE EXECUTION TIME OF OUR IMPLEMENTATION OF EPLL WITH AND WITHOUT PROPOSED ACCELERATIONS. EXPERIMENT CONDUCTED ON 40 IMAGES EACH OF SIZE 481 × 321 IN DENOISING SETTING. PROFILING WAS CARRIED OUT USING MATLAB (R2014B) ON A PC WITH INTEL(R) CORE(TM) I7-4790K CPU @4.00GHZ AND 16 GB RAM. EXECUTION TIMES ARE REPORTED AS AVERAGE NUMBER OF SECONDS PER IMAGE (S) AND PERCENTAGE OF THE TOTAL TIME (%).

	Step	Without accelerations		With the proposed accelerations
	(Gaussian selection) 38.27s (Patch estimation) 0.95s (Patch extraction) 0.43s (Patch reprojection) 0.15s Others 0.52s Total 40.32s	2 % 1 % 0 % 1 %	95 %	0.21s 0.04s 0.02s 0.01s 0.03s 0.31s	13 % 6 % 4 % 10 %	67 %

TABLE II PSNR

 II , SSIM AND EXECUTION TIME ON THE BSDS TEST SET (AVERAGE ON 60 IMAGES OF SIZE 481×321), AND ON SIX STANDARD IMAGES (EACH OF SIZE 512×512 AND RESULTS AVERAGED OVER 10 INDEPENDENT NOISE REALIZATIONS) FOR THE PROPOSED FEPLL AND FEPLL , EPLL (WITH TIMING GIVEN FOR BOTH EPLLM[START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] AND OUR EPLLC), BM3D[8], CSF[START_REF] Schmidt | Shrinkage fields for effective image restoration[END_REF], ROG[START_REF] Rosenbaum | The return of the gating network: combining generative models and discriminative training in natural image priors[END_REF] AND DNCNN[START_REF] Zhang | Beyond a Gaussian Denoiser: Residual Learning of Deep CNN for Image Denoising[END_REF] WITH 3 DIFFERENT LEVELS OF NOISE.

	σ	Algo.	Berkeley	Barbara	Boat	Couple	Fingerprint	Lena	Mandrill
					PSNR/SSIM
	5 20 60	FEPLL FEPLL EPLL RoG BM3D CSF 3×3 DnCNN FEPLL FEPLL EPLL RoG BM3D CSF 3×3 DnCNN FEPLL FEPLL EPLL RoG BM3D CSF 3×3 DnCNN	36.8 /.959 37.1 /.962 37.3 /.963 37.1 /.959 37.3 /.962 36.8 /.952 34.5 /.915 29.1 /.812 29.3 /.831 29.5 /.836 29.4 /.828 29.4 /.824 29.0 /.805 30.0 /.847 24.5 /.614 24.5 /.620 24.8 /.631 24.6 /.623 24.8 /.637 22.0 /.489 25.3 /.665	36.9 /.958 36.6 /.930 36.6 /.944 37.1 /.959 36.7 /.933 36.8 /.946 37.6 /.962 36.8 /.933 37.3 /.950 36.7 /.954 36.5 /.920 37.2 /.946 38.3 /.964 37.3 /.939 37.4 /.949 37.0 /.955 36.7 /.929 37.1 /.945 35.5 /.937 34.4 /.878 34.5 /.893 29.0 /.852 30.2 /.802 29.9 /.813 29.5 /.866 30.3 /.814 30.1 /.825 29.8 /.872 30.6 /.821 30.4 /.834 28.4 /.838 30.5 /.815 30.3 /.827 31.7 /.904 30.8 /.824 30.7 /.842 28.3 /.821 30.2 /.802 29.9 /.812 31.0 /.896 31.1 /.831 31.0 /.848 23.6 /.636 25.5 /.644 25.1 /.629 23.8 /.648 25.5 /.646 25.1 /.636 24.0 /.660 25.8 /.659 25.4 /.649 23.3 /.626 25.6 /.654 25.2 /.640 26.3 /.757 25.9 /.671 25.6 /.666 21.4 /.490 22.7 /.502 22.5 /.505 24.9 /.717 26.4 /.686 26.0 /.683	35.4 /.984 35.5 /.984 36.4 /.987 36.4 /.987 36.5 /.987 36.2 /.986 32.9 /.968 27.6 /.908 27.8 /.916 28.3 /.924 28.3 /.922 28.8 /.928 28.0 /.916 28.8 /.930 22.2 /.722 22.4 /.745 22.6 /.755 22.4 /.739 23.8 /.801 21.2 /.727 23.1 /.781	38.2 /.941 35.1 /.959 38.3 /.943 35.2 /.960 38.6 /.944 35.2 /.960 38.3 /.939 35.1 /.956 38.7 /.944 35.3 /.959 38.2 /.938 34.8 /.953 36.3 /.907 32.3 /.899 32.3 /.863 26.4 /.784 32.4 /.864 26.5 /.802 32.6 /.869 26.7 /.807 32.5 /.865 26.5 /.794 33.0 /.876 26.6 /.794 31.8 /.837 26.1 /.778 33.4 /.883 27.0 /.815 27.4 /.742 21.4 /.468 27.3 /.733 21.6 /.499 27.6 /.747 21.7 /.506 27.4 /.747 21.5 /.482 28.2 /.778 21.7 /.500 23.4 /.511 20.2 /.466 28.5 /.790 22.1 /.555
					Time (in seconds)
		FEPLL FEPLL EPLLc EPLLm RoG BM3D CSF 3×3 CSF gpu 3×3 DnCNN DnCNN gpu	0.27 0.96 43.82 82.68 1.16 1.60 0.87 0.38 2.09 0.28	0.38 1.28 71.14 145.15 1.92 2.52 1.09 0.39 3.57 0.63	0.38 1.28 71.24 143.68 1.92 2.68 1.09 0.39 3.55 0.63	0.37 1.28 71.31 144.30 1.93 2.59 1.13 0.39 3.60 0.64	0.36 1.26 71.28 144.13 1.91 2.17 1.09 0.39 3.60 0.63	0.38 1.27 71.23 144.15 1.93 2.61 1.09 0.39 3.60 0.64	0.38 1.28 71.31 143.86 1.92 2.64 1.08 0.39 3.61 0.64
		22.1 / .368	30.6 / .872 (1.68s)			30.6 / .873 (44.88s)	30.2 / .862 (0.36s)
	(a) Reference x	(b) Noisy image y	(c) BM3D result	x	(d) EPLLc result	x	(e) FEPLL result	x

  Table II shows the quantitative performances of FEPLL on the denoising task compared to EPLLm [43],

	30
	29.8
	29.6
	29.4
	29.2
	29
	28.8
	10 -1	10 0	10 1	10 2
		(a)
	0.85
	0.84
	0.83
	0.82
	0.81
	0.8
	0.79
	10 -1	10 0	10 1	10 2

TABLE III PSNR

 III , SSIM AND EXECUTION TIME ON THE BSDS TEST SET (AVERAGE OF 60 IMAGES OF SIZE 481×321), AND ON STANDARD IMAGES (AVERAGE OF 6 IMAGES OF SIZE 512×512) FOR THE PROPOSED FEPLL AND FEPLL , CSF[START_REF] Schmidt | Shrinkage fields for effective image restoration[END_REF], ROG[START_REF] Rosenbaum | The return of the gating network: combining generative models and discriminative training in natural image priors[END_REF] AND CHEN ET AL.'S[START_REF] Chen | Insights into analysis operator learning: From patch-based sparse models to higher order mrfs[END_REF] METHOD THAT IS CALLED iPiano IN THEIR IMPLEMENTATION. THE BLUR KERNEL USED IS THE ONE PROVIDED IN ALONG WITH IPIANO IMPLEMENTATION AND NOISE IS SET TO σ = 0.5.

	Algo.	Berkeley	Classic
		PSNR/SSIM	Time (s)	PSNR/SSIM	Time (s)
	iPiano CSFpw RoG FEPLL FEPLL'	29.5 / .824 30.2 / .875 31.3 / .897 33.1 / .928 33.2 / .930	29.53 0.50 (0.14*) 30.5 / 0.870 29.9 / .848 1.19 31.8 / .915 0.40 32.8 / .931 1.01 33.0 / .933	59.10 0.47 (0.14*) 2.07 0.46 1.82

Supplementary for "Accelerating GMM-based patch priors for image restoration:

Three ingredients for a 100× speed-up"

Shibin Parameswaran * , Charles-Alban Deledalle * † , Loïc Denis ‡ and Truong Q. Nguyen * .

August 8, 2018

This document gives additional results and analysis of FEPLL algorithm described on "Accelerating GMM-based patch priors for image restoration: Three ingredients for a 100× speed-up".

Additional analysis

Flat tail approximation based acceleration First, we demonstrate the superiority of the proposed flat tail approximation over the commonly used approach of ignoring the least significant eigendirections. In Supplemental Figure 1, we show that the naïve approach of zeroing out (or ignoring) coefficients of the least significant eigendirections is inferior to the proposed approach of replacing these coefficients by the mean. Our approach can provide reasonably good performance even for very low values of ρ. For example, when ρ = 0.20 our approach leads to a PSNR of 27.5 compared to 25.1 obtained by the zeroing out strategy.

In a second experiment, we analyze the influence of the parameter ρ ∈ (0, 1] used in the flat tail approximation. Supplemental Figure 2 shows the curve of PSNR as a function of speed-up for varying values of ρ. This experiment is repeated twice either with or without enabling the other two accelerations. Visual results highlight that as ρ decreases to zero, residual noise starts appearing around salient structures. From these experiments, the choice of ρ = 0.95 leads to a good trade-off in terms of speed and visual quality for a drop of PSNR lower than 0.2dB.

Search tree acceleration

In the second supplementary experiment, we analyze the influence of the clustering method used to successively collapse the GMM model at a given level into a smaller model at a lower level. We compare the proposed approach based on Multiple Traveling Salesmen Problem [2] (MTSP) using the symmetric KL divergence as a semi-metric to a hierarchical GMM clustering algorithm similar to K-means based on KL divergence as proposed in [1] and the standard hierarchical agglomerated clustering (HAC) using symmetric KL divergence. Supplemental Figure 3 shows that MTSP leads to a well-balanced tree, with a height of 7 (almost a binary tree except for the last level due to 200 mixture components). In comparison, using [1]'s K-Means inspired algorithm provides a tree of height 7 but is not balanced and HAC leads to a tree of height 59 with comb structured branches. Please note that a shorter tree is preferred for a faster Gaussian selection step. The tree built using the proposed MTSP strategy not only leads to better PSNR/SSIM, but also provides a more stable computation time for all images of a fixed size (irrespective of content) due to its balanced structure. This is confirmed in Supplemental Figure 4 Supplemental Figure 1: Comparison between the naïve approach of zeroing out (ignoring) coefficients of the least significant eigendirections versus the proposed flat-tail approximation approach of replacing the coefficients with mean. Performance is compared in terms of PSNR obtained while denoising with different values of ρ (percent energy captured) used for thresholding singular values. obtained from computation time statistics based on five independent runs of the algorithm on 40 different images of the BSDS dataset (all images have the same size). In contrast to our MTSP based strategy, the trees built using [1]'s method and HAC strategy lead to computation times that vary drastically depending on the image content.

Stochastic patch sub-sampling The last analysis focuses on the influence of the period 1 s √ P used in the stochastic patch sub-sampling. In addition, we also compare against regular sub-sampling. Supplemental Figure 5 shows the curve of PSNR as a function of speed-up for varying values of s where P = 8. The experiment has been performed in both cases where either the other two accelerations are disabled or enabled. Visual results highlight that as s increases to √ P , blocky Supplemental Figure 4: Timing scatter obtained over 5 runs for 40 images using Gaussian trees built using different clustering methods we tested. "MTSP" refers to the Multiple Traveling Salesman (proposed), "K-Means" is the method introduced by [1] and "HAC" is a simple Hierarchical agglomerative clustering approach.


Luis Sanchez
email: luis.sanchezgu@usach.cl

Erika Rosas
email: erosas@inf.utfsm.cl

Nicolas Hidalgo
email: nicolas.hidalgoc@mail.udp.cl

Crowdsourcing under Attack: Detecting Malicious Behaviors in Waze

Keywords: Online Social Networks, Malicious Behaviors, Sybil attack, Sybil detection

Social networks that use geolocalization enable receiving data from users in order to provide information based on their collective experience. Specifically, this article is interested in the social network Waze, a real-time navigation application for drivers. This application uses methods for identifying users that are open and free, where people are able to hide their identity by using a pseudonym. In this context, malicious behaviors can emerge, endangering the quality of the reports on which the application is based. We propose a method to detect malicious behavior on Waze, which crawls information from the application, aggregates it and models the data relationships in graphs. Using this model the data is analyzed according to the size of the graph: for large interaction graphs, we use a Sybil detection technique, while for small graphs we propose the use of a threshold-based mechanism to detect targeted behaviors. The results show that it is complex to use the large-scale Sybil attack detection techniques due to parameter tuning. However, good success rates can be achieved to tag users as honest and malicious if there are a small number of interactions between these groups of users. On the other hand, for small graphs, a straightforward analysis can be performed, since the graphs are sparse and the users have a limited number of connections between them, making clear the presence of outliers.

Introduction

In recent years, Online Social Networks (OSN) usage has increased exponentially and some of them have become part of most people's everyday life. Facebook 4 , Instagram 5 , Twitter6 and Waze7 are only a few examples of OSNs accounting for millions of users that interact daily creating and sharing information. In large OSNs, threats are just around the corner, not only menacing users' private data, but also the whole network goals. Identity theft, malwares, fake profiles (or Sybils) are common examples of threats present in this type of networks [START_REF] Fire | Online social networks: Threats and solutions[END_REF]. In this work, we are particularly interested in Sybil-based attacks under the well-known OSN for drivers, Waze.

Waze is a crowdsourcing application to assist drivers by providing online information on traffic and road conditions. It was created in 2008 and to-date, it has approximately 50 million users [START_REF] Goel | Maps that live and breathe with data[END_REF]. Waze creates an online report of traffic conditions for a given route based on the information collected from and reported by users: current speed, position, origin and destination, police controls, traffic jams, accidents, etc. One of Waze's main features is user engagement to contribute to the common good, i.e., Waze is not just crowdsourcing, but personal participation [START_REF] Goel | Maps that live and breathe with data[END_REF]. Waze's success is directly related to the good-will of users; therefore, malicious behaviors such as Sybil attacks can seriously compromise the application's precision and success.

In the last years, some of these attacks have been reported in [START_REF] Carney | Angry LA residents are trying to sabotage Waze data to stop side-street congestion[END_REF][START_REF] Sinai | Exploiting Social Navigation[END_REF]. In the former, researchers were able to generate fake users and mobile devices using virtual machines in Android and they created fake traffic jams by setting low speed configurations to their fake devices. In the latter, people coordinated to emulate a traffic jam in a residential area, in order to reduce car passing in their neighborhood. Attacks like these can seriously compromise the behavior of the application and detecting this can be challenging in the presence of millions of users dynamically interacting online with the presence of anonymous users in the system. Sybil detection in similar environments such as Twitter have been studied; however, in the context of Waze, the application of state-of-the-art Sybil detection mechanisms requires modelling the malicious behaviors in terms of the data collected from the network.

In this work, we propose three models that attempt to characterize three different behaviors on which we focus our study: (1) Collusion for traffic jams: people collude in order to simulate a traffic jam by not moving [START_REF] Carney | Angry LA residents are trying to sabotage Waze data to stop side-street congestion[END_REF]; (2) Driver speed attack: a coordinated group of Sybils simulate slow driving so that Waze declares a false traffic jam [START_REF] Sinai | Exploiting Social Navigation[END_REF]; and (3) False event attacks: a coordinated group of Sybils vote for a false event, that obscures honest users.

The main contributions of this work are the models associated to the three malicious Sybil behaviors previously exposed. Our models were tested with real Waze traces and our results that show that malicious behaviors can be detected using a state-of-the-art Sybil attack detection mechanism and a threshold-based mechanism. In our experiments, we have exploited SybilDefender [START_REF] Wei | SybilDefender: Defend against sybil attacks in large social networks[END_REF] and a threshold-based mechanism to detect abnormal behaviors. The former is applied over large interaction graphs and the latter, over small ones.

The remainder of this work is organized as follows: Section 2 briefly introduces the identity problems and their relationship with Sybil attacks. Also, we present some literature works that attempt to tackle this problem. Section 3 presents the main contribution of this work, the graphs that model the malicious behaviors we attempt to detect. In Section 4, we evaluate the proposed malicious behavior models using SybilDefender as the mechanism to detect Sybils and a thresholdbased mechanism for those small interactions graphs. Finally, our conclusions and future work is stated in Section 5.

Background and Related Work

Identity Attacks in Social Networks

An identity in a social network is the set of characteristics of a particular person or group (entity) that distinguishes it from others in the network. Different to real-life identities, such as identification cards or a passport that are shown by a person and that can be confirmed comparing a picture and the biometric indexes, in the online world it is more difficult to establish the link between a physical entity and the online identity that represents it.

This problem has been widely discussed because it is easy to change one's identity in several OSNs, whereas in real life, this is a complex process. Friedman and Resnick have called this type of identity cheap pseudonyms [START_REF] Friedman | The Social Cost of Cheap Pseudonyms[END_REF], and allows a person interacting anonymously to constantly change identifiers or to maintaining a persistent identity.

In this context, one unique entity can build a set of pseudonyms in the system, which makes it appear as different entities. What we call a Sybil attack occurs when one physical entity creates and uses a set of identities in the system in order to perform malicious behaviors [START_REF] Douceur | The sybil attack[END_REF]. The malicious behaviors may vary according to the online environment and it may range from exploiting more resources than allowed to performing active attacks that hamper the veracity of the information exchanged in the network.

The problem of cheap pseudonyms is that they reduce the number of accountable actions in the system, and in the case of Waze, one user that has multiple identities in the system or multiple users may collude to spread false information for other drivers. This network shows relevant information to the users based on their localization and false information may show false traffic jams, which may produce longer routes for other drivers.

This kind of problem has been studied by [START_REF] Zhang | You are where you have been: Sybil detection via geo-location analysis in OSNs[END_REF] in the context of the social network Dianping. They have found that some user accounts make positive comments about places that are very far from each other in time intervals that are impossible to achieve. A few users control these accounts that give good rating to some places and bad rating to the competence, in exchange for money.

In the case of Waze, Sinai et al. [START_REF] Sinai | Exploiting Social Navigation[END_REF] coordinated a Sybil attack by creating multiple identities using multiple virtual machines in Android that ran the Waze application. They simulated slow driving in all the identities on a specific street so the system detected a false traffic jam. They have proved that it is possible to control routes, which may produce important problems for other drivers. A collusion attack has also been documented by Carney [START_REF] Carney | Angry LA residents are trying to sabotage Waze data to stop side-street congestion[END_REF] in LA, USA. In this case, the neighbors colluded and activated Waze outside their houses in order to simulate a false traffic jam, in order to force Waze to not recommend that neighborhood to drivers.

Managing Sybil Attacks and Collusions

In the context of large-scale systems, we can find two types of approaches to prevent the Sybil attack and collusion: detection and tolerance. The detection of malicious behaviors is focused on detecting identities that are acting with malice in the system, and evict them from the system. However, in presence of cheap pseudonyms, there is no problem for them to obtain a fresh new identity. For this reason, in large-scale networks, a more common approach is to tolerate malicious behaviors, for example, by avoiding using information generated by suspicious identities.

Community detection techniques have been used to detect suspicious entities in the presence of the Sybil attacks, assuming that the number of interactions of these types of identities with real users is limited, that there is at least one honest known user in the network, and the honest region is densely connected. Several Sybil detection algorithms have been proposed [START_REF] Danezis | SybilInfer: Detecting Sybil Nodes using Social Networks[END_REF][START_REF] Wei | SybilDefender: Defend against sybil attacks in large social networks[END_REF][START_REF] Yu | Sybilguard: Defending against sybil attacks via social networks[END_REF] that classify identities as Sybils or normal.

In this work, we used SybilDefender [START_REF] Wei | SybilDefender: Defend against sybil attacks in large social networks[END_REF], but any other Sybil detection system can be used in its place. SybilDefender proposes four algorithms; the first one obtains statistics from the neighborhood of an honest node identifying J judges from its vicinity and performing R random walks of length L. The second one identifies a suspicious node as Sybil or non-Sybil using the results of the first algorithm. This is performed through R random walks of length L from the suspicious node, and comparing the values of recurrent nodes in the results of algorithm 1. The third and fourth algorithms enable detecting a Sybil region around a node classified as Sybil. A detailed description of the algorithms is presented in [START_REF] Wei | SybilDefender: Defend against sybil attacks in large social networks[END_REF].

Recently, Sybil attacks have been studied in the context of Vehicular Ad-hoc Networks (VANETs) [START_REF] Feng | A method for defensing against multisource sybil attacks in vanet[END_REF][START_REF] Gu | Vehicle driving pattern based sybil attack detection[END_REF][START_REF] Gu | Support vector machine (svm) based sybil attack detection in vehicular networks[END_REF]. In [START_REF] Feng | A method for defensing against multisource sybil attacks in vanet[END_REF], the authors built an event-based reputation system that feeds a trust management system that restricts the dissemination of false messages. In [START_REF] Gu | Vehicle driving pattern based sybil attack detection[END_REF][START_REF] Gu | Support vector machine (svm) based sybil attack detection in vehicular networks[END_REF], the authors use driving patterns of vehicles and detect Sybils using classifiers, such as minimum distance classifier and support vector machines. In location-based social networks, [START_REF] Xu | Towards efficient detection of sybil attacks in location-based social networks[END_REF] states that it is not normal the appearance of continuous gatherings, and use the detection of these events to identify Sybils. However, we argue that traffic jams may produce continuous gathering in urban zones. Other graph-based solution has been proposed for mobile online social networks in [START_REF] Chinchore | Intelligent sybil attack detection on abnormal connectivity behavior in mobile social networks[END_REF], where the authors use a connection analysis to differentiate honest versus fake nodes. Finally, in the context of mobile crowdsourcing, recent work proposes a passive and active checking scheme that verify traffic volume, signal strength and network topology [START_REF] Chang | Protecting mobile crowd sensing against sybil attacks using cloud based trust management system[END_REF], differentiating nodes using an adaptive threshold.

Modeling Sybil behaviors

The goal of the model is to identify Sybil attacks in Waze. The key contribution is the way of modeling the data of Waze in order to detect Sybils or collusion. Figure 1 shows the proposed pipeline that is detailed in the following subsections. In general, data is captured from the LiveMap of Waze, reordered, and aggregated in order prevent redundant information. Then, we generate graphs that model the interactions between the users according to the malicious behavior we target. Then, an analysis of the graphs is performed using state-of-the-art Sybil detection mechanisms.

Crawling and indexing

We crawled data from Waze using the endpoint that feeds the LiveMap 8 of Waze. This is public data that is delivered in JSON format of the current state of the requested area, defined by coordinates. The data captured was requested every 1 minute, obtaining a snapshot of the map at that time. We categorize data in three types:

-Alerts: Alerts are events explicitly reported by Waze users, such as vehicle accidents, police locations, traffic jam reports, etc. These alerts are characterized by a point on a map and the number of votes received by other users that corroborate the information. -Jams: Traffic jams are events that Waze reports using the location data of users. They are created when large traffic is detected on a street and are represented by geographic coordinates that build a line, and other data such as the severity of the traffic jam and current user speed, among others. -Users: The data of users on the map shows their current location and is represented by an identifier, geographic coordinates, and current speed, among others.

We have indexed the data by user identifier and event identifier. The information about traffic jams is not used in this study, since it does not contain user information, which is the main focus of this study.

Data Aggregation

We mainly produce two structures that facilitate graph modeling: In order to build the relationship between users that interact with a real event, we fuse the votes that are close in time and space.

We have set the time and space parameters for 30 minutes and 200 meters respectively. -Routes generation: In order to relate users to a route trip they have followed, we have generated routes from the geographic coordinates of the users. The routes are built from the user geographic points that are relatively close in time. Since the data is captured every 1 minute, it is not easy to generate a route from two geographically different points. We have used the process called map matching first, in order to deal with vagueness of the GPS location, so that the points are inside a street. Then, we have used Google Roads API, to correct the coordinates and infer the route that the user took.

Graph modelling

As we attempt to exploit state-of-the-art Sybil detection mechanisms, we must model the data associated to the different malicious behaviors that we attempt to detect. In this work, we focus on three problems, some of them already reported in literature. The malicious behaviors we target are:

-Collusion for traffic jam: A group of users is detained on a street so that Waze declares a false traffic jam. This behavior is described in [START_REF] Carney | Angry LA residents are trying to sabotage Waze data to stop side-street congestion[END_REF]. -Driving speed attack: A coordinated group of Sybils that simulate slow driving so that Waze declares a false traffic jam. This behavior is described in [START_REF] Sinai | Exploiting Social Navigation[END_REF]. -False event attacks: A coordinated group of Sybils that vote for a false event that obscure the honest users.

The generated interaction graphs are non-directed defined as G = (V, A), where V is the set of vertexes that represent users, and A indicates the set of edges that represent interactions between two users. The weight of the edge is related with the amount of interactions the users had in time.

If the modeled graphs are fast mixing and the number of connections from the malicious to the honest area is limited, then, we can apply a Sybil Detection algorithm. With these properties, we guarantee that random walks can iterate all over the honest area, and it is difficult to walk outside this area.

Traffic jam graph This graph is built with the users that do not move at the same time and within a close distance. Then, an interaction between two users a, b is defined as the number of times they were standing still at the same time in any of their trips v. Equation (1) shows how the weight of an edge between two users is computed: the sum of the number of times they may be colluded. ST is the set of trips where the users had the same origin and destination, meaning that they did not move from the beginning of the trip.

Equation ( 2) shows when a possible collusion is detected: when in a time window t the users were within a distance of less than d meters. We consider that the function distance(v 1 , v 2 ) gives the distance in meters from the position of the trip v 1 to the position of the trip v 2 and time(v) gives the time when the travel v happened. If the users collude to stand still in the streets in order to produce a traffic jam, then they will appear more connected in the graph, and we would like to detect this malicious area.

wg1(a, b) = vua ,v ub ∈ST col(vua, v ub ) (1) col(v1, v2) =      1, if distance(v1,
Driving speed graph The goal of this graph is to detect inconsistencies in the user behaviors that share part of their travel routes. We compare speed and temporality in order to check if their characteristics validate each other. The weight of an edge in this graph is high when two users share routes close in time and their driving speeds were similar. A lower weight of an edge indicates that their driving speeds were very different. Equation (3) shows the weight of the edge between users a and b: In the sum of all the times there is a similarity found between their trips T , normalized by the minimum number of trips of the users. We call v ua a trip of the user a and v ub a trip of the user b. Equation ( 4) shows what we consider a similarity in this case. We take into account three properties of the trips, the time, the routes and the speed. Route similarity s routes (v 1 , v 2 ) computes the number of segments shared in the trip v 1 and v 2 , divided by the number of segments of the longest trip between v 1 and v 2 . This route similarity is modified by factor α and β so to favor a homogeneous distribution, incrementing the similarity to 1 in cases there is a high similarity in speed and time, and generating a medium effect when there is a regular similarity in speed and time.

The temporal similarity s time (v 1 , v 2 ) is equal to zero if v 1 and v 2 were separated in time more than a value min and equal to 100% if v 1 and v 2 occurred at the same time. Temporal distances in between are proportionally computed. The speed similarity s speed (v 1 , v 2 ) indicates whether the segments that share v 1 and v 2 have a similar speed. This is computed with the average speed of these segments.

wg2(a, b) = vua,v ub ∈T sim(vua, v ub ) min(∪ i=a,b number of user trips i) (3) sim(v1, v2) =                sroutes(v1, v2) × α, if s speed (v1,
α = 1 max(∪ n i=0 wg2) (5)
β = 1 + ( 1 -α 2 ) (6)
In extreme cases of malicious behaviors, the users that share temporality and space will be strongly connected. Honest users would have low weighted edges.

False event graph

The main goal of this graph is to identify user groups that vote on the same events. In this case, the weight of each edge w g3 (a, b) that links the users a and b is the sum of the events that the users a and b have voted for the same event.

Malicious users that also vote for events that honest users have voted on may create a graph that is not fast mixing and hinder Sybil detection. However, they require the effort of creating strong links with other users. Figure 2 shows a random subgraph built with the experimental data. Each vertex represents a user and the links represent the interaction the users had in time. The weight of these edges is determined by the number of votes of the users on the same events.

Malicious Behavior Detection

The malicious behavior was detected by analyzing the graphs and applying a threshold based approach or a Sybil detection algorithm, according to the characteristics of the graph. We will tag a user as malicious when: Honest users that do not collude with others should present a small number of connections. -In the second graph, strong relations are given when users share their behavior, which can be honest or malicious. We have to start the process with a previously known honest user to identify the regions. -The third graph is similar to the second; we assume that users do not vote in group for an event, so the number of connections between them are going to be smaller compared to malicious Sybil users.

Experimentation

The data was crawled in a timespan of 6 months, from July 2015 until January 2016, with some data missing in October 2015 due to a server failure. We used the LiveMap API of Waze that was consulted every 1 minute. The area consulted was around the city Santiago, Chile with the coordinates -33.2 North, -33.8 South, -70.87 East, -70.5 West.

In total, we considered 1, 667, 400 events, that where generated by 223, 031 users. We captured 4, 547, 887 users with locations inside the coordinates, which is a large number, considering the size of the city of Santiago; however, Waze uses new identifiers for anonymous users, which explains the number of users.

In the data aggregation step, 30% of the events were fused. Figure 3 shows the distribution of the events per hour of the day. We are able to observe the peak time of the day in the figure, which is normal for a city like Santiago, at times when people go to work and when they return home.

In the data aggregation step, we have also obtained user trips with the trace we built from their locations that are close in time (maximum 7 minutes of difference between each subsequent pair). We consider trips that have at least 5 consecutive location points. These values were experimentally chosen, since most of the trips have an average of 1 minute between consecutive points. The result is 192, 248 trips of 184, 992 users (mostly anonymous users). Figure 4 and 5 show the distribution of the duration and the distance traveled in each trip. Most trips are short in time and distance.

False Event Graph

The obtained graph has 223, 030 vertexes and 1, 452, 261 edges, with an average degree of 13. This is a non-connected graph, and has 59, 815 components. In order to perform the experiments we selected the largest component of the graph that has 160, 894 vertexes and 1, 499, 717 edges. The remaining graphs have a size that is negligible for this experiment.

In order to obtain a fast mixing graph, we removed vertexes with small degrees. Figure 6 shows how the mixing time changes with the minimum degree of the graph. This means that we can apply a Sybil detection mechanism when the graph is dense, with a minimum degree of 64. The resulting graph has 9, 743 vertexes and 720, 152 edges with an average degree of 73.

We connected a Sybil area to this graph in order to observe how the algorithm behaves. The Sybil area was created using the Erdõs-Rényi model to create random and sparse graphs. This model is often invoked to capture the structure of social networks [START_REF] Dobrescu | Large Scale Networks: Modeling and Simulation[END_REF] and is defined as a set of N nodes connected by n edges chosen random from the (N (N -1))/2 posible edges.The Sybil area had 523 users and 1, 080 edges, with a mixing time of 41. The number of vertexes corresponds to 0.03% of the number of honest nodes in the event graph. To connect both regions we created random attack edges between each region.

SybilDefender works mainly with two algorithms: the first identifies J judges from an honest known node. Then, from each judge, R random walks are computed of length L, counting how many nodes appear more than T times in the random walks. Then, the average and the standard deviation are computed for each element of the list L. In our case, we set the number of judges to T = 50, the number of random walks to R = 100, L = 100, 200 . . . 1000 and T = 5, considering the size of the experiments presented in [START_REF] Wei | SybilDefender: Defend against sybil attacks in large social networks[END_REF].

The second algorithm identifies if a node is Sybil or not, using the results of the first algorithm. This is done computing R random walks of length L from the suspicious node, counting how many nodes exceed T repetitions, which is called m. Then, a comparison is done to determine if the node is Sybil or not: media -m > deviation × α.

SybilDefender, with 20 attack edges joining the Sybil and honest areas, found a 99.1% of successfully detected Sybils and 100% successfully detected honest users using SybilDefender algorithms. However, if the number of attack edges We modified the parameters of the algorithms, in order to observe the influence on the results. When using a T = 7 and changing the value α we obtained the results of Figure 7. The best results are obtained when using α = 55, when we obtained a 76% of Sybils detected and 96% of honest users.

Furthermore, if we modify the length of the random walks, we obtain the results of Figure 8. In this case, the algorithm always detects all the honest users and the best case for Sybils occurs at 350 maximum length where it finds 83% of Sybils.

Driving Speed Graph

This graph was built with parameters min and max set on 30 and 70, respectively.

The resulting graph has 149, 492 vertex and 7, 048 edges. This is a non-connected graph with 142, 582 components. Unlike the previous graph, this one has 138, 810 components with one vertex, 2, 479 with two vertexes and 1, 293 with three or more vertex, with the largest component with 24 vertexes. The size of the graph is too small to apply Sybil defender. Figure 9 shows the value of the edges in the largest component. Analyzing the edges values, we located an edge with a similarity value of 44,51. This is because they had a very different speed in the shared segments of the trip.

We have drawn all the travels involved in the generation of the graph, and obtained Figure 10. It is clear that the graph show honest users that have similarities in the highways of the city.

Traffic Jam Graph

In this case, the parameter d was set to 1km and t equals to 30 minutes. With this parameters we obtained a graph for this case has 38, 455 vertex and 8, 592 edges, with an average degree of 0.44 (maximum degree 14). This is a non-connected disperse graph that has 31, 967 connected components. The largest part has 711 vertex and 1, 448 edges, which is small for applying Sybil defender. Figure 11 shows the locations of the users of the graph, mostly located in the northeastern area of Santiago. We analyzed the distribution of the edges weight and found that from 711 users, there are 643 that have only one stand still trip, 58 have only two stand still trips, 9 have between 3 and 8 trips and only one user have 27 stand still trips. Thus, we conclude that colluded users may be detected by their degree in the graph, which in our case can be set around 10 to consider a user as suspicious.

Conclusion

In this work, we have proposed a model to process Waze traces and detect Sybil behaviors. The model consists of five steps, standing out the modeling of targeted behaviors on an interaction graph. We provide three of these models that characterize the three different behaviors which we focus on in our study: (1) Collusion for traffic jam, where people collude in standing still to simulate a traffic jam and divert traffic; (2) Driver speed attack, where a coordinated group of Sybils simulate slow driving to trick the application into assuming a traffic The general model was tested with real Waze traces and our results that show that malicious behaviors can be detected using a state-of-the-art Sybil attack detection mechanism and a threshold-based mechanism. In our experiment, we have exploited SybilDefender [START_REF] Wei | SybilDefender: Defend against sybil attacks in large social networks[END_REF] and a threshold-based mechanism to detect abnormal behaviors. The former is applied on large interaction graphs and the latter, over small ones where the application of a large-scale detection mechanism is not necessary.

Our results show that it is complex to use the large-scale Sybil attack detection techniques due to the parameter tuning. However, good success rates can be achieved to tag users as honest and malicious, if the number of interactions between those groups of users is small. On the other hand, for small graphs, a straightforward analysis can be performed since the graphs are sparse and the users have a small number of connections between each other, making clear the presence of unusual behaviors.

Fig. 1 .

 1 Fig. 1. Overview Steps of the Solution

  v2) < d and |time(v1) -time(v2)| < t 0, in any other case (2)

Fig. 2 .

 2 Fig. 2. Events graph

Fig. 3 .Fig. 4 .

 34 Fig. 3. Events per hour of the day

Fig. 5 .Fig. 6 .

 56 Fig. 5. Distribution of total distance of trips

Fig. 7 .Fig. 8 .

 78 Fig. 7. Sybil defender results modifying α values

Fig. 9 .Fig. 10 .

 910 Fig. 9. Largest component of driver speed graphs

Fig. 11 .

 11 Fig. 11. Location of users of traffic jam graph

www.facebook.com

www.instagram.com

www.twitter.com

www.waze.com

Acknowledgment

This work is partially supported by the University of Santiago research project PMI-USA 1204 and FONDEF Idea ID15I10560, CONICYT, Chile.