Portrait Quality Assessment using Multi-Scale CNN

In this paper, we propose a novel and standardized approach to the problem of camera-quality assessment on portrait scenes. Our goal is to evaluate the capacity of smartphone front cameras to preserve texture details on faces. We introduce a new portrait setup and an automated texture measurement. The setup includes two custom-built lifelike mannequin heads, shot in a controlled lab environment. The automated texture measurement includes a Region-of-interest (ROI) detection and a deep neural network. To this aim, we create a realistic mannequins database, which contains images from different cameras, shot in several lighting conditions. The ground-truth is based on a novel pairwise comparison technology where the scores are generated in terms of Just-Noticeable-differences (JND). In terms of methodology, we propose a Multi-Scale CNN architecture with random crop augmentation, to overcome overfitting and to get a low-level feature extraction. We validate our approach by comparing its performance with several baselines inspired by the Image Quality Assessment (IQA) literature.

Introduction
With the rapid rise of camera technology and social media photography, short videos and selfie photos gained a huge interest between users, especially young adults and teenagers. Most of these shots are captured with the front/selfie camera of a smartphone, requiring manufacturers to enhance the quality of their cameras and pushing forward the competition on the social media market. For that, brands compete to deliver the best rendering of portraits, which constitute the main subject of a selfie camera.
The quality of an image can be evaluated on several attributes: target exposure, dynamic range, color (saturation and white balance), texture, noise and different artifacts [1]. For portrait images, the main focus is to evaluate the render quality of the face. On this matter, characteristics like skin tone, bokeh, texture details and skin smoothness are mostly of interest. These attributes should be measured in a way that reflects the human perception. In this study, we aim at evaluating camera capabilities to preserve fine texture details on the face. This study falls within a more general topic that we refer to as portrait quality assessment.
When evaluating a camera device, the analysis usually focuses on its capabilities in low-light, zoom or shallow depth-offield simulation. One standard way to evaluate the quality difference between camera devices is to compare their output in a controlled environment and on the same visual content, namely a chart. Using a chart to evaluate the camera quality has both benefits and drawbacks over uncontrolled scenes. On the one hand, it ensures a repeatable measurement and ease of interpretation of the comparison. Since it provides a consistent visual content, environment variance and bias can be eliminated, thus associating the score inconsistency to the camera's capability alone. On the other hand, one can argue that charts are not very good in reflecting the human perception. In fact, lab experiments do not reflect real-world conditions and the visual content of the chart may differ from the content that the final user is interested in.
Cameras have been traditionally evaluated with techniques based on an explicit estimation of the Modulation Transfer Function (MTF) corresponding to the optical system. These approaches can be employed only in the case of synthetically generated visual charts. MTF-based methods suffer from important drawbacks. These methods were originally designed for optical systems that can be modeled as linear. Consequently, non-linear processing such as multi-image fusion or deep learning-based image enhancement, may lead to inaccurate quality evaluation [2]. Moreover, these methods assume that the norm of the device transfer function is a reliable measure of texture quality. However some recent works have shown that the magnitude of image transformations do not always coincide with the perceived impact of the transformations [3]. Therefore, we claim that comparing camera devices on a synthetic visual content is not sufficient to capture the complex behaviour of modern imaging systems.
We favour an explicit contribution of human judgment in the texture quality measurement process. On this matter, deep learning approaches have been proposed recently [4]. In this paper, the authors introduced a convolutional neural network that was trained on a custom database, composed of an arranged collection of objects with different textures and annotated by image quality experts. That said, the focus of this database is to simulate a natural scene environment and not portraits. Also, since it only contains charts, it has strict framing conditions, which does not respond to the natural framing of selfies and human portraits.
Responding to the limitations of the previous techniques, we follow a new and standardized approach for assessing portrait quality preservation. Our contributions are threefold. First, we introduce a novel laboratory setup based on custom-built lifelike mannequins to evaluate the portrait quality preservation of cameras. Using this setup, we establish a portrait database of mannequin images, shot in a controlled lab environment with nonstrict framing and annotated by image quality experts. That said, the data we use in the experiments are not available along with this paper. Second, we adopt a deep learning approach based on a multi-scale convolutional neural network (MSCNN) to measure the portrait quality. The scores should reflect the subjective quality judgment; thus, we employ a recent approach of generating subjective image quality scores based on pairwise comparisons to build our ground-truth. Finally, we conduct an extensive study that shows that our mannequin setup combined with the learningbased method can standardize portrait quality assessment of cameras on specific attributes, like details preservation, and performs better than existing methods for general purpose texture quality evaluation.

MTF-based methods
MTF-based methods suppose that the camera can be modeled as a linear system that produces an image y as a convolution of the point spread function h and the incoming radiant flux x. In the frequency domain, Y ( f ) = H( f )X( f ) + N( f ), where we also consider additive noise N. The modulation transfer function, MT F( f ) = |H( f )| is commonly used to characterize an optical acquisition device.
These MTF-based methods assume that the noise-free content, referred to as x, is available in order to estimate the transfer function of the system. This implies that they work the best with synthetic visual charts and cannot be employed for real-world images. Early methods use charts containing a blur spot or a slanted edge for this computation. In [5], Loebich et al. propose a method using the Siemens-Star. Cao et al. propose to use the Dead-Leaves model [6], and introduce an associated method in [7], which is shown to be more appropriate to evaluate fine detail preservation since the texture is more challenging for camera devices. In the rest of this paper, we refer to this chart as the Dead-Leaves chart (DL). In this chart, the reference image consists of occluding disks generated with a random center location, radius and grey-scale value.
Importantly, digital camera systems present high-frequency noise, which affects the MTF estimation by dominating signal in the higher frequencies. Consequently, estimating the noise power spectral density (PSD) is key to obtain an accurate acutance evaluation, and this task is not easily performed on the textured region. One approach is to estimate the noise PSD on a uniform patch with a known reference color, which then can be used to estimate the MTF. In modern digital cameras, this approach is affected by denoising algorithms, which not only interfere with the noise PSD, but also perform differently between textured and uniform regions. To bypass this limitation, other approaches [8,9] propose to compute the MTF using the crosspower spectral density between the reference and target image.
As a conclusion, MTF-based methods might give a good estimation of texture quality. However, it has been shown that MTF does not always reflect the human perception of qualilty. This observation pushes for more modern solutions, mainly learningbased methods that aim at estimating texture quality as perceived by human subjects.

Learning-based methods
In opposition to the MTF-based approach, learning-based methods usually require annotated datasets and can be separated into classical and deep learning approaches. The first automated approaches [10,11] were based on a combination of hand-crafted features and a regression technique such as support vector regression (SVR). However, these methods were surpassed later by convolutional neural networks (CNN) [12] that proved their superiority in image processing and computer vision applications. Consequently, we focus on deep learning solutions and divide them into two categories, based on their main usage: image quality oriented and camera quality oriented.

Image quality oriented methods
Early datasets (LIVE [13], CSIQ [14], and TIDs [15,16]) consist of noise-free images processed with several artificial distortions and annotated based on subjective preference. These distortions aim to describe compression or transmission scenarios and most of them are not relevant to the problem of camera evaluation. Besides, they fail to capture the complexity of modern camera systems with non-linear processing pipelines. Other datasets such as KonIQ10k [17] and LIVE In the Wild [18] consist of media-gathered images with unknown distortions. The latter claims that the images were captured using a representative variety of modern mobile devices. For both datasets, annotations were collected using subjective preference crowdsourcing, which we refer to as subjective image quality assessment.
These datasets are large enough to conduct a deep learning solution for image quality evaluation [19,20,21]. However, because of their wild nature, uncontrolled environment and their focus on evaluating any input image, they do not form a strong background to evaluate the quality of camera devices, which we are most interested in.

Camera quality oriented methods
Since evaluating camera quality focuses on fixed visual content comparisons, there is a need to define a different training setting, relying on specific datasets. Tworski et al. [4] introduced, to the best of our knowledge, the first of its kind deep learning camera quality evaluation solution. They adopt a regression formulation and train a network to estimate the camera capacity to preserve texture, by comparing images of a common perceptual chart taken with different devices. They introduce a mechanism to identify the chart regions that are the most suited to evaluate quality. Nevertheless, this solution assumes an effective registration of the chart, thus posing a risk of distorting the visual content. Besides, they try to estimate texture preservation quality on a simulated natural scene, which does not respond to our problem of evaluating texture quality on portraits.

Portrait quality assessment methods
Despite the lack of methods that solely focus on evaluating portrait quality preservation of cameras, there exist many attempts to evaluate face quality to help face detection, or even utilize portraits as a subject of image quality assessment. For instance, in [22], the authors propose a method to define a standard portrait image, which can be later used to evaluate colour rendering and other attributes between cross-media. In this work, we can see a first attempt to use a standard portrait as a subject for image quality assessment. Many other approaches to evaluate the face quality as a support for face recognition exist in the literature. In [23], Patrick Grother et al. from the National Institute of Standards and Technology (NIST), intend to support accurate face recognition by establishing specifications for face image quality assessment algorithms as well as evaluate the performance of these algorithms. Other papers [24,25] propose different methods to evaluate the face quality in the context of face recognition.
Although these solutions deal directly with the problem of face quality, they try to solve a distinct problem, which is not related to camera evaluation. That said, they can still be used in the future as a support for portrait quality preservation measurement, if we need to extend the protocol.
Responding to the limitations and lack of similar work in the literature, we introduce the first standardized texture preservation quality evaluation on portraits, based on deep learning.

Realistic mannequins
In this work, we introduce a new portrait quality database with the main focus on repeatability and reliability. We are primarily interested in texture quality, although it is possible to extend to other attributes like noise, skin tone, and sharpness. To ensure a repeatable measurement, it is important to have a consistent visual content, specifically a chart. However, using a chart similar to what has been used previously, as for example a printed pristine portrait image, will not capture complex behaviour of camera devices. This is especially true when rendering the face's shape and hair's texture, which are not well reflected in a printed image. For this reason, we employ a dataset of custom-built lifelike mannequins, referred to as realistic mannequins (RM).
The dataset contains images of silicone heads, shot at a fixed distance, in different lighting conditions, and with several smartphones' front camera. Apart from the strict distancing, the framing is loose and is not fixed to a strict angle, contrarily to the usual chart setups. This configuration helps to simulate a real self-portrait (selfie) scenario as shown in Figure 1a.

Ground truth & region of interest
Many attributes can be extracted from the realistic mannequins, in particular texture details preservation. Practically, to ensure a reliable measurement, we need to focus on a specific region of interest (ROI). In the case of texture measurement, we focus on the eye and eyebrows region, as shown in figure 1b. We choose the highest resolution in the dataset and upscale the other crops to match it, keeping the same aspect ratio.
To obtain a ground truth for the dataset, we need to collect subjective preference evaluations. However, evaluating the image quality in an absolute manner is not a simple task, and is still a hot topic in research. One of many approaches is pairwise comparison, where one can infer quality scores out of a comparison matrix. The main problem with this approach is its quadratic growth, which means that the cost and difficulty increase in a quadratic manner with the size of the dataset. Following ITU-T [26] and ITU-R [27] recommendations, a minimum of 15 full comparisons O(n 2 ) (i.e. 15 annotators to compare all n(n − 1)/2 pairs) is required to generate reliable results, which can be very costly and time consuming. Usually it is not possible to achieve a complete design (a full passage on all data points for each observer), which implies the need for a sampling strategy or active sampling as it's referred to in the literature.
We adopt the research of Mikhailiuk et al. [28] in order to efficiently utilize observers' time and obtain the most possible accurate scale. In this algorithm, we select the next comparison to deliver the most information, i.e., the one that has the highest impact on the posterior distribution of the scores. For that, separately for each image, we obtained the distribution of quality scores after every performed comparison using Expectation Propagation (EP), assuming Thurstone case V [29], and then estimated the distribution of the scores assuming every possible future comparison. The comparison maximizing the Kulback-Leibler divergence between the current distribution and the result after every possible comparison was chosen to be performed next.
The analysis and pairwise comparisons were conducted in a controlled environment with a fixed viewing condition. We adopt our setting so that our viewing condition is aligned with that of a human eye, with a cutoff frequency ν cut = 30cpd. Hence, we use a 32" 4k monitor with a pixel pitch of 0.185 , and we fix the eye to screen distance to 65cm. We use a calibrated display (D65 whitepoint with luminance ≥ 75cd/m 2 with no direct illumination of the screen and a background illumination with a lighting panel set to D65 / 15% for reducing eye stress.

Method
In this section, we detail the proposed method for estimating texture quality on portraits. We formulate this task as a regression problem and suppose that a perfect registration of the images is not needed.

ROI extraction
The first step in evaluating the texture quality on portraits is to correctly detect the ROI. Usually, the ROI is fixed and its detection assumes an effective registration of the chart. However, because of the inaccurate framing of portraits, registration is a bad choice, since it can introduce face distortions and alter the quality of the image. Therefore, using a manually extracted reference ROI, with a sufficiently high resolution, we could determine the crop dimensions in the other images. Thoroughly, we detect a set of facial landmarks to localise the eyes region. Then, we fix the width of the crop to cover the ROI horizontally and compute the height using the same reference aspect ratio. Finally, we upscale the crop to the reference resolution.

Batch creation with random crops
In our problem, the dataset is content-specific, meaning we have the same content in every image (same face). Consequently, there is a high possibility of overfitting when using a complete image. To avoid this, we use random crops as a data augmentation technique. We create a batch of N images, where every image is represented by n patches, assigned with the same score as the source image; the total batch size is then N * n, where each patch is considered as an independent image. Also, By using random crops, we can fix the input size without the need to resize the original image, which maintains the quality and solves the problem of limited memory. Moreover, with random crops there is a possibility of advanced extensions, as variance-based random crops and attention-based crops, as in [30,31]. Nonetheless, random crops can increase the variance and limit the training capacity, since some of the crops, especially skin crops, are not relevant and can alter the average score. To tackle this problem, we use the Huber loss as described in [4].

Multi-scale CNN
This paper aims to provide an automatic, all-in-one protocol to measure the portrait quality preservation capabilities of cameras. Therefore, we do not intend to introduce state-of-theart neural network architectures, but we prefer to adopt proven working solutions.
Inspired by the work of Chen et al. [32], we adopt a multiscale architecture for our convolutional neural network model, referred to as MSCNN. This technique aims to solve a problem with previous random patch methods, by paying equal attention to global and local features. The model consists of three convolutional blocks. Each convolutional block extracts specific features from one scale. Features are then concatenated into one vector and fed to a fully connected block (FC) with four layers. The final output is a single float value, representing the final quality value of the input list.
Since we have limited resources, we are not able to train a model that big from the ground up. For instance, the large parameter set can be easily trained on our relatively small dataset, however, the model will overfit because it is pretty complex for a small dataset. An intuitive solution to the problem is fine-tuning. We use the convolutional backbone of a well known CNN architecture that was pretrained on a large dataset, to extract relevant features from the image. The fine-tuning process significantly decreases the time required to train and optimize a new model, especially for a smaller task, as it is already trained on a larger dataset, and is assured to give robust results. Finally, since we consider random crops as augmentation of the data, we optimize the loss over individual patch scores instead of the average image score.

Experiments
In this section, we perform an extensive experimental study of the proposed solution. The study consists of two tests: first, we test the performance of the measure on a mannequin database, which we create separately from the training database, referred to as the ground-truth. Second, we compare the measure with an MTF-based general purpose texture quality measure. All the data are internally collected and might not be available outside of the working team.

Training details
We implemented our deep learning model with PyTorch, which was then trained on a 48 devices dataset of realistic mannequin images, shot in three lighting conditions, a total of 144 images. All the images are in fact ROI crops that were upscaled to the same reference resolution with a fixed aspect ratio, using bicubic upsampling. We choose the reference resolution as the highest crop resolution in our dataset. We adopt a transfer learning technique using a MobileNetV3 [33] pre-trained on Im-ageNet [34]. We employ Adam optimizer, with an initial learning rate of 1·10 −5 and a decay of 10% every 5 epochs for a total of 50 epochs. To limit the effects of noise, contrast, exposure, saturation and over sharpening on the texture score, we add 7 different augmentations that slightly alter these attributes. After visually checking many examples of the augmented data, we conclude that the small alteration will not change the texture quality, thus we keep the same score for each augmentation. In addition to that, we add horizontal and vertical flips as well as a grayscale variation. In total, our training database contains 1584 images.

Test datasets
To ensure that the new method performs better than previous solutions, we test its performance on a well-established ground truth, which was carefully annotated by human observers. Also, we compare the correlation of our scores with the MTF based methods; we aim to conclude whether the MTF-Based methods are sufficient to evaluate the texture preservation quality of a selfie camera, or do we absolutely need a solution that is inspired by and respects human observations.
As there is no well-established reference dataset for our problem, we collected and annotated our own data on two different charts/visual content.

New-Resolution
First, we create a reference dataset of realistic mannequin images, shot at a fixed distance, in two different lighting conditions. We collect data from 23 smartphone front cameras, a total of 46 images. This dataset is referred to as New-Resolution (NR). The procedure to collect and annotate the images is similar to what has been done for the training dataset (Section Dataset). We choose to separate the annotations of the RM and the NR datasets, in order to ensure that no hidden correlation exists between the scores. Additionally, no device has been shared between both datasets.
To ensure a good variability in our texture evaluation, we shoot mannequins at two distances: 30cm and 55cm. For the 30cm, we only provide female mannequins, while for the 55cm we only provide male mannequins with a consistent beard. To evaluate the texture, we use eyebrow crops for the female, and eyebrow and beard crops for the male. Finally, for each setup, we evaluate the models that were trained on its respective training set, as well as evaluate the average outputs of both setups.

Gray-DL
Second, we employ the Dead-Leaves chart proposed in [7] (Section MTF-based methods). In our experiment, we refer to this dataset as Gray-DL. We use the same two lighting conditions and 23 devices as in the New Resolution dataset. In the case of the Dead-Leaves charts, since the charts are unnatural images, human perceptual annotation is problematic. Since we do not have a proper ground truth for this chart, we choose to evaluate it on the ground truth of the New Resolution dataset, similarly to what we do on our deep learning model.
Nevertheless, there exists one limitation that concerns us when assigning the NR ground truth to the Gray-DL dataset, which we try to explain and solve next. In fact, the Dead-Leaves chart is shot with a specific framing, and depending on the camera's field of view (FOV), the lens to chart distance can lie between 45cm and 55cm. Although the DL chart is theoretically scale invariant and thus the measurement should not depend on the distance, the device's range of focus will drastically alter the results. For instance, one device might feature an Autofocus (AF), while the other might have a fixed range of focus. Usually, this range is between 30 and 60cm, but can be freely positioned by the manufacturer. For this reason, only the 55cm distance is relevant to evaluate the Gray-DL scores. To solve this problem, we choose to test on two different device sets. The first set includes all 23 devices (46 images) for the 55cm distance, referred to as the Full set, while the second only includes 13 devices (26 images) that have AF and/or a wide range of focus, covering both 30cm and 55cm distances; we refer to this set as the Limited set.

Quantitative study Metrics
Since we adopt a regression formulation, we need to use relevant metrics that comply with an image quality score. That said, this score is not absolute and is usually relative to the approach of origin. Consequently, MTF-based methods predict a quality score that is not directly comparable to the score provided by human annotators. Also, since we annotate the RM and NR datasets separately, their scores lie on separate scales, which makes a direct comparison irrelavant. One alternative is to explore the linear correlation between the output and the ground truth. However, assuming that the predictions correlate linearly with the ground truth may not hold and bias the experiment. Therefore, we adopt two metrics that are commonly used in the IQA literature, and that are based on the correlation of the rank-order. First, we adopt the Spearman Rank-Order Correlation Coefficient (SROCC) defined as the linear correlation coefficient of the ranks of between two variables. Second, we report the Kendall Rank-Order Correlation Coefficient (KROCC) defined by the difference between concordant and discordant pairs divided by the number of possible pairs. The key advantage of the second metric lies in its robustness to outliers.

Results
We compare the MTF-based measurement on the Dead-Leaves chart with the ground truth of the New-Resolution dataset for the 30cm and 55cm distances, as well as for the average score over both datasets. For the realistic mannequin texture evaluation using the MSCNN, we compare the output of each model with its respective NR dataset. Additionally, we compare the average of both outputs to the average score of the NR datasets. All the results are shown in Tables 1, 2.
As we explained in Section Gray-DL, we proceed with two different evaluations based on the focus characteristics of each device. To ensure a correct ground truth, we compute the correlation between the NR-55cm and NR-30cm for the Limited set. If the value of this correlation is large enough, we can suppose that we have the same device ranking regardless of the distance. For the Limited set with 13 devices, we get srocc = 0.935 & krocc = 0.796, which ensures a similar ranking across both distances.
On both device sets we notice that the MSCNN based methods give larger correlations over the MTF-based methods. Moreover, we can see that the model performs better on the 30cm distance. That said, we should take the results with a grain of salt. Since the Limited set's size is relatively small, we cannot extract a statistically significant conclusion. However, when we back these results up with the Full set, we have a higher significance to support our new method. Because of the fixed range of focus, the evaluation of texture can change drastically between distances, which explains why we do not include the comparison for the 30cm distance on the Full set. However, in average, over the two distances, we can get an idea of the texture preservation capabilities of a smartphone front camera, which is expressed here in the NR-AVG scores.
Because of the ever-growing image quality enhancement technologies, new camera systems incorporate specific pipelines for certain scenes, especially portraits. One advantage of the RM-based approach is that a mannequin face triggers these pipelines, which helps to correctly assess the capabilities of the selfie camera. This scenario is not possible with a normal chart. Therefore, we see that the RM-based method can give a more accurate ranking and a better evaluation of the texture preservation quality of a selfie camera, where portraits are the main subject.

Conclusion
In this paper, we proposed a novel and standardized approach to the problem of camera quality assessment on portrait scenes. We developed a repeatable setup as well as a deep learning-based method to test the texture preservation quality on portraits, inspired by human observations. Our results show that a deep learning approach, can outperform the classical measurements and is more relevant to human preference. One limitation of our method is that it only works on specific mannequins, and not on real human portraits. Consequently, as future work, we plan to extend our database to cover human portraits, and multiple human portraits, in order to highlight several complementary discriminant features and better measure the intrinsic qualities of a smartphone front camera, in a real-life scenario.


Distributed Evaluation of Graph Queries using Recursive Relational Algebra

We present a system called Dist-$\mu$-RA for the distributed evaluation of recursive graph queries. Dist-$\mu$-RA builds on the recursive relational algebra and extends it with evaluation plans suited for the distributed setting. The goal is to offer expressivity for high-level queries while providing efficiency at scale and reducing communication costs. Experimental results on both real and synthetic graphs show the effectiveness of the proposed approach compared to existing systems.

I. INTRODUCTION
With the proliferation of large scale graphs in various domains (such as knowledge representation, social networks, transportation, biology, property graphs, etc.), the need for efficiently extracting information from these graphs becomes increasingly important. This often requires the development of methods for effectively distributing both data and computations so as to enable scalability. Efforts to address these challenges over the past few years have led to various systems such as MapReduce [1], Dryad [2], Spark [3], Flink [4] and more specialized graph systems like Google Pregel [5], Giraph [6] and Spark Graphx [7]. While these systems can handle large amounts of data and allow users to write a broad range of applications, they still require significant programmer expertise. The system programming paradigms and its underlying configuration tuning must be highly mastered. This includes for example figuring out how to (re)partition data on the cluster, when to broadcast data, in which order to apply operations for reducing data transfers between nodes of the cluster, as well as other platform-specific performance tuning techniques [8].
To facilitate large-scale graph querying, it is important to relieve users from having to worry about optimization in the distributed setting, so that they can focus only on formulating domain-specific queries in a declarative manner. A possible approach is to have an intermediate representation of queries (e.g. an algebra) in which high level queries are translated so that they can be optimized automatically. Relational Algebra (RA) is such an intermediate representation that has benefited from decades of research, in particular on algebraic rewriting rules in order to compute efficient query evaluation plans.
A very important feature of graph queries is recursion, which enables to express complex navigation patterns to extract This research has been partially supported by the ANR project CLEAR(ANR-16-CE25-0010). useful information based on connectivity from the graph. For instance, recursion is crucial for supporting queries based on transitive closures. Recursive queries on large-scale graphs can be very costly or even infeasible. This is due to a large combinatorial of basic computations induced by both the query and the graph topology. Recursive queries can generate intermediate results that are orders of magnitude larger than the size of the initial graph. For example, a query on a graph of millions of nodes can generate billions of intermediate results.
Therefore being able to optimize queries and reduce the size of intermediate results as much as possible becomes crucial.
Several works have addressed the problem of query optimization in the presence of recursion, in particular with extensions of Relational Algebra [9]- [11]; and with Datalog-based approaches [12] such as BigDatalog [13]. Recently, µ-RA [11] proposed logical optimization rules for recursion not supported by earlier approaches. In particular, these rules include the merging and reversal of recursions that cannot be done neither with Magic sets nor with Demand Transformations that constitute the core of optimizations in Datalog-based systems [11]. However, µ-RA is limited to the centralized setting.
In this paper, we present Dist-µ-RA a new method and its implementation for the optimized distributed evaluation of recursive relational algebra terms. Specifically, our contribution is threefold: 1) a new method for the optimization of distributed evaluation of queries written in recursive relational algebra.
Since it uses a general recursive relational algebra, it can be of interest for a large number of mainstream RDBMS implementations; and it can also provide the support for distributed evaluation of recursive graph query languages. For example Dist-µ-RA provides a frontend where the programmer can formulate queries known as UCRPQs [14]- [17] 1 ). This method provides a systematic parallelisation technique by means of physical plan generation and selection. These plans automatically repartition data in order to reduce data transfer between cluster nodes and communication costs during recursive computations.
2) a classification of graph queries by the means of six query classes. Each class characterizes queries with a particular feature: for example a recursion with a filter, or concatenated recursions. A fundamental idea of this classification is that a query may belong to one or more classes. Whenever a query belongs to several classes this means that it requires the optimization techniques of all the corresponding classes, together with a technique capable of combining them. Therefore, the more classes a query belongs two, the harder is its optimization. 3) a prototype implementation [18] of the system on top of Apache Spark. Specifically, Dist-µ-RA can use plain Apache Spark or Apache Spark with PostgreSQL as a DBMS backend. Dist-µ-RA is experimentally evaluated using queries covering the different classes, and using datasets (both real and synthetic) of various sizes. Experimental results show that Dist-µ-RA is more efficient than state-of-the-art systems such as BigDatalog [13] in all query classes except one where their performance is similar. This indicates that Dist-µ-RA represents a better alternative for the distributed evaluation of recursive graph queries.
The outline of the paper is as follows: we first present a summary of useful preliminary notions in Section II. Section III presents the whole architecture of Dist-µ-RA. In Section IV, we show how µ-RA terms are distributed and how physical plans are generated. Experimental evaluation is presented in Section V and related works are discussed in Section VI.

II. PRELIMINARIES
The µ-RA algebra [11] is an extension of the Codd's relational algebra with a recursive operator whose aim is to support recursive terms and transform them when seeking efficient evaluation plans. The syntax of µ-RA is recalled from [11] in Fig. 1. It is composed of database relation variables and operations (like join and filter) that are applied on relational tables to yield other relational tables. µ is the fixpoint operator. In µ(X = Ψ), X is called the recursive variable of the fixpoint term. Like in RA, the data model in µ-RA consists of relations that are sets of tuples which associate column names to values. For instance, the tuple { src → 1, dst → 2 } is a member of the relation S of Fig. 2.
Let us consider a directed and rooted graph G, a relation E that represents the edges in G, and a relation S of starting edges (a subset of edges in E that start from the graph root nodes), as represented in Fig. 2.
The following examples illustrate how µ-RA algebraic terms can be used to model graph operations, such as navigating through a sequence of edges in a graph: Example 1. The term π c (ρ c dst (S) ρ c src (E)) returns pairs of nodes that are connected by a path of length 2 where the first element of the pair is a graph root node. For that purpose, the relation S is joined ( ) with the relation E on the common column c, after proper renaming (ρ) to ensure that c represents both the target node of S and the source node of E. After the join, the column c is discarded by the anti-projection ( π c ) so as to keep only the two columns src, dst in the result relation.
Example 2. Now, the recursive term µ(X = S ∪ π c (ρ c dst (X) ρ c src (E))) computes the pairs of nodes that are connected by a path in G starting from edges in S.
The subterm ϕ = π c (ρ c dst (X) ρ c src (E)) computes new paths by joining X (the previous paths) and E such that the destinations of X are equal to the sources of E.
The fixpoint is computed in 4 steps where X i denotes the value of the recursive variable at step i:

B. Semantics and properties of the fixpoint
The semantics of a µ-RA term is defined by the relation obtained after substituting the free variables in the term (like E and S in example 2) by their corresponding database relations. The notions of free and bound variables and substitution are formally defined in [11].
As an slight abuse of notation, we sometimes use a recursive term Ψ (i.e. a term that contains a recursive variable X) as a function R → Ψ(R) that takes a relation R and returns the relation obtained by replacing X in the term Ψ by the relation R. In the above example Under this notation, µ(X = Ψ) is defined as the fixpoint F of the function Ψ, so Ψ(F ) = F . Let us consider the following conditions (denoted F cond ) for a fixpoint term µ(X = Ψ) (as also considered in [11]).
• positive: for all subterms ϕ 1 ϕ 2 of Ψ, ϕ 2 is constant in X (i.e. X does not appear in ϕ 2 ); • linear: for all subterms of Ψ of the form ϕ 1 ϕ 2 or ϕ 1 ϕ 2 , either ϕ 1 or ϕ 2 is constant in X; • non mutually recursive: when there exists a subterm µ(Y = ψ) in Ψ, then any occurence of X in this subterm should be inside a term of the form µ(X = γ).
Proposition 2. Every fixpoint term µ(X = Ψ) that satisfies F cond can be written like the following: µ(X = R ∪ ϕ) where R is constant in X and ϕ(∅) = ∅. R is called the constant part of the fixpoint and ϕ the variable part.
In Example 2, S is the constant part and π c (ρ c dst (X) ρ c src (E)) is the variable part. In the rest of the paper, we only consider fixpoint terms satisfiying the conditions F cond in their decomposed form µ(X = R ∪ ϕ) since their existence is guaranteed thanks to proposition 1.

C. Evaluation of the fixpoint
A fixpoint term can be evaluated with the algorithm: The fixpoint is obtained by evaluating ϕ repeatedly starting from X = R until the fixpoint is reached. In this algorithm, we apply ϕ on the new results only (obtained by making a set difference between the current result and the previous one) instead of the entire result set. This is possible thanks to the property of ϕ stated in proposition 1, which implies that Evaluating recursive computations on the new iteration results is well known in the context of Datalog [12] and transitive closure evaluation [19] with the semi-naive (or differential) approach (see Section VI on related works for more discussion on these aspects).

III. ARCHITECTURE OF THE DIST-µ-RA SYSTEM
The Dist-µ-RA system takes a query as input parameter, translates it into µ-RA, optimizes it, and then performs the evaluation in a distributed fashion on top of Spark. For this purpose the Dist-µ-RA system is composed of several components, as illustrated in Fig. 3.

PhysicalPlanGenerator
PgSQLExecutor SparkExecutor Spark PgSQL The Query2Mu component translates UCRPQs into µ-RA terms. UCRPQ syntax is given in [11] and we give an example below. For sure, Dist-µ-RA also accepts more general µ-RAterms that are not expressible as UCRPQs 2 , as long as they satisfy the triple condition F cond mentioned in Section II.
From a given input µ-RA term, the MuRewriter explores the space of semantically equivalent logical plans by applying a number of rewrite rules. In addition to the rewrite rules already known in classical relational algebra, MuRewriter applies a set of rules specific to the fixpoint operator. These rules and the conditions under which they are applicable are formally defined in [11].
The evaluation cost of these terms are estimated by the CostEstimator component, based on data cardinality estimations proposed in [20], and that outputs a best estimated logical plan for recursion. From a given recursive logical plan, the PhysicalPlanGenerator generates a physical plan for distributed execution. It can also generate a centralized plan for PostgreSQL, for the dual purpose of (i) efficiently executing certain subqueries in a centralized manner and for (ii) performance comparison with the centralized case, in particular with [11].
The PhysicalPlanGenerator is an original contribution of this paper and it is described in more detail in Section IV.
We describe below steps that sample queries go through for evaluation in the Dist-µ-RA system. Consider for instance the following UCRPQ composed of a conjunction of two RPQs: ?a,?b,?c ←?a wasBornIn/IsLocatedIn+ Japan, The first RPQ computes the people ?a that are born in a place that is located directly or indirectly in Japan. The query is first translated into µ-RA by Query2Mu so that MuRewriter can generate semantically equivalent plans. We describe below the rewrite rules specific to fixpoint terms leveraged from [11] that can apply in MuRewriter, and we give the intuition of their effect on performance: • Pushing filters into fixpoints: with this rule, the query ?x isLocatedIn+ Japan is evaluated as a fixpoint starting from ?x such as ?x isLocatedIn Japan, which avoids the computation of the whole isLocatedIn+ relation followed by the filter Japan. • Pushing joins into fixpoints: let us consider the query ?x isMarriedTo/knows+ ?y. Instead of computing the relation knows+ and joining it with isMarriedTo, this rule rewrites the fixpoint such that it starts from ?x and ?y that verify ?x isMarriedTo/knows ?y.
The application of this rule is beneficial in this case because the size of the isMarriedTo/knows relation is usually smaller than the size of the knows relation. • Merging fixpoints: instead of computing both fixpoints separately then joining them, this rule generates a single fixpoint that starts with isLocated/dealsWith then recursively appends either isLocatedIn to the left or dealsWith to the right. • Pushing antiprojections into fixpoints: this rule gets rid of unused columns during the fixpoint computations. For instance, the query ?y ?x isLocatedIn+ ?y (this query asks for ?y only) is evaluated by starting only from the destinations ?y of the isLocatedIn relation and by recursively getting the new destinations, thus avoiding to keep the pairs of nodes ?x and ?y then discarding ?x at the end. • Reversing a fixpoint: the fixpoint corresponding to the relation a+ can either be computed from left to right by starting from a and by recursively appending a to the right of the previously found results, or from right to left by starting from a and appending a to the left. Reversing a fixpoint consists in rewriting from the first form to the other or vice versa. This rule is necessary to account for all possible filters and joins that can be pushed in a fixpoint. For instance, a filter that is located at the left side of a+ can only be pushed if the fixpoint is evaluated from left to right.
After these transformations, the best (estimated) recursive logical plan selected by CostEstimator is given as input parameter to PhysicalPlanGenerator that is in charge of generating the best physical plans for distributed execution, and which we now further describe in the next section.

IV. DISTRIBUTED EVALUATION
Non-recursive µ-RA expressions are directly translated into Spark [3] operations using the Dataset API. Datasets are used to represent relational data in Spark. The optimization of these expressions is delegated to Spark's Catalyst internal optimizer [21] before execution. We now describe how fixpoint terms are evaluated in a distributed manner, first by explaining the principles and then how physical plans are generated.

A. Principles of the distributed evaluation of the fixpoint
In this section we present two ways of distributing the fixpoint computation on a Spark cluster: 1) Global Loop on the Driver (P gld ): In order to execute the fixpoint in a distributed setting, a simple way is to distribute the computations performed at each iteration of the Algorithm 1. Fig. 4 (left side) illustrates this execution in Spark.The so-called driver 3 performs the loop and, at each iteration, instructions at lines 4 and 5 are executed as Dataset operations that are distributed among the workers. We call this execution plan P gld . On Spark, ∪ is executed as Dataset union followed by a distinct operation. This means that in P gld , at least one data transfer (shuffle) per iteration is made to perform the union.
2) Parallel Local loops on the Workers (P plw ): Another way to distribute the fixpoint is based on the following observation: Under the conditions F cond , we have: which is a consequence of proposition 1. This proposition means that a fixpoint whose constant part is a union of two datasets can be obtained by making the union of two fixpoints, each with one of these datasets as a constant part. Thanks to this property 3, the fixpoint can be executed by distributing Figure 4. P gld (left) and P plw (right) execution on the cluster. Colored arrows show data transfers that occur at each iteration of the fixpoint. the constant part R among the workers, then each worker i executes a smaller fixpoint µ(X = R i ∪ ϕ) locally starting from its own constant part R i . We call this execution plan P plw . Execution is illustrated in the right side of Fig. 4. As opposed to P gld , P plw incurs only one data shuffle at the end to make the distinct union between the local fixpoints.
In Example 2, if we split the start edges S among two workers by giving the first (1, 2) and (10, 11) and the second (1, 4) and (10, 13), after executing P plw , the first worker will find the paths (1, 3), (10,5), (10,6) and (10,12) and the second will find (1, 5), (1,6), and (10, 12). Data distribution for P plw : there are cases where the final data shuffle incurred by P plw can also be avoided by properly repartitioning input data among workers. We first give the inuition, then we give the proof below.
We look for a column col (or a set of columns) in X that is left unchanged by ϕ. In other words, a tuple in R having a value v at column col will only generate tuples having the same value at this column throughout the iterations of the fixpoint. So if we put all tuples in R having v at column col in one worker, no other worker will generate a tuple with this value at that column.
For this, we use the stabilizer technique defined in Definition 10 of [11] and used to push filters in fixpoint expressions. It consists of computing "the set of columns which are not altered during the fixpoint iteration". For instance, 'src' is a stable column in the fixpoint expression of example 2 meaning that tuples in the fixpoint having 'src' = 1 can only be produced from tuples in S having 'src' = 1, which implies that filtering tuples having 'src' = 1 before or after the fixpoint computation lead to the same results. However, this is not true for the column 'dst' which is not stable.
To summarize, when the constant part of the fixpoint is repartitioned by the stable column (or columns) prior to the fixpoint execution, we know that there will be no duplicate across the workers (so we can avoid calling distinct at the end of the computation).
Proof. Let c be a stable column of µ(X = R ∪ ϕ), which means that ∀e ∈ µ(X = R ∪ ϕ) ∃r ∈ R e(c) = r(c) [11]. In µ-RA, an element r in R is a mapping (tuple), which means that it is a function that takes a column name and returns the value that r has at that column.
Let us consider a partitioning R 1 , ..., R n of R by the column c which verifies the following This statement means that there are no two elements of R at different partitions that share the same value at column c. We next show that this statement is also true for the fixpoint term.
In conclusion, the sets µ(

B. Physical plan generation
We first present the different physical plans for the µ-RA operators, and then we explain how they are selected. a) The fixpoint operator: As mentioned earlier, there are two different ways of evaluating the fixpoint operator in a distributed manner. P plw has two implementations (physical plan alternatives): • P pg plw : The local fixpoints are executed on PostgreSQL. The fixpoint operator is performed as a Spark mapPartition operation where each worker performs a part of the fixpoint computation on PostgreSQL. A PostgreSQL instance runs on each worker. The part of data assigned to each worker is represented as a view in the PostgreSQL instance running on this worker, and the µ-RA expression (that computes the fixpoint) is translated to a PostgreSQL query that is executed using this view as the constant part of the fixpoint. Each PostgreSQL executor returns its result as an iterator which is then processed by Spark. • P s plw : The fixpoint computation is implemented using a loop in the driver that uses Spark operations to compute the recursive part of the fixpoint. These Spark operations are written in such a way each worker performs its own fixpoint independently, so without data being sent across the cluster. For example, to perform a join, the Spark broadcast join is used. To perform the union (or setdifference), a special union (set-difference) operation is used that computes the union (set-difference) partitionwise. These special union and set-difference operations are implemented as part of the SetRDD API. SetRDD is a special RDD 4 where each partition is a set. This SetRDD is used to store the value of the recursive variable X at each iteration. This means that each partition of X holds the intermediate results of the local fixpoint performed by the worker to which this partition has been assigned. This SetRDD technique is inspired from BigDatalog [13]. P gld is implemented using a global loop in the driver. At each step of the recursion a distinct operation is performed. Intermediate data is stored in a Spark Dataset. b) Other operators: An operator in µ-RA has a local version implemented using PostgreSQL and a distributed version implemented using the Spark Dataset API. Some operators may have more than one distributed execution plan. For instance, for the join operator, we choose which argument (if any) to broadcast in order to guide Spark on whether to use broadcast join or another type of join. c) Plan selection: while generating the physical plan for the fixpoint operator, we check whether there is a stable column (Sec. IV-A2). If so, we first shuffle data according to this column and execute the plan P plw . If not, the P gld plan is selected.
When P s plw is chosen, joins are executed as broadcast joins (all relations in the variable part of the fixpoint apart from the recursive relation are broadcasted), antiprojections are executed without the need of applying the distinct operation. If P pg plw is chosen, then the operators in the fixpoint expression are executed using the local PostgreSQL instances. 4 RDD is an abstraction that Spark provides to represent a distributed collection of data. An RDD is split among partitions which are assigned to workers.

V. EXPERIMENTS
We evaluate the performance of a prototype implementation of the Dist-µ-RA system on top of the Spark platform. We extensively compared the performance against other state-ofthe-art systems on various datasets and queries. We report below on these experiments.

A. Experimental setup
Experiments have been conducted on a Spark cluster composed of 4 machines (hence using 4 workers, one on each machine, and the driver on one of them We use real and synthetic datasets of different sizes and topological properties, as summarized in Table I. We consider the following real graphs: • Yago 5 : A knowledge graph extracted from Wikipidia and other sources [22]. • The Facebook, Reddit, DBLP and Live Journal real word datasets from the Snap network dataset collection [23].
We consider the synthetic graphs defined as follows: • rnd_n_p: random graphs generated with the Erdos Renyi algorithm, where n is the number of nodes in the graph and p the probability that two nodes are connected.
• uniprot_n: a benchmark graph of n nodes generated using the gMark benchmark tool [25]. It models the Uniprot database of proteins [26]. • tree_n: a random tree of n nodes generated recursively as follows: tree_1 is a tree of 1 node, and then tree_i + 1 is a tree of i + 1 nodes where the i th + 1 node is connected as a child of a randomly selected node in tree_i. • Other graphs derived from rnd_p_n by adding a set of predefined labels randomly. They are used to compute terms that require labeled graphs like the a n b n query and the concatenated closure queries (Sec. V-D).

C. Systems
We compare Dist-µ-RA with the following systems: • BigDatalog [13] available at [27]: a large-scale distributed Datalog engine built on top of Spark. • GraphX [7]: a Spark library for graph computations. It exposes the Pregel API for recursive computations. In order to compare our system with GraphX we need to convert UCRPQs to GraphX programs 6 . Specifically, we compute a regular graph query by making each node send a message to its adequate neighbors in such a way that the query pattern is traversed recursively from left to right. This means that for a query that starts by selection (?x ← A pattern ?x), only the node A sends a message at the start of the computation. • Myria [28]: A distributed big data management system that supports Datalog and exposes MyriaL, a query language similar to Datalog. We were not able to run this system on our cluster because the software stack has evolved and Myria is not maintained anymore. However, we were able to test it on a local configuration of four workers on a single machine. For this reason, comparison against Myria is performed on this local configurations on smaller test files and on a subset of test programs. • Centralized µ-RA [11]: Implementation of µ-RA on PostgreSQL that runs locally on a single machine.

D. Queries
Knowledge graphs like Uniprot and Yago are queried using UCRPQs. We also provide additional experiments using more general µ-RA terms not expressible by UCRPQs because they are not regular (e.g. same generation example).
Queries found in practice may contain various forms of recursion. To ensure that tested queries cover all forms of recursion, we rely on a classification of recursive queries in six classes: C 1 − C 6 . Each class regroups queries with a particular recursive feature as follows, where we give for each class an example in UCRPQ: 6 In the GraphX framework, a recursive computation is composed of "supersteps" where, in each superstep, graph nodes send messages to their neighbor nodes, then a merge function aggregates messages per recipient and each recipient receives its aggregated messages in order to process them. A computation is stopped when no new message is sent.
• C 1 : Single recursion, e.g. ?x, ?y ← ?x a+ ?y • C 2 : Filter to the right of a recursion, e.g. ?x ← ?x a+ C • C 3 : Filter to the left a recursion, e.g. ?x ← C a+ ?x • C 4 : Concatenation of a non recursive term to the right of a recursion, e.g. ?x, ?y ← ?x a+/b ?y • C 5 : Concatenation of a non recursive term to the left of a recursion, e.g. ?x, ?y ← ?x b/a+ ?y • C 6 : Concatenation of recursions, e.g. ?x, ?y ← ?x a+/b+ ?y Each class requires specific optimizations. For instance, the optimization of queries of classes C 2 and C 3 requires pushing a filter in the fixpoint (in two different directions). Queries of classes C 4 and C 5 can be subject of a rewrite by pushing a join in the fixpoint. C 2 and C 4 require reversing the fixpoint before applying rewritings. Queries of C 6 can be rewritten by merging fixpoints or by pushing a join in a fixpoint.
A query may belong to one or more classes. Whenever a query belongs to several classes this means that it requires the optimization techniques of all the corresponding classes, together with a technique capable of combining them. Therefore, the more classes a query belongs two, the harder is its optimization. For example, the query ?x ← C a/b+ ?x belongs to C 3 because there is a filter to the left of the recursion b+ and also belongs to C 5 because there is a concatenation to the left of the recursion b+.
To ensure that a wide range of recursive queries is covered in the experiments (see Figures 5 and 6), there is, for each class C i , at least a query that belongs to C i alone. The rest of the tested queries belong to different combinations of these classes. a) Yago queries: Fig. 5 lists the UCRPQs evaluated on the Yago dataset along with their classes. Queries Q 1 to Q 7 are taken from [29], Q 8 , Q 9 from [30], and Q 10 , Q 11 from [31]. We have added the rest of the queries to illustrate larger transitive closures.  Figure 5. Queries for the YAGO dataset 7 . 7 "isL" stands for "IsLocatedIn", "dw" for "dealsWith", "haa" for "hasAcademicAdvisor", "SA" for "Shannon_Airport", "JLT" for "John_Lawrence_Toole", and "wce" for "wikicat_Capitals_in_Europe". b) Concatenated closures: are queries of the form a 1 +/a 2 +/.../a n + where 2 ≤ n ≤ 10. They all belong to C 6 . c) µ-RA queries: the tested µ-RA terms are the following: • a n b n : computes the pairs of nodes connected by a path composed of a number of edges labeled a followed by the same number of edges labeled b. It is expressed by the following µ-RA term: µ(X = πm(ρ m trg (σ pred=a (R)) ρ m src (σ pred=b (R))) ∪ πm( πn(ρ m trg (σ pred=a (R)) ρ n trg (ρ m src (X)) ρ n src (σ pred=b (R))))) • Same Generation: computes the pairs of nodes that are of the same generation in a graph. The graph edges R represent the parent relation. We use the following term to express it: µ(X = πm(ρ m src (R) ρ m src (R)) ∪ πm( πn(ρ m src (R) ρ n trg (ρ m src (X)) ρ n src (R)))) • Reach: computes the reachable nodes in a graph from a source node N (chosen randomly from each test file).
All of these queries belong to C 1 . d) Uniprot queries: In order to test the comparative behavior of Dist-µ-RA on the Uniprot benchmark graphs Each tested query ( fig. 6) has at least one of those shapes. For instance, the query ?x,?y ← ?x interacts/(encodedOn/-encodedOn)+ ?y has the shape ?x,?y ← ?x a/b+ ?y.

E. Results
We report and comment on experimental results. Absence of a time in a figure means that the query evaluation has failed (the system has crashed). 1) Evaluation of UCRPQs on Yago:  Fig. 7 show that the implementation using SetRdd is faster. Therefore, we use this implementation as a basis of comparisons in the rest of the experiments. b) Comparison with other systems: Fig. 9 shows the execution times. The comparison between Centralized µ-RA and Distµ-RA shows that 17 out of the 25 tested queries execute faster when distributed. The 17 queries have all in common a large size of intermediate results. Queries handling large amounts of data tend to execute faster using Dist-µ-RA.
In Dist-µ-RA with P gld , all fixpoints in the queries are executed with the P gld plan. Comparison with Dist-µ-RA shows the impact of communication cost reduction performed by P plw . As explained in Sec. IV, P gld requires communications between the workers at each step of the recursion while P plw does not.
We notice that all the queries where Dist-µ-RA is significantly faster than BigDatalog (Q 9 and Q 22 to Q 25 ) have a transitive closure of size exceeding 20M records and belong to classes other than C 1 . Q 9 for instance belongs to C 2 , so it requires reversing the fixpoint first before pushing the filter "Kevin Bacon". This fixpoint reversal is not supported by the Magic Sets optimization technique of Datalog (see Sec. VI for more details). A similar observation can be made for Q 25 of class C 6 and where two fixpoints are merged in Dist-µ-RA. The other queries either have transitive closures of sizes less than 20M or belong to C 1 like Q 21 . For these queries, Dist-µ-RA and BigDatalog times are very close. This confirms that Dist-µ-RA performs more optimizations and shows that their impact on performance is significant when the size of the data processed by the query is large.
GraphX is much slower overall except for queries like Q 10 and Q 24 where filtering is performed at the beginning of the query (as mentioned in Sec. V-C). We believe that this lack of performance is due to the fact that, in the GraphX Pregel model, each node has to keep track of its ancestors that satisfy a given regular path query (or a part of it) and transmit this information to their successors in order to get the pairs of nodes satisfying the whole query. So while GraphX has been proven to be efficient for a number of graph algorithms [7], it can be less suitable for this kind of queries.
2) Concatenated closures: We evaluate concatenated closure queries on the graph obtained from rnd_100k_0.001 by randomly assigning 10 different labels to its edges. Results are shown in Fig. 10. Dist-µ-RA is faster on all queries (which all belong to C 6 ). The time difference between Dist-µ-RA and other systems for a query with n concatenations (a 1 +/.../a n +) gets larger when n increases. BigDatalog fails for queries where n ≥ 5 and Centralized µ-RA times out for all queries, while GraphX fails (crashed) on all queries. The plans that were selected in Dist-µ-RA for the execution of these queries apply a mixture of the rewritings that "push joins" and "merge fixpoints" (see Sec. III). The observed results show significant improvements brought by these rewritings on the performance of class C 6 queries.
3) µ-RA queries: Execution times are shown in Fig. 11. Distµ-RA and BigDatalog have comparable execution times. This shows that for a variety of queries in C 1 , Dist-µ-RA and Bigdatalog have a similar performance. The comparison between Myria and Dist-µ-RA is done on the same generation query. It can be observed in Fig. 12 that the relative perfomance of Dist-µ-RA and Myria increases with the dataset size in favor of Dist-µ-RA. Myria even crashes for rnd_10k_0.001. 4) Evaluation of UCRPQs on Uniprot: the execution times (on the cluster) are shown in Fig 13. Dist-µ-RA answers all of the queries and is faster on all queries belonging to C 2−6 (apart from Q 42 where the size of the transitive closure is small). This confirms the results observed for the Yago dataset (Sec. V-E1).
The comparison with Myria (on a single machine and on a smaller Uniprot file) is shown in Fig 14. We observe that all queries where Myria is slower than Dist-µ-RA (or where it fails) have a transitive closure of size exceeding 500,000 records (and none of the other queries contain a transitive closure greater than this size). Execution times on all other queries are very similar except for Q 42 . Note that Dist-µ-RA spends approximatey 10s to set the Spark context (and this time is included in the measured Dist-µ-RA execution time).
We perform further scalability tests by measuring Dist-µ-RA and BigDatalog execution times for each Uniprot query on the generated uniprot_n graphs with varying sizes of 1M, 5M and 10M edges. Fig. 8 shows the execution times for Dist-µ-RA in comparison to BigDatalog. BigDatalog fails in 44 cases out of a total of 75 query evaluations. Dist-µ-RA answers all of them and scales better.
For more comprehensive testing, queries and graph sizes have been selected so as to cover a wide range of result sizes. Q 40 is one of the queries with the smallest result size (14K for uniprot_10M) and Q 46 is one with the largest (around 1.5B for uniprot_10M, which is 150 times the size of the graph).

F. Summary
In summary, when it comes to queries involving large amounts of data, in particular as intermediate results, Dist-µ-RA is more efficient than centralized µ-RA. In addition, Dist-µ-RA is significantly more efficient when compared to GraphX

VI. RELATED WORKS
In order to evaluate expressive queries (such as UCRPQs) over graphs, it is essential for a system to be able to: (i) support recursion and the optimization of recursive terms, and (ii) provide distribution of data and computation. We examine and compare to the closest related works along these aspects below.

A. Logical recursion optimization
The ability to express recursive queries has been studied in various formalisms including RA, Datalog, and automata as an ad-hoc technique to evaluate regular queries like UCRPQs on graphs.
The system we present in this paper is based on the µ-RA algebra, an extension of RA with a fixpoint operator introduced in [11] which offers a good balance between expressivity, possible optimizations, and yields state-of -the-art performance for evaluating UCRPQs in a centralized setting with PostgreSQL [11].
Earlier works tried to extend RA with recursion, they are either less expressive than µ-RA or are more expressive but support less optimizations than µ-RA. See [11] and references thereof for more details. α-extended RA [9] introduces an operation α that expresses transitive closure, which means it has the same expressivity as UCRPQs. LFP-RA [10] introduces a least fixpoint operation and has the same expressivity as Datalog with stratified negation. Various fragments of LFP-RA have also been studied, notably the fragment that is restricted to linear recursion and is as expressive as linear Datalog. Finally, the WHILE [12] language is at least as expressive as LFP-RA. It is shown in [11] that µ-RA with the restrictions F cond (Sec. II-B) is as expressive as linear Datalog. So µ-RA is more expressive than α-extended RA, as expressive as the linear version of LFP-RA, and less expressive than the other formalisms.
In Datalog, Magic Sets [32], [33], recently improved by Demand Transformation [34], is a well-known optimization technique for recursive Datalog programs. The optimizations provided by Magic Sets or Demand Transformations are equivalent to pushing selections and projections in µ-RA. However, there is no equivalent to merging fixpoints. Furthermore, depending on the way the Datalog program is written, some optimizations may or may not be applied. For instance a left-linear DL program (e.g. P (x, y) ← P (x, z), R(z, y)) cannot push filters that are applied on the right side (on y in the example). The optimization framework has then to be coupled with a technique for reversing DL programs as proposed in [35]. Since Datalog engines use heuristics to combine optimization techniques, optimizations are not always performed as observed in [11]. In all cases, they lack the ability to merge fixpoints. In summary, Dist-µ-RA relies on µ-RA which enables optimizations at the logical stage (before distribution) that other systems are unable to achieve.

B. Optimizations in distributed data management frameworks
The seminal systems MapMeduce and Dryad are known to be inefficient for iterative applications [4].
Spark [3] and Flink [4] were introduced to improve upon these systems and became prevalent for large scale and dataparallel computations. Work in [36] proposes a technique that improves Spark task scheduling and thus performance for iterative applications. This system-level optimization is transparent for Spark applications and thus Dist-µ-RA can directly benefit from it.
Systems specifically designed for large-scale graph processing include Google's Pregel [5], Giraph [6] an open-source system based on the Pregel model, GraphLab [37] and Powergraph [38]. GraphX [7] is a Spark library for graph processing that offers a Pregel API to perform recursive computations. Pregel is based on the Bulk Synchronous Parallel model. A Pregel program is composed of supersteps. At each superstep, a vertex receives messages sent by other vertices at the previous iteration and processes them to update its state and send new messages. Computation stops when no new message is sent.
Is is not straightforward to evaluate UCRPQs in Pregel. An automata like algorithm needs to be written to know which stage of the regular query each processed path has reached. The idea is to traverse the paths in the graph (by sending messages from vertices to their neighbors) while traversing the regular query. [39] proposes a system that implements RPQ queries on GraphX and proposes optimizations to reduce communications between nodes. In all these systems, selections can be pushed in one direction only. For instance, if the program traverses the regular query from left to right, the execution of the program naturally computes the filters and edge selections occuring before a recursion first, thereby pushing these operations in the recursion. Selections which occur after the recursion cannot be pushed. Additionally, communications between workers of the cluster happen in every superstep of the recursion, which is avoided by the P plw plan in Dist-µ-RA. SociaLite [44] is an extension of Datalog for social network graph analysis. Its distributed implementation runs queries on a cluster of multi-core machines in which workers communicate using message passing. SociaLite does not offer a distribution plan equivalent to P plw where recursion can be executed without communication between workers at every step.
Myria [28] is a distributed system that supports a subset of Datalog extended with aggregation. Queries are translated into query plans that are executed on a parallel relational engine. Myria supports incremental evaluation of recursion and provides both a synchronous and an asynchronous mode. It does not support logical optimizations of the generated query plan involving the recursive operator like pushing joins in fixpoints and merging fixpoints. Myria does not either propose a distribution plan equivalent to P plw .
RaSQL [45] proposes an extension of SQL with some aggregate operations in recursion. Queries are compiled to Spark SQL to be distributed and executed on Spark. This system does not propose rules to push selections in the fixpoint operator nor to merge fixpoints. RaSQL also proposes a decomposable plan for recursion similar to the one in BigDatalog but has no automated technique to distribute data. The RASQL implementation is not available for testing.
BigDatalog [13] is a recursive Datalog engine that runs on Dist-µ-RA BigDatalog Graphx Figure 13. Running times on uniprot_1M.  Spark. It uses the Datalog GPS technique [46] that analyses Datalog rules to identify decomposable Datalog programs and determine how to distribute data and computations. These ideas are tied to Datalog and are not applicable to the relational algebra. The present work proposes a new method specifically designed for recursive relational algebraic terms. It uses the µ-RA filter pushing technique to automatically repartition data. Compared to BigDatalog, Dist-µ-RA is superior because it supports optimizations that BigDatalog is unable to provide. For instance, as mentioned earlier, BigDatalog cannot perform Dist-µ-RA's optimization that merges fixpoints. In fact, it is not possible to do so in the Datalog framework. This is known as an intrinsic limitation of Magic Sets and Demand Transformations at the core of Datalog optimizations. Another example is that BigDatalog cannot push all filters and joins that can be pushed, without requiring the support of techniques to reverse fixpoints. This is notably hard and BigDatalog does not implement such techniques. In comparison, Dist-µ-RA supports more logical optimizations of terms (and regardless of their initial form) and is thus capable of generating logical query plans that are beyond reach for BigDatalog. In practice, Dist-µ-RA provides superior performance on a wider range of query classes, as reported by the experiments in Section V. The only technique that Dist-µ-RA reuses from BigDatalog is a Spark optimized implementation of RDDs: SetRDDs which is a specialization of the Spark RDD for sets.

VII. CONCLUSION
The Dist-µ-RA system combines various techniques for the optimization and distribution of recursive relational queries in a single framework. Regarding the logical/algebraic aspect, it integrates well with the relational algebra. It inherits RA's advantages including the fact that queries are optimized regardless of their initial shape and translation in the algebra. With respect to distribution, different strategies for evaluating recursive algebraic terms in a distributed setting are provided. These strategies are implemented as plans with automated techniques for distributing data in order to reduce communication costs. Experimental results on both real and synthetic graphs show the effectiveness of the proposed approach compared to existing systems.