Hermano Lustosa
email: hermano@lncc.br

Fabio Porto
email: fporto@lncc.br

Pablo Blanco
email: pjblanco@lncc.br

Patrick Valduriez
email: patrick.valduriez@inria.fr

Database System Support of Simulation Data

Supported by increasingly efficient HPC infra-structure, numerical simulations are rapidly expanding to fields such as oil and gas, medicine and meteorology. As simulations become more precise and cover longer periods of time, they may produce files with terabytes of data that need to be efficiently analyzed. In this paper, we investigate techniques for managing such data using an array DBMS. We take advantage of multidimensional arrays that nicely models the dimensions and variables used in numerical simulations. However, a naive approach to map simulation data files may lead to sparse arrays, impacting query response time, in particular, when the simulation uses irregular meshes to model its physical domain. We propose efficient techniques to map coordinate values in numerical simulations to evenly distributed cells in array chunks with the use of equi-depth histograms and space-filling curves. We implemented our techniques in SciDB and, through experiments over real-world data, compared them with two other approaches: row-store and column-store DBMS. The results indicate that multidimensional arrays and column-stores are much faster than a traditional row-store system for queries over a larger amount of simulation data. They also help identifying the scenarios where array DBMSs are most efficient, and those where they are outperformed by column-stores.

INTRODUCTION

Scientific applications produce ever increasing amounts of data that need to be managed and analyzed. A numerical simulation (simulation for short) is a type of simulation, based on numerical methods, used to represent the evolution of a physical system quantitatively. It has been extensively adopted in many disciplines, such as: astronomy; medicine, biology, and oil and gas. It enables the investigation of phenomena where observed data is not available, such as the state of our universe after the big bang, or even when a complex and dangerous intervention, as in medicine, needs to follow some trials.

A simulation solves a set of differential equations, representing the physics involved in the studied phenomenon and its evolution in space and time. The simulation adopts a spatial representation of the phenomenon domain in the form of a geometrical mesh, composed of triangles or tetrahedrons, that guides the computation through a discretization of the space on its elements such as: vertexes, edges and faces. The solution of the equations finds the values for the predictive variables, such as velocity and pressure, and their variation along the discretized space and time. Moreover, a simulation may combine one dimension (1D) with 3 dimensional (3D) meshes. The latter requires more intensive computation and is left for areas of the simulation domain that demand more precise calculation.

Simulations need to have their initial conditions tuned in order to improve their quality. Aiming at finding the right parameter composition, modelers run parameter sweeps, that iterate the simulation execution with different parameter values. Thus, the combination of precise 3D space model and extended time-step simulation interval, with parameter sweeps leads to large simulation output files.

In this work, we are interested in managing the data produced by simulations. Such data can be viewed as a recurring point cloud that repeats itself along many time steps and simulation trials ( see Figure 1). Each point on a cloud is associated to predictive variables, whose values vary on each time step and simulation trial, forming a set of point clouds.

Typically, the simulation output is written to disk as raw files. The latter may be input into a visualization tool, such as Paraview, which allows the modeler to visually assess the simulation quality. Similarly, modelers may apply quantitative analysis on the predictive variable results by writing specific programs to process the simulation files. This method becomes inefficient and hard to scale as the volume of data increases. As an example, our colleagues at HeMoLab have written programs that load in memory large simulation datasets to capture the simulation results in a zone of interest, such as looking for the behavior on a brain aneurysm region [START_REF] Cebral | Efficient pipeline for image-based patient-specific analysis of cerebral aneurysm hemodynamics: technique and sensitivity[END_REF], which requires a query on a spacetime slice of the dataset. Traversing thousands of simulation datasets and comparing results in a particular region and time interval can be efficiently computed in a parallel DBMS, if the parameter sweep computation output is well stored in a distributed system.

Furthermore, due to the lack of isolation between the applications and the data, any eventual changes to the data (i.e. improving the simulation) might require the software to be changed as well. In this context, providing the benefits of a declarative query language for simulation results analysis and a system that scales to support large numerical simulation data will bring state-of-the art data management techniques to this area.

Figure 2 illustrates an analysis used to evaluate the simulation output. It comprises the calculation of the mean squared error between the simulation data (denoted by F ) and a reference function (denoted by G), in order to determine how much the resulting predictive values deviate from a pre-calculated reference value in every point. If the data is stored in a DBMS, this analysis can be expressed by a declarative (SQL-like) query, instead of requiring additional coding of another application.

Figure 2: Mean Squared Error Analysis Example

Traditional row-store DBMSs, usually tailored for OLTP workloads, are considered inadequate for the analytical work-load required by scientific applications [START_REF] Stonebraker | Requirements for science data bases and scidb[END_REF]. Stonebreaker has pointed [START_REF] Stonebraker | One size fits all: An idea whose time has come and gone[END_REF] the inability of current relational DBMS products to satisfy the requirements of a wide range of domains, and argues for specific database engines for different kinds of data centric applications with distinct data models.

Specific data models can offer not only a more convenient representation but may improve query performance through more adequate algorithms. An array DBMS, such as SciDB [START_REF] Brown | Overview of scidb: Large scale array storage, processing and analysis[END_REF], offers an interesting multidimensional model to represent simulation data, as it leads to a very intuitive mapping of the geometrical portion of simulation data to dimensions and cells in the array. Specifically, the space-time domain dimensions and the simulation trials are mapped to array dimensions, while predictive variable values are mapped to cells of the multidimensional array.

However, some data characteristics might pose difficulties for an array representation. Unstructured meshes, designed to represent complex physical domains, have an irregular distribution of points in space. Mapping such irregular distribution of points onto arrays can produce a sparse multidimensional cube, and incur into performance penalties during array transfer from disk into memory. As our experiments show, different mappings may produce variations in performance of more than 700%. For this reason, we investigate different techniques for managing the data produced by simulations, such as those coming from HeMoLab, by using an array DBMS.

Moreover, recent advances in in-memory column-stores for analytical workloads in scientific data [START_REF] Cijvat | Genome sequence analysis with monetdb: a case study on ebola virus diversity[END_REF][START_REF] Ivanova | Monetdb/sql meets skyserver: the challenges of a scientific database[END_REF], such as Mon-etDB [START_REF] Idreos | Monetdb: Two decades of research in column-oriented database architectures[END_REF], may also switch the balance from a choice based on data representation to an efficient relational DBMS. Therefore, we implemented our techniques in the SciDB array DBMS and compared them, through experiments over realworld data, with two other techniques based on PostgreSQL (regular row-store RDBMS) and MonetDB (column-store RDBMS). The experiments are carried out according to a benchmark composed by range queries and real life analysis used for improving the simulation.

The experimental results indicate that whenever the complete dataset fits in memory, MonetDB outperforms the other techniques. For large datasets, SciDB equipped with our techniques shows better results.

The remainder of this paper is organized as follows. Section 2 describes the cardiovascular system simulation, developed at LNCC, that motivated this work. Section 3 presents the multidimensional array data model. Section 4 describes our techniques to represent simulation data using an array DBMS. Section 5 describes the representation of simulation data using a relational DBMS. Section 6 presents our experimental evaluation. Section 7 discusses related work. Finally, Section 8 concludes.

NUMERICAL SIMULATION DATA

In order to run a simulation, one must first build the mathematical model that captures the behavior of the studied phenomenon. Commonly, mathematical models are expressed by differential equations that are solved by a variety of numerical methods, such as FDM (Finite Difference Method) and FEM (Finite Element Method). The software component that solves a set of differential equations applying a particular numerical method is called a solver. The adopted numerical method might require the discretization of the physical domain in a form of a mesh.

Meshes have topological and geometrical representations. Geometrical aspects are related to shapes, sizes and absolute positions of their elements, such as vertexes or points. The topology representation captures the relationships among elements, like their neighborhood or adjacency, without considering their position in time and space.

The distinction between these two concepts, i.e. mesh geometry and topology, is important in the context of our work. We believe that, although the topological representation might require a more elaborated solution [START_REF] Howe | Gridfields: Model-driven Data Transformation in the Physical Sciences[END_REF], both arrays and relations can be sufficient to represent the geometry of a mesh. Indeed, important queries can evaluate the simulation output by assessing the values of predictive variables associated to mesh points in time, irrespectively of topological constraints. Thus, in this paper, we constrain our discussions to geometric aspects of meshes.

Vertexes in the geometrical representation of a mesh can be viewed as forming a point cloud. A mesh guides solvers during simulation computation, such that, the predictive variables computed by the model have their values assigned at each vertex point, for each time step. Moreover, runs with different parameters of the same model may use the same geometry mesh as space discretization. The resulting output data, shown in Figure 1, can be viewed as a recurring point cloud appearing at different time steps and simulations.

Use Case

The HeMoLab group at LNCC has developed computational tools to simulate the human cardiovascular system [START_REF] Hemolab | Hemodynamics modeling laboratory[END_REF]. The simulation adopts a geometrical representation of the arteries in the human body in the form of a mesh. This representation is used in 1D and 3D simulations, working together as a coupled model [START_REF] Blanco | On the potentialities of 3d-1d coupled models in hemodynamics simulations[END_REF].

One dimensional (1D) models represent the human cardiovascular system in a simplified form. Main arteries of the human body are modeled as a set of lines. Three dimensional (3D) models are used in order to achieve a higher level of detail about the behavior of the blood flow inside an artery.

Figure 3 shows how the 1D and 3D models are arranged in space. Both 1D and 3D meshes actually contain points in 3D with x, y and z coordinates. However, the 1D model contains a linear representation of the arteries in the entire body, while the 3D model has a much more detailed representation of a single artery. Combined, 1D and 3D meshes have over 90 thousand points, most of them (about 99%) belong to the 3D model and the remaining 1% to the 1D model.

Figure 3: 3D-1D coupled model

The numerical solver, which is in charge of executing the simulation, computes for each position in the mesh, and for many time steps, physical quantities; such as: pressure, blood flow velocity and mesh displacements in the case of dealing with deformable domains. The solver reads a set of input files that describe the parameters of the model and stores the results in raw files, which requires scientists to write specific applications for any kind of analysis needed. In this context, storing simulation data in a DBMS would allow researchers to use a declarative query language to execute analyses more efficiently.

These analyses could be as simple as retrieving data from a subset of simulations or time steps, or even obtaining the result from an aggregation function for a spatial window (range query). In our benchmark we include more complex real life analyses, such as calculating the mean squared error between the data in the simulation and a reference function (see Figure 2), and obtaining the time step required for an attribute to reach its maximum value on each point.

The Human Cardiovascular simulation and the analytical queries mentioned above will form the basis of our benchmark to evaluate strategies for managing point cloud data.

MULTIDIMENSIONAL ARRAYS

In this section, we describe the multidimensional array data model, and its implementation in SciDB. We chose to use SciDB as the array DBMS in our evaluation because it directly implements the model instead of just adding array processing capabilities to a relational system. Also, SciDB is a distributed system, tailored to the requirements of some scientific application domains.

Multidimensional arrays are defined by a set of named dimensions. A set of indexes for all dimensions identifies a cell in the array. Cells can have many attributes in which the values are stored just like tuples in a relational DBMS. Cells can be non-existing or empty, and arrays are said to be sparse if they have many empty cells.

SciDB partitions the array into regular tiles, or chunks [START_REF] Brown | Overview of scidb: Large scale array storage, processing and analysis[END_REF]. Chunks are physical units of distribution of arrays on computer nodes. The array is partitioned by ranges of indexes of fixed size on each dimension. The product of the number of partitions on each dimension gives the number of chunks the array will be split into.

Figure 4 shows a bi-dimensional array with dimensions i and j. The chunk size for dimension i is 4 while its total length is 16. The chunk size for dimension j is 5, and its total length is 20. This configuration results in a total of 16 chunks. SciDB allocates chunks to worker nodes that are responsible for maintaining portions of the arrays.

We believe that arrays are a convenient model for representing simulation data. Different experiments, time steps and the spatial coordinate system can be mapped to dimensions in the array, while the physical quantities can be stored as attributes. The DBMS takes advantage of the array structure in order to answer range queries for dimension values. Dimensions work like multivalued indexes in a relational DBMS, which means there is a fast access path to a set of cells when they are specified by index ranges. In practice, portions of the array for each dimension or combination of dimensions can be retrieved efficiently by the DBMS.

Regular tiling or chunk partitioning also plays an important role in improving performance. Storing subarrays on disk as chunks yields an efficient representation of data, i.e., cells that are close to each other on the array, will likely [START_REF] Brown | Overview of scidb: Large scale array storage, processing and analysis[END_REF]. The advantage of this representation is that when a subarray is referenced by a query, the data will probably be sitting on a limited amount of disk pages where the chunks containing the data of interest are stored, instead of being scattered across many different pages. Obviously, a reduced amount of disk pages read from disk results in a faster query processing time.

ARRAY REPRESENTATION

There is an intuitive simplicity regarding the mapping of simulation data to multidimensional arrays. However, a naive mapping of coordinate values onto indexes of an array may lead to inefficient access due to the following data characteristics:

• the subjacent geometry mesh coordinate system may adopt a non-integer representation.

• the physical geometry of a phenomenon may be modeled by an irregular mesh.

Considering the mapping of dimensions involved in a simulation, the time-step and the simulation trial dimensions are both integers and do not include holes in their series, which leads to a more straightforward mapping to indexes in the array. Conversely, a mesh whose points are referenced by non-integer coordinates, when naively mapped onto integer array indexes, produces huge gaps between neighboring cells within the array. In this scenario, a space preserving mapping strategy would induce huge variable density chunks due to an uneven distribution of points.

Uneven distributed data can directly impact on the cells distribution per chunk. As the latter is used as I/O access unit in SciDB, the performance during query processing is affected by this irregular data distribution. In the following sections we describe three different techniques for mapping the coordinate values to integer indexes. We start by formalizing the terminology we will use.

• Definition 1. A geometry mesh is defined as M = (P, E), such that P = {p1, p2, • • • , pn} is a set of vertexes and E ⊆ P × P . Moreover, each vertex pi = (id, Xn), in which Xn is a n dimensional coordinate index and id is an identifier.

• Definition 2. An array A is defined as

A = (D, C), such as D =< D1, D2, • • • , Dm > is a list of m di- mensions and C = {c1, c2, • • • , cv} is a set of v at- tributes. Each dimension Di is indexed from 1 to an integer maxindex, in one increments. An index i =< i1, i2, • • • , i k >, with k ≤ m, and 1 ≤ ij ≤ maxindexj, for 1 ≤ j ≤ k , defines a subarray of A.
To deal with the problems caused by the irregularity of data distribution, we classify the mapping techniques into two groups: space preserving and non-space preserving.

Precision Elimination Method

We denote by precision elimination, (P.E. for short), the naive method for mapping coordinate values with finite precision to integer values. It consists of multiplying each different coordinate value by 10 P r , where P r is the precision in which the coordinate values are specified.

Suppose a mesh with dimensions [X1, X2, ... , Xn] containing the point p(x1, x2, ... , xn) specified with P r decimal places of precision. The array used in this case also contains the dimensions [D1, D2, ... , Dn], and the data related to the point p(x1, x2, ... , xn) will be stored on the cell identified by the indexes (x 1 , x 2 , ... , x n ) defined as

x i = xi * 10 P r .
The method is efficient and enables spatial window queries to be expressed easily. Consider a region of interest R in the original mesh representation, specified by their bordering points: lower bounds (lr1, lr2, ... , lrn) and upper bounds (ur1, ur2, ..., urn), where lri/uri is the lower/upper bounds for all i dimension, 1 ≤ i ≤ n. The subarray S containing only the cells related to the points within R can be defined with the lower bounds lsi = lri * 10 P r and the upper bounds usi = uri * 10 P r , for all dimensions i, 1 ≤ i ≤ n.

However, this method has some drawbacks. Multiplying the coordinate values in order to eliminate the precision may be undesirable in some applications. The index values will need to be truncated in case the coordinate values are specified with many decimal places. It might lead to information loss if the original mesh representation is not explicitly stored somewhere else.

In addition, multiplying coordinate values is equivalent to scaling up the points. The higher the precision, the higher the scale factor is, and the farther apart adjacent cells will be in the array. The highly sparse array obtained can be hard to deal with, as the cells are distant from each other, and unevenly distributed throughout array regions. Defining the partitioning in this case is challenging, and some chunks will likely be excessively filled while other chunks will contain just a few cells.

As we show in our experimental evaluation, unbalanced sparse chunks can negatively affect performance. For this reason, we propose the following two methods in order to avoid the sparse and unbalanced data distribution.

Non Space Preserving Methods

The naive mapping described above is a space preserving method, i.e., the general form of the mesh is preserved in the final array representation. However, in order to cope with the irregular data distribution one must rearrange the points altering the form of the mesh. We denote by non space preserving methods, the techniques in which the alteration of coordinate values does not keep the original spatial distribution of mesh points.

The first step for such methods is to create a more compact data representation. Dead spaces between coordinate values are removed without changing the relative positions of the points to each other. This process reduces the spacing between adjacent points and creates a less sparse array.

A sorted list Li is created for each X1, X2, ... , Xn dimensions, with 1 ≤ i ≤ n. Each list Li contains the coordinate values ai0, ai1, ai2, ... , aim that appear in some point on the mesh for the dimension i. Every referenced coordinate value aij has its position in the compacted representation given by a function Li(aij). A point p(x1, x2, ..., xn) on the mesh is replaced by p(l1(x1), l2(x2), ..., ln(xn)). This means that the existing space among points in the mesh have been removed, thus producing a more compact array representation.

Figure 5.2 shows a resulting compact representation. Note that although the shape of the mesh is deformed, because the spacing between adjacent coordinate values has been removed, the relative position of the points is maintained. This property is important, since the compact form must allow spatial window querying on the original representation to be translated into the new representation. Given a window R enclosing a set of points in the mesh, it must be possible to define a window R on the compact representation that encloses the exact same set of related points.

This statement can be verified as a coordinate value that appears on many points will be consistently replaced by the same integer value wherever it appears. Not only this, but the integer values are attributed in a fashion that maintains the ordering of the coordinate values. As all points are changed in the same manner, they maintain their relative position to each other, making it possible to translate spatial window queries.

Equi-depth Histogram Method

The equi-depth histogram based partitioning method is a non space preserving mapping strategy. We adopt it in oder to minimize the irregularity in data distribution and generate more balanced chunks. Equi-depth histograms are similar to equi-width histograms. The main difference is that regular histograms have bins for equally sized portions of the data domain, but containing a different number of occurrences on each one, while equi-depth histograms have bins for portions of the domain with variable sizes, but containing roughly the same amount of occurrences on each one. Thus, equi-depth histograms help splitting the mesh into regions containing a close number of points on each one. These regions will be used to define the chunk partitioning scheme of the multidimensional array.

Figure 5 depicts the steps executed to obtain the mapping in a 2D projection for the HeMoLab's 3D model mesh. In Figure 5.1, we have the original points in the mesh projected in a 2D surface. Figure 5.2 shows the resulting representation after the first step comprising the removal of dead spaces.

A set of equi-depth histograms, one for each dimension, is created in the following step (Figure 5.3). The combination of histograms for each dimension creates variable sized regions, which enclose a variable amount of points. The difference in the amount of points contained in the most populated and the least populated regions tends to decrease as we increase the granularity of the partitions (decrease the bin sizes for the histograms). Thus, we can vary the sizes of The regions created by the combination of the histograms will be used as the base for the definition of the chunk configuration. Ideally, chunks in SciDB should occupy a few megabytes on disk [START_REF] Brown | Overview of scidb: Large scale array storage, processing and analysis[END_REF]. In other words, chunks should contain between 500 thousand to 1 million cells for arrays with 8-bit double precision floating point attributes. Therefore, the amount of points on most regions should be large enough so the chunks have enough cells. The higher the granularity, the smaller is the amount of points on each region. Thus, we search for the more granular partitioning that still enables to define sufficiently populated chunks.

As stated previously, the regions depicted in Figure 5.3 have variable sizes, so they cannot be directly used as the chunk configuration for an array in SciDB, due to the fact that it only allows regular tiling. To overcome this limitation, points are translated so they become distributed in equally sized regions of space.

Figure 5.4 depicts the mesh after the points have been translated. The process to define the new position for a point is as follows. Every histogram Hi has a set of bins Bi0, Bi1, ... , Bi1 for every dimension i. A coordinate value xi on a point p of the mesh is contained in the range of values associated to bin Bij (the j th bin for the histogram associated with the dimension i). The index j of a bin Bij associated to a coordinate value xi is denoted by Bi(xi). We define as IBi(xi) the initial coordinate of a bin, i.e., the smallest value contained in a bin Bi that also contains the coordinate value xi.

The size of the new equally sized regions in the dimension i will be given by the largest bins LBi for the histograms Hi. Note that in Figure 5.4, the new regions of the space have the length in dimension i of the largest bin for the histogram Hi. Every point p(x0, x1, ..., xn) is replaced by a new point p (x 0 , x 1 , ..., x n ). The new coordinate values for p are given by the formula:

x i = LBi * Bi(xi) + (xi -IBi(xi)) for 1 ≤ i ≤ n.
The array used to store data mapped with the use of this method has the dimensions [X1, X2, ... , Xn] and the coor-dinate values of every p will be used as indexes to allocate the points in the array. The LBi values are used to define the chunk sizes for each dimension. It is important to note that even after moving the points, we still are able to translate spatial window queries, as long as the mapping between the original coordinate values and the indexes used for the points is explicitly stored.

The final array representation is still sparse, but the cells will be distributed much more evenly among the chunks, which is the main advantage of this technique. This method not only maps the coordinate values into indexes, but also helps defining the chunk configuration, which is an important aspect in designing multidimensional arrays in SciDB.

However, the process for creating this mapping might be costly as we potentially need to create many histograms and evaluate the partitions given by their combinations. Additionally, every spatial window query specified based on the original mesh representation must be translated into a region in the array.

This translation can be described as follows. Suppose a region of interest R in the original mesh representation, specified as lower bounds (lr1, lr2, ... , lrn) and upper bounds (ur1, ur2, ..., urn). We define a function Mi(xoi) → x i that maps the original real coordinate value to the corresponding index stored in the array. The values for this function are pre-calculated and explicitly stored, thus, evaluating the function requires a lookup in the structure used to store the mapping. The subarray S containing only the cells related to the points within R can be defined with the lower bounds lsi = Mi(lri) and the upper bounds usi = Mi(uri) for every dimension i.

We consider the use of this technique under two main assumptions about the underlying mesh:

• It has low dimensionality (usually 2D or 3D).

• It is a relatively small portion of the entire dataset.

In our use case, the main dataset contains data from 20 simulations, each one possessing 240 time steps. The number of data records is equal to the product of the number of points in the mesh (close to 90 thousand), the number of simulations and time steps, which amounts to over 430 million. Since the number of points in the mesh is much smaller than the amount of records in the entire dataset, maintaining and querying the mapping, that is as large as the amount of points in the mesh, should be much cheaper than maintaining and querying the entire dataset.

Space-Filling Curve Method

Although a more balanced data distribution is achieved with the use of the second method, the resulting array is still sparse and contains irregularly filled chunks. Our third method to map coordinate values to array indexes relies on the use of space-filling curves to create a dense data representation.

Space-filling curves can be understood as a path that defines a order in which points are visited. They enable the mapping of multidimensional points into a single scalar value or code, that depends on the position of the point in the curve. Space-filling curves are defined in such manner that adjacent codes will be attributed to points that are relatively close to each other in space, because the next point the curve passes by is likely to be on the neighborhood of the last visited point. This characteristic of space-filling curves makes them useful in indexing points, because it helps maintaining a coherent data representation. With this method, the original mesh dimensions X1, X2, ... , Xn are represented by a single dimension D in the array. A single index is attributed for points according to a space-filling curve, in such a manner that adjacent indexes are given to points that are close in space.

Spatial window queries, previously defined as a set of upper and lower bounds for different dimensions, now have to be expressed with multiples ranges of values for a single dimension. The points enclosed in a window R of the original mesh representation will be distributed in a set of continuous intervals of code values [L1, U1], [L2, U2], ..., [Ln, Un]. The best, and very unlikely, scenario is when all points are contained in a single interval. The worst scenario is when the number of intervals is equal to the number of points within R, which means that the codes are scattered in noncontiguous intervals.

In a multidimensional array, fewer intervals for a range query means data will probably be sitting on a limited number of chunks. Fewer chunks, means faster retrieval time, and faster query processing time. It is common that just a portion of the data contained in the chunk is necessary to answer a query. Therefore, minimizing the amount of chunks read from disk to memory is also going to minimize the amount of unnecessary data loaded along the data of interest.

The Hilbert Space-Filling curve has very interesting clustering properties [START_REF] Moon | Analysis of the clustering properties of the hilbert space-filling curve[END_REF] that help creating an efficient data representation in a single dimension. We use compact Hilbert indexes [START_REF] Hamilton | Compact hilbert indices for multi-dimensional data[END_REF] to allocate the positions of the points in the array. Figure 6 shows a mesh with vertexes colored by their respective compact Hilbert indexes. Note that large clusters of points in a same region of the mesh have similar colors, indicating their respective indexes have also close values. The process of attributing indexes to points goes as follows. First, we create a compact integer representation of the mesh as explained in section 4.2. Then, for every point p(x0, x1, ..., xn) in the compact representation, we obtain chi(p), that is the compact Hilbert index for p. All points p are ordered according to their chi(p) value and are added to a list L chi . Finally, the index for the point in the array is attributed according to the order of the points in L chi . This is done because not necessarily the chi(p) values will be contiguous. There may be intervals of indexes that are not attributed to any point p. These intervals would create an undesirable sparse distribution, so they are removed once the indexes are attributed according to the order in which the points appear in L chi .

In order to enable querying and rebuilding the original space distribution, an auxiliary one-dimensional array must be created, containing a dimension D and an attribute set [Xa0, Xa1, ..., Xan] for every dimension in the original mesh representation where each attribute holds values for a given mesh dimension. This array contains the mapping between the original coordinate values representing a point and the associated indexes. Spatial window queries are expressed with joins between the main array and the auxiliary array, which is slightly less efficient than directly expressing the boundaries of the spatial window.

Even though range queries might not be as efficient in this representation, the overall performance of queries that represent more elaborated data analyses are expected to be improved for two main reasons. The resulting array is dense, and the data distribution per chunk is perfectly balanced. The array has a smaller dimensionality that might affect positively the performance for executing array operations that are the building blocks of complex queries.

RELATIONAL REPRESENTATION

This section briefly describes the relational structure used to represent simulation data in relational DBMSs. In our evaluation, we use both a row-store DBMS and a columnstore DBMS, respectively PostgreSQL and MonetDB.

Schemas in both systems are roughly the same. The data is stored in a single table containing attributes to represent the physical quantities, the mesh coordinate system, the time and the simulation identifiers. A single table offers a simple representation of data and does not require costly joins in order to answer range queries.

Different from the array representation that requires a mapping between coordinate values to integer indexes, the relational representation naturally allows the original mesh coordinate values to be directly stored as double precision floating point numbers in columns of the relation. Any range query can be expressed as a logical predicate using the SQL operator BETWEEN for defining intervals.

The primary index for the relation is composed by the set of attributes containing the combination of the coordinate values for the points in the mesh and the time and simulation identifiers. This combination of values uniquely determines another set of attributes that represent the physical quantities, such as pressure, velocity and displacement related to the points in a instant in time.

Our preliminary tests showed that independent secondary indexes have better performance for the workload under consideration in our benchmark. Thus, we use secondary indexes for every attribute on the primary key, enabling the DBMS to decide which index or combination of indexes will be more efficient. The definition of secondary indexes is necessary only for PostgreSQL, since MonetDB implements the Imprints index structure [START_REF] Sidirourgos | Column imprints: A secondary index structure[END_REF] that is created during query processing.

EXPERIMENTAL EVALUATION

Environment

The benchmark was executed in a virtual machine environment (single core, 12GB of memory and 250 GB of storage), with Ubuntu Server 14.04 LTS as both the guest and the host system. The choice of a single core VM aimed at reducing uncontrolled internal system parallelism that could affect results. We used versions 14.12, 5.0 and 9.4 for SciDB, MonetDB and PostgreSQL respectively. For the comparative tests, all DBMSs were running with a single instance, while for the distributed tests we evaluate SciDB alone using 4 and 8 instances distributed respectively in 2 and 4 nodes.

The Benchmark

The benchmark is designed considering typical quantitative analysis performed while evaluating the quality of simulation or supporting the interpretation of the modeled phenomenon. It includes the data from the cardiovascular system simulation implemented into SciDB, MonetDB and PostgreSQL and four set of analytical queries.

Data Definition

We defined the schemas in SciDB using the DDL for the AQL (Array Query Language). AQL is one of the two languages supported in SciDB. The following CREATE AR-RAY commands define the schemas for arrays that implement the P.E., Histograms Based Partitioning and Hilbert Space-filling Curve methods respectively (main and auxiliary mapping arrays). Attributes are listed between angle brackets and dimensions are specified between square brackets. Each dimension is defined with a name and boundaries, followed by two numbers representing the chunk size and the overlap between chunks. Chunks sizes are defined according to each strategy. Larger chunks are specified for the very sparse PrecisionElimination array, while more compact chunks are defined by the Histogram Based Partitioning strategy for the Histogram array. The Hilbert array contains dense chunks whose sizes on each dimension are defined in such manner that the total amount of cells is 1 million (5000 * 50 * 4).

CREATE ARRAY P r e c i s i o n E l i m i n a t i o n

Tables for the relational schema are defined with the following CREATE TABLE command. The only difference between the schemas is the addition of secondary indexes in PostgreSQL.

CREATE TABLE [ p o s t g r e S Q L T a b l e / monetDBTable ] ( s i m u l a t i o n

BIGINT , t i m e s t e p BIGINT , x DOUBLE, y DOUBLE, z DOUBLE, p r e s s u r e DOUBLE, PRIMARY KEY ( s i m u l a t i o n , t i m e s t e p , x , y , z ) ) ;

The dataset consists of 20 simulations comprising 240 time steps and over 90 thousand points for the 3D-1D mesh. The total amount of records (rows and cells) is bigger than 430 million.

Queries

Four sets of queries are used in our benchmark. Two sets are composed of spatial window queries (range queries) evaluating the results of an aggregation function on a region of space. The regions of space selected are the complete 1D model, and a portion of the 3D model representing an aneurysm formed on the artery.

The 1D model region of the mesh is very sparse and has a relatively small amount of points compared to the entire mesh. The aneurysm region contains much more points, that are arranged in a denser region of space. These two opposite regions of the mesh are chosen to evaluate the performance of each strategy in different scenarios.

We denote by Hemolab 1 the third set of queries in our benchmark. Hemolab 1 comprises the calculation of the mean squared error between the data in the simulation and a reference function. The values for this reference function are precalculated and stored in the database. The final result for this query is a single value for each point for every simulation under consideration. The value is the mean squared error between the values on the simulation and the reference function on each point of the mesh. This query evaluates the execution of a join operator based only on dimension values.

The hemolab 1 query in AFL and SQL is: The fourth set of queries is denoted by Hemolab 2. This analysis is executed in order to obtain the time required for attributes in the simulation to reach their maximum values for the first time. Maximum values of the attributes on every point are calculated with use of an aggregation function.

After that, values are joined with the original dataset so they can be paired up with the time steps in which they occur. Once we have all time steps in which the attribute value is at its peak, we simply return the lowest time step, i. e., the first to present the maximum value. The result for this query is also a single value for every point, that represents the time necessary for the attribute to reach its maximum value. This query contains a join based on attribute values, and it evaluates performance for queries that need to perform this kind of operation.

The hemolab 2 query in AFL and SQL is:

Experiments with SciDB

This section presents the experimental results for different methods of representing simulation data in SciDB arrays. The charts contain the average execution time for 30 runs. Queries were executed with three configurations, containing 1, 4 and 8 SciDB instances. Tags below the horizontal axis of the charts in this section indicate the number of instances related to the execution time. All queries in this section evaluate data from all time steps and simulations. This is true even for spatial range queries.

Figure 7 shows the results for the 1D model range query. The 1D model mesh contains a relatively small portion of the points. In this case, the sparse representation generated by the P.E. method was more efficient. The irregular data distribution in this case creates a set of almost empty chunks containing data for the points on the 1D model only. These chunks do not contain the dense data for the 3D model, and thus, the DBMS does not read and process unnecessary data.

The exact opposite happens for the Histogram Base Partitioning method. In this case, the more balanced chunks contain data from both 3D and 1D models, which means SciDB needs to read chunks filled with data unnecessary to answer the query, making the total execution time many times higher. Figure 8 shows the results for the range query evaluating the region encompassing the aneurysm in artery. This is a subset of the data from the 3D model only. The results differ from those of Figure 7. In general, there is no chunk configuration that is good for all range queries. Ideally, a dense and well distributed configuration, such as the one provided by the Hilbert space-filling curve method, would minimize the variations between distinct queries, improving the overall processing time for different queries.

The irregular precision elimination method is two times slower than other strategies for the aneurysm range query. Since the entire 3D model is distributed in a small portion of space, the related points are also sitting on a limited amount of chunks. In order to retrieve data for just a portion of the 3D model enclosing the aneurysm, the system needs to read chunks filled with cells for the entire 3D model. The Histogram Based Partitioning method has the best performance in this case, because the granular and balanced distribution of 3D model points enabled SciDB to access chunks that encompass the data of interest more closely, thus diminishing the need to read from disk regions of the array that do not intersect with the range query. Figure 11 shows the time required to load data into the DBMSs in each case. A COPY command was used to load the entire dataset at once for relational systems, and specific array operators, such as LOAD and REDIMENSION, were used to load the data into SciDB.

MonetDB is the most efficient because all it does is to append data to files for each column on the relation. Both PostgreSQL and SciDB need to execute costly operations during data loading. Operations responsible for ordering data and dividing it into chunks take a longer time than just adding records to data files. Also, PostgreSQL needs to maintain the secondary indexes for each attribute on the primary key, which slows down the row insertions. PostgreSQL shows better results than most other methods, except P.E., as the amount of data under consideration is very small, and the system is able to use a combination of indexes to retrieve a minimum amount of disk pages to answer the query. However, as the next results show, Post-greSQL is not able to maintain fast execution time when we increase the amount of data.

Figure 13 shows the result obtained for the range query evaluating a region of the 3D model enclosing an aneurysm. For this query, and all following analyses, PostgreSQL has the worst performance by a great margin. Both MonetDB and SciDB are tailored for analytical workloads, and offer MonetDB is very efficient in the second analysis. Even though no index is previously defined, it is able to execute the analyses almost as rapidly as SciDB with the Histogram Based Partitioning Method. Another point to highlight is that the Hilbert Space-filling curve method obtained intermediate results in both range queries. The need to execute a join operation slightly increases execution time, but the dense array representation with lower dimensionality makes it still a good alternative, even for range queries.

Figures 14 and15 show the results for the queries Hemolab 1 and Hemolab 2 respectively. Each figure contains two charts, the upper chart gives the results for all systems and the bottom chart for MonetDB and SciDB.

The hemolab 1 results for analyses with 1 and 10 simulations are favorable to MonetDB. However, the Hilbert Space-filling curve yields a better performance when the query is scaled up to 20 simulations.

MonetDB has better performance for Hemolab 2 in all cases. Hemolab 2 in SciDB requires a join to be performed between attributes of the array instead of dimension indexes, which is inefficient. MonetDB column-store model, along with its dynamic index creation, is able to completely outperform SciDB in this scenario.

The Hilbert Space-filling curve method offers the best option among the multidimensional array representations. It is 60% faster for Hemolab 2 than P.E method and the Histograms Base Partitioning method. The dense and balanced data distribution allied with the lower dimensionality enables analyses to be executed more efficiently.

Comment on Execution Profile

MonetDB execution model heavily relies on the materialization of intermediate query results. When these intermediate results do not fit into main memory, they must be swapped to disk, resulting in performance loss. SciDB, on the other hand, has an execution model based on pipelining, and intermediate arrays do not always need to be completed materialized during query execution. Therefore, SciDB performance was less dependent on the amount of memory present on the system.

The differences between SciDB and MonetDB execution models is further evidenced in Figure 16

RELATED WORK

In this work, we are interested in managing scientific data generated by numerical simulations computed according to a space discretization modelled as an irregular mesh. There are few relevant works that discuss this topic.

A notorious example of work in mesh data management is Howe's [START_REF] Howe | Gridfields: Model-driven Data Transformation in the Physical Sciences[END_REF] Gridfields data model, developed to represent computational meshes. The model is not bounded to any data structure, and separates mesh topology and mesh geometry. Another model designed to represent unstructured meshes is ImG-Complex [START_REF] Rezaei Mahdiraji | Img-complex: graph data model for topology of unstructured meshes[END_REF]. This model represents meshes by using incidence multi-graphs. Both works deal with the topological representation of meshes and its underlying data. In our use case, we are concerned about executing analyses depending only on the geometric aspects, and thus, we do not require a specific model for dealing with topology yet.

Histograms have long been used in database management, mainly in query optimization. A technique to generate histograms for multidimensional data taking into account its sparsity is proposed in [START_REF] Furfaro | Compressed hierarchical binary histograms for summarizing multi-dimensional data[END_REF]. The use of histograms in query optimization and load balancing for cloud environments is discussed in [START_REF] Shi | Hedc: A histogram estimator for data in the cloud[END_REF]. These works are only examples of how histograms have been widely used in data management. The main point to highlight is that instead of using histograms in query optimization, we use equi-depth histograms as the basis for the multidimensional array partitioning strategy over distributed nodes so that each chunk holds approximately the same number of points, obeying the ordering on each dimension.

There is an extensive literature [14][11] in multidimensional indexing structures for databases using space-filling curves. In [START_REF] Martinez-Rubi | Benchmarking and improving point cloud data management in monetdb[END_REF], there is a discussion involving the use of space-filling curves in point cloud data management with MonetDB. The main difference with our work is that instead of having a dataset of simple point cloud data, we have a recurring point cloud structure associated with other dimensions, and we evaluate the use of multidimensional arrays and a space-filling curve to represent it.

CONCLUSION

In this paper, we proposed the use of array DBMSs to manage the data produced by simulations. We considered multidimensional arrays as it nicely models the dimensions and variables used in numerical simulations. We presented methods for the efficient mapping of variable values in simulations to evenly distributed cells in array chunks with the use of equi-depth histograms and space-filling curves.

We also considered the use of a row-store DBMS, Post-greSQL, and MonetDB a column-store DBMSs designed to provide better performance for analytical workloads. Unlike the multidimensional array data model, the relational model does not require any modification in the original data.

Our experiments show that using our techniques in an array DBMS produces a consistent allocation of simulation data within chunks with respect to a query workload. Each technique shows better performance results in comparison to MonetDB and PostgreSQL, at least on a particular class of the benchmark queries. Furthermore, our solution becomes less sensitive to restrictions on available memory than MonetDB. Nevertheless, both SciDB and MonetDB showed significant better performance than PostgreSQL, indicating that row-store DBMSs are not a good solution for analytical queries. Conversely, MonetDB offers a much simpler data representation with generally good performance.

We observe that the choice of DBMSs, as expected, impact on our results. For instance, column-stores adopting an in-memeory pipeline execution model may not incur in the problems observed with MonetDB for large datasets. Similarly, the mapping stucture needed by the Histogram-based method could have been implemented in a relational system with spatial indexing support, improving the query resolution performance.

Figure 1 :

 1 Figure 1: Recurring Point Cloud Data on Different Time Steps and Simulations

Figure 4 :

 4 Figure 4: Chunk configuration

Figure 5 :

 5 Figure 5: Equi-depth Histogram Based Partitioning Method

Figure 6 :

 6 Figure 6: 3D Model Mesh Points Colored by Their Hilbert Codes

  <p r e s s u r e : d o u b l e > [ x = 0 : * , 3 5 0 0 0 0 , 0 , y = 0 : * , 3 5 0 0 0 0 , 0 , z = 0 : * , 3 5 0 0 0 0 , 0 , t i m e s t e p = 0 : * , 5 0 , 0 , s i m u l a t i o n n u m b e r = 1 : * , 4 , 0 ] ; CREATE ARRAY H i s t o g r a m <p r e s s u r e : d o u b l e > [ x = 0 : * , 1 6 0 9 1 , 0 , y = 0 : * , 2 4 5 6 5 , 0 , z = 0 : * , 2 5 1 8 6 , 0 , t i m e s t e p = 0 : * , 5 0 , 0 , s i m u l a t i o n n u m b e r =1: * 4 , 0 ] ; CREATE ARRAY HistogramAux <x : d o u b l e , y : d o u b l e , z : d o u b l e , x map : d o u b l e , y map : d o u b l e , z map : d o u b l e > [ i : * , 1 0 0 0 0 0 0 , 0 ] ; CREATE ARRAY H i l b e r t <p r e s s u r e : d o u b l e > [ h i l b e r t C o d e = 0 : * , 5 0 0 0 , 0 , t i m e s t e p = 0 : * , 5 0 , 0 , s i m u l a t i o n n u m b e r = 1 : * , 4 , 0 ] ; CREATE ARRAY H i l b e r t A u x <x : d o u b l e , y : d o u b l e , z : d o u b l e > [ h i l b e r t C o d e = 0 : * , 5 0 0 0 , 0 ] ;

  ( [ARRAY] , p r e s s u r e ) , n u l l , n u l l , [ SIM NUM ] , n u l l , n u l l , [ SIM NUM ] ) a s a , [ REFERENCE ARRAY ] a s b , a . t i m e s t e p , b . t i m e s t e p , a . s p a c e , b . s p a c e ) , prod , ( a . p r e s s u r eb . p r e s s u r e ) * ( a . p r e s s u r eb . p r e s s u r e ) ) , sum ( p r o d ) , s i m u l a t i o n n u m b e r , s p a c e ) SELECT SUM( ( a . p r e s s u r eb . p r e s s u r e ) * ( a . p r e s s u r eb . p r e s s u r e ) ) FROM (SELECT p r e s s u r e , s i m u l a t i o n , x , y , z , t i m e s t e p FROM [ TABLE ] WHERE s i m u l a t i o n = [ SIM NUM ] ) AS a INNER JOIN (SELECT p r e s s u r e , s i m u l a t i o n , x , y , z , t i m e s t e p FROM [ REFERENCE TABLE ] ) AS b ON a . x = b . x AND a . y = b . y AND a . z = b . z AND a . t i m e s t e p = b . t i m e s t e p GROUP BY a . x , a . y , a . z , a . s i m u l a t i o n ; "

  APPLY(BETWEEN(PROJECT ( [ARRAY] , p r e s s u r e ) , n u l l , n u l l , 1 , n u l l , n u l l , 1 ) , t i m e , t i m e s t e p ) AS a , AGGREGATE(BETWEEN( [ARRAY] , n u l l , n u l l , 1 , n u l l , n u l l , 1 ) , max ( p r e s s u r e ) AS max , s i m u l a t i o n n u m b e r , s p a c e ) AS b , a . s i m u l a t i o n n u m b e r , b . s i m u l a t i o n n u m b e r , a . s p a c e , b . s p a c e ) , a . p r e s s u r e = a . max ) , min ( a . t i m e ) , s i m u l a t i o n n u m b e r , s p a c e ) AS a , PROJECT(APPLY ( [ ARRAY MESH ] , norm , SQRT( x * x + y * y + z * z ) ) , norm ) AS b , a . s p a c e , b . s p a c e ) SELECT A . s i m u l a t i o n , a . x , a . y , a . z , min ( a . t i m e s t e p ) , s q r t ( a . x * a . x + a . y * a . y + a . z * a . z ) FROM (SELECT s i m u l a t i o n , x , y , z , t i m e s t e p , p r e s s u r e FROM [ TABLE ] WHERE s i m u l a t i o n <= [ SIM NUM ] ) AS a INNER JOIN (SELECT s i m u l a t i o n , x , y , z , max ( p r e s s u r e ) AS max FROM [ TABLE ] WHERE s i m u l a t i o n <= [ SIM NUM ] GROUP BY s i m u l a t i o n , x , y , z ) AS b ON a . s i m u l a t i o n = b . s i m u l a t i o n AND a . p r e s s u r e = b . max AND a . x = b . x AND a . y = b . y AND a . z = b . z GROUP BY A . s i m u l a t i o n , a . x , a . y , a . z ; " "

Figure 7 :

 7 Figure 7: 1D Model Range Query Results

Figure 8 :

 8 Figure 8: Aneurysm Range Query Results Figures 9 and 10 show the results for queries Hemolab 1and Hemolab 2 respectively. The Hilbert Space-filling curve method offers the best option among the multidimensional array representations. For instance, it is up to 60% faster for Hemolab 2 than the Precision Elimination method and the Histograms Base Partitioning method. The dense and balanced data distribution together with the lower dimensionality enable analyses to be executed more efficiently.

Figure 9 :

 9 Figure 9: Hemolab1 Query Results

Figure 10 :

 10 Figure 10: Hemolab2 Query Results

Figure 11 :

 11 Figure 11: Loading Time Figure 12 shows the results for the range queries evaluating data of the 1D model in a distributed settings. The results observed in the sequential scenario discussed in Section 6.3 remain valid. In this context, the PE method is again more efficient and the Histogram Base Partitioning method reads unnecessary data.PostgreSQL shows better results than most other methods, except P.E., as the amount of data under consideration is very small, and the system is able to use a combination of indexes to retrieve a minimum amount of disk pages to answer the query. However, as the next results show, Post-greSQL is not able to maintain fast execution time when we increase the amount of data.Figure13shows the result obtained for the range query evaluating a region of the 3D model enclosing an aneurysm. For this query, and all following analyses, PostgreSQL has the worst performance by a great margin. Both MonetDB and SciDB are tailored for analytical workloads, and offer

Figure 12 :

 12 Figure 12: 1D Model Range Query Results

Figure 13 :

 13 Figure 13: 3D Model Aneurysm Range Query Results

Figure 14 :

 14 Figure 14: Hemolab 1 Query Results

  . The charts show the CPU usage and Disk I/O accomplished during execution of Hemolab 1 query. SciDB and MonetDB exhibit completely different patterns. MonetDB's CPU usage varies much more, and is interleaved with high disk I/O peaks. SciDB maintains a steady CPU usage throughout most of the execution time, and performs much less disk I/O than MonetDB.

Figure 15 :Figure 16 :

 1516 Figure 15: Hemolab 2 Query Results

Acknowledgment

This work has been partially funded by CNPq, CAPES, FAPERJ and Inria (MUSIC projects) and the European Commission (HPC4E H2020 project) and performed (for P. Valduriez) in the context of the Computational Biology Institute (www.ibc-montpellier.fr) and for (F. Porto, H. Lustosa and P. Blanco) in the context of the INCT-MACC project (http://macc.lncc.br).


Marco Congedo
email: marco.congedo@gmail.com

Ronald Phlypo

Alexandre Barachant

A Fixed-Point Algorithm for Estimating Power Means of Positive Definite Matrices

Keywords: Power Mean, Geometric Mean, High Dimension, Riemannian Manifold, Symmetric Positive-Definite Matrix. I

The estimation of means of data points lying on the Riemannian manifold of symmetric positive-definite (SPD) matrices is of great utility in classification problems and is currently heavily studied. The power means of SPD matrices with exponent p in the interval [-1, 1] interpolate in between the Harmonic (p = -1) and the Arithmetic mean (p = 1), while the Geometric (Karcher) mean corresponds to their limit evaluated at 0. In this article we present a simple fixed point algorithm for estimating means along this whole continuum. The convergence rate of the proposed algorithm for p = ±0.5 deteriorates very little with the number and dimension of points given as input. Along the whole continuum it is also robust with respect to the dispersion of the points on the manifold. Thus, the proposed algorithm allows the efficient estimation of the whole family of power means, including the geometric mean 1 .

INTRODUCTION

The study of means (centers of mass) for a set of symmetric positive definite (SPD) matrices has recently attracted much attention, driven by practical problems in radar data processing, image and speech processing, computer vision, shape analysis, medical imaging (especially Diffusion Magnetic Resonance Imaging and Brain-Computer Interface), sensor networks, elasticity, mechanics, quantum entanglement, numerical analysis and machine learning (e.g., [START_REF] Barachant | Multi-Class Brain Computer Interface Classification by Riemannian Geometry[END_REF][START_REF] Bhatia | Positive Definite Matrices[END_REF][START_REF] Congedo | Riemann Geometry: an Universal BCI Classification Framework[END_REF][START_REF] Fillard | A Riemannian Framework for the Processing of Tensor-Valued Images[END_REF][START_REF] Li | Riemannian Distances for Signal Classification by Power Spectral Density[END_REF][START_REF] Moakher | The Riemannian Geometry of the Space of Positive-Definite Matrices and Its Application to the Regularization of Positive-Definite Matrix-Valued Data[END_REF]). In many applications the observed data can be conveniently summarized by SPD matrices, for example, some form of their covariance matrix. In others, SPD matrices arise naturally as kernels, density matrices, elements of a search space, etc. Regardless the nature of the data generating them, SPD matrices may be treated as points on a smooth Riemannian manifold in which natural geometrical notions such as distance along the geodesic between two points and the center of mass (mean) between two or more points are naturally defined [START_REF] Bhatia | Positive Definite Matrices[END_REF]. In turn, these notions allow useful operations such as interpolation, smoothing, filtering, approximation, averaging, signal detection, and classification. In the latter application, a simple minimum distance to mean (MDM) classifier has been shown to perform as state-of-the-art classifiers in the challenging field of braincomputer interfaces [START_REF] Congedo | Riemann Geometry: an Universal BCI Classification Framework[END_REF]. Thence, given a cloud of data points, their mean on the manifold plays an important role as the natural point where a tangent space to manipulate them should be constructed, allowing the application of sophisticated non-Riemannian classification algorithms within the Riemannian framework. For these reasons, among many others, the numerical estimation of means on Riemannian manifold is currently a subject of intense investigation.

A. The Manifold of Symmetric Positive-Definite matrices

In differential geometry, a smooth Riemannian manifold M is a topological space equipped with an inner product on the tangent space defined at each point and varying smoothly from point to point. The tangent space TGM at point G is the Euclidean vector space containing the tangent vectors to all curves on M passing through G (Fig. 1).

Figure 1: Schematic representation of the SPD manifold, the geometric mean G of two points and the tangent space at G. Consider two points (e.g., two covariance matrices) C1 and C2 on M. The geometric mean of these points is the midpoint on the geodesic connecting C1 and C2, i.e., it minimizes the sum of the two squared distances δ 2 (C1, G)+δ 2 (C1, G). Now construct the tangent space TGM at G. There exists one and only one tangent vector ζ1 (respectively ζ2) departing from G and arriving at the projection of C1 (respectively C2) from the manifold onto the tangent space; we see that the geodesics on M through G are transformed into straight lines in the tangent vector and that distances are mapped logarithmically; the map from the manifold (symmetric positive definite matrices S++) to the tangent space (symmetric matrices S) is of logarithmic nature. The inverse map from the tangent space to the manifold is of exponential nature [START_REF] Bhatia | Positive Definite Matrices[END_REF].

B. The Geodesic

The SPD manifold has nonpositive curvature and is complete; for any two points C1 and C2 on M, a unique path on M of minimal length (at constant velocity) connecting the two points always exists. The path is named the geodesic and the points along it satisfy

    1 1 1 1 2 2 2 2 1 t 2 1 1 2 1 1 # , t 0,1 t C C C C C C C   , (1)
with t the arc-length parameter [START_REF] Bhatia | Positive Definite Matrices[END_REF] (Fig. 1). With t=0 we are at C1, with t=1 we are at C2 and with t=1/2 we are at the geometric mean of the two points (denoted G in Fig. 1). As a special case we note I#tC=C t . Geodesic equation (1) verifies C1#tC2=C2#1-tC1 and (C1#tC2) -1 =C1 -1 #tC2 -1 . The points along the geodesic can be simply understood as the means of C1 and C2 weighted by t according to the Riemannian metric, in analogy with the weighted mean according to the Euclidean metric given by (1-t)C1+tC2.

C. The Distance

Given two matrices (points) C1 and C2 dimension N·N on M, their Riemannian distance is the length of the geodesic (1) connecting them. It is given by [START_REF] Bhatia | Positive Definite Matrices[END_REF] (Fig. 1)

      11 22 22 1 2 1 2 1 n n , F C C Ln C C C tr Ln Λ ln       , (2)
where  is the diagonal matrix holding the N eigenvalues λ1,…,λN of matrix CC  . This distance has a remarkable number of properties, some of which are listed in [START_REF] Congedo | Approximate Joint Diagonalization and Geometric Mean of Symmetric Positive Definite Matrices[END_REF]. Besides the obvious symmetry and positivity, particularly useful in signal processing are the following invariance properties (for any invertible B):

Congruence     2 1 2 1 , , TT BC B BC B C C   , ( 3
)
Self-Duality     11 12 2 1 , , C C C C    . ( 4
  2 k k k k kk ; arg min F G G w w C G w C     C A , (5)
    1 2 1 1 1 k k k k kk ; arg min F G G w w C G w C         C H . (6)

E. The Geometric Mean of a Matrix Set

In M the w-weighted Geometric Mean GG (w;C) is the point realizing the minimum of Σkwkδ 2 (Ck, GG), where the Riemannian distance function δ(,) acting on M has been defined in [START_REF] Bhatia | Positive Definite Matrices[END_REF]. The geometric mean is the unique point on M satisfying non-linear matrix equation [START_REF] Moakher | A differential geometric approach to the arithmetic and geometric means of operators in some symmetric spaces[END_REF][START_REF] Lawson | Weighted means and Karcher equations of positive operators[END_REF]  

11 22 kk k 0 w Ln G C G    GG . (7)
In general, it has closed-form solution only for K=2, in which case it is indeed equal to C1#½C2 (indicated shortly as C1#C2)see [START_REF] Barachant | Multi-Class Brain Computer Interface Classification by Riemannian Geometry[END_REF] and Fig. 1 -furthermore, for K=2 it is the unique solution to Riccati equation (C1#C2)C2 -1 (C1#C2)=C1 [START_REF] Bhatia | Positive Definite Matrices[END_REF] and is also equal to B -1 D1 ½ D2 ½ B -T for whatever joint diagonalizer B of C1 and C2, i.e., for whatever matrix B satisfying BC1B T =D1 and BC2B T =D2, with D1, D2 invertible diagonal matrices [START_REF] Congedo | Approximate Joint Diagonalization and Geometric Mean of Symmetric Positive Definite Matrices[END_REF].

The geometric mean enjoys all 10 properties postulated in the seminal work [START_REF] Ando | Geometric means[END_REF]. These properties can be found listed in [START_REF] Lawson | Weighted means and Karcher equations of positive operators[END_REF][START_REF] Lawson | Monotonic properties of the least squares mean[END_REF]. We remind here the congruence invariance ( 8) and self-duality ( 9), inherited directly from the corresponding properties of its associated distance function (see ( 3) and ( 4)):

    1 K 1 K ; , , ; , , T T T BG w C C B G w BC B BC B  GG , (8)     1 1 1 1 K 1 K ; , , ; , , G w C C G w C C     GG . (9)

F. Power Means

In the same spirit of section D. and E., given again C={C1,…,CK} and w={w1,…,wK}, we can generalize to SPD matrices the power mean of real numbers (generalized means) in the continuum p[-1,1] such as [START_REF] Lawson | Weighted means and Karcher equations of positive operators[END_REF][START_REF] Lawson | Monotonic properties of the least squares mean[END_REF][START_REF] Lim | Matrix Power means and the Karcher mean[END_REF], [START_REF] Ando | Geometric means[END_REF] where, for any SPD matrices Ck, GP #pCk with p(0,1] is the mean of GP and Ck weighted by p (1), C -1 ={C1 -1 ,…,CK -1 } and GG (w;C) is the geometric mean of section E. GP (w;C;p) is named the w-weighted power mean of order p [START_REF] Lawson | Weighted means and Karcher equations of positive operators[END_REF][START_REF] Lawson | Monotonic properties of the least squares mean[END_REF][START_REF] Lim | Matrix Power means and the Karcher mean[END_REF]. Following [START_REF] Ando | Geometric means[END_REF], the pair of power means obtained at symmetric values of p around zero are the dual of each other; for a negative value of p the mean is defined as the inverse of the mean for -p as applied on the inverted input matrices C -1 . Power means interpolate continuously in between the harmonic mean (p = -1), the geometric mean (p = 0, evaluated at the limit) and the arithmetic mean (p = 1). Thereupon, the family of the power means encompasses and generalizes all Pythagorean means we have encountered so far. All of them enjoy the congruence invariance as the geometric mean (8), but their duality (second line of (10)) coincides with the selfduality property (9) only for p = 0. The numerous properties of the power means can be found in [START_REF] Lim | Matrix Power means and the Karcher mean[END_REF].

                kp k 11 ; ; p # , p 0,1 ; ; p ; ; p , p 1,0 ; ;0 , , p = 0 k G w w G C G w G w G w G w                C CC CC PP PP PG ,

II. ALGORITHMS FOR POWER MEANS

A. Motivation

We sought a general algorithm for computing the w-weighted power mean of order p, with p(-1, 1)\{0}. We are also interested in an effective algorithm for estimating the geometric mean (p=0). The most popular algorithm for computing the geometric mean is its gradient descent flow with fixed step size, however its convergence rate deteriorates rapidly as the SNR decreases. The same is true for the method based on approximate joint diagonalization in [START_REF] Congedo | Approximate Joint Diagonalization and Geometric Mean of Symmetric Positive Definite Matrices[END_REF]. Second order methods have complexity growing very fast with the size of the input matrices, thus they are little useful in practical applications [START_REF] Jeuris | A survey and comparison of contemporary algorithms for computing the matrix geometric mean[END_REF]. The algorithm proposed in [START_REF] Zhang | A Majorization-Minimization Algorithm for the Karcher Mean of Positive Definite Matrices[END_REF] has high complexity per iteration and slow linear convergence rate. For a review see both [START_REF] Jeuris | A survey and comparison of contemporary algorithms for computing the matrix geometric mean[END_REF] and [START_REF] Congedo | Approximate Joint Diagonalization and Geometric Mean of Symmetric Positive Definite Matrices[END_REF]. In contrast, the complexity per iteration of the algorithm we propose is low and its (linear) convergence rate is fast and almost uniform for p =  ½. Here uniform convergence rate means that the number of iterations does not increase with the size nor with the number of the input set matrices, a very useful property in practical applications. Importantly, the convergence rate is robust with respect to the dispersion of the points (noise) in the whole interval p(-1, 1)\{0}, as we will see.

B. A General Multiplicative Fixed-Point Algorithm

Hereafter it will be convenient to lighten notation; let us denote the weighted power mean of order p as P, which by ( 10) is equal to GP (w;C;p) if p(0, 1] or to (GP (w;C -1 ;-p)) -1 if p[-1,0). This way we will need to handle only one expression for whatever value of p(-1, 1)\{0}, such as

  ** ; ; p P G w  C P ( 11
)
where *=sgn(p) is the dual operator. As initialization we take as customary the closed form solution of the mean in the case when all matrices in set C all pair-wise commute. It has been proved in [START_REF] Lim | Matrix Power means and the Karcher mean[END_REF] (Property 1, p. 1502) that the mean is then

  1 0 p p kk k P w C   . ( 12
)
Let us now turn to the algorithm. We write ( 11) using ( 1) to obtain

  * 2 2 2 2 p ** kk k P P w P C P P          . ( 13
)
Equation ( 13) as a unique SPD solution [START_REF] Lim | Matrix Power means and the Karcher mean[END_REF]. Numerical experiments show that iterating this expression as it is (hereafter referred to as "naive fixed-point") results in a rather slow convergence rate, maximal for abs(p)=½, but slower and slower as abs(p) get closer to 0 or to 1. In order to fasten convergence we design a multiplicative algorithm as it follows: post-multiplying both sides of (13) by P -*/2 and taking the inverse at both sides we obtain

From [START_REF] Jeuris | A survey and comparison of contemporary algorithms for computing the matrix geometric mean[END_REF] we see that upon convergence H=I. H here plays the role of the origin in the SPD manifold M for data linearly transformed by P -*/2 . In particular, the identity matrix I is the point of symmetry in M corresponding to 0 in the Euclidean space; as P -1/2 is a whitening matrix for the arithmetic mean (p=1), so P -*/2 is a whitening matrix for the whole family of power means along continuum p[-1, 1]. We wish to proceed by multiplicative updates according to [START_REF] Zhang | A Majorization-Minimization Algorithm for the Karcher Mean of Positive Definite Matrices[END_REF]. Rather than to P * itself, we thus seek an algorithm converging to P -*/2 , which is its inverse square root for * = 1, i.e., when p(0, 1] and its square root for * = -1, i.e., when p[-1,0). The numerical stability of fixed-point iterates ( 14) is ensured by the fact that H converges toward I. Moreover, using our update rule any update matrix with form H -φ in ( 14) is equivalent to H -1 upon convergence. We have observed that by using a fixed exponent φ from the first iteration we can reach convergence rate faster than the one achieved by the naive fixed point in the whole interval p(-1, 1)\{0}. For doing so, we have to take φ inversely proportional to p as per

  1 1 2 φ ε p , 2 ε 1,    , ( 16
)
where ε is a constant eccentricity parameter for hyperbola ( 16) delimiting an acceptable interval for φ as a function of |p|. At the lower limit ε = 1 we have a unit hyperbola φ = ½ / abs(p). By increasing ε toward 2 we obtain faster convergence, up to a certain value, which depends upon the signal to noise ratio. In this study we take ε as 4/3 and we keep it fixed in all simulations, obtaining quasi-uniform convergence rate for p =  ½ in a wide range of situations, as we will show next by means of simulations.

Algorithm MPM (Multiplicative Power Means)

INPUT: p(-1, 1)\{0}, K positive weights w={w1,…,wK} such that Σkwk=1 and K SPD matrices C * ={C1 * ,…,CK * }, with *=sgn(p). OUTPUT: P, the w-weighted Power Mean of order p. BEGIN Initialize X as the principal square root inverse of ( 12) if p(0,1] or as its principal square root if p[-1,0). Set ζ equal to a small floating precision number (e.g., 10 -10 ) Set φ = 0.375/abs(p)

REPEAT   p * w T kk k H XC X      X H X    UNTIL 1 N F HI   RETURN     p 0, 1 , p 1,0 T T X if P Y Y Y X if          END.

III. SIMULATION STUDIES

In many engineering applications, the matrix condition number of the SPD matrices summarizing the data (observations, recordings,…) tends to be positively correlated with the number of sensors. Also, the dispersion in the manifold of the matrices is proportional to the noise level. The following generative model for input data matrices {C1,…,CK} of size NN is able to reproduce these properties:

  TT k k k k k C UD U V E V I     , (17)
where -the signal part is given by UDkU T , where U is a matrix with elements drawn at random at each simulation from a uniform distribution in [-1,1] and then normalized so as to have columns with unit norm and Dk are K diagonal matrices with diagonal elements dk,n randomly drawn at each simulation from a squared Gaussian distribution with expectation=1/2 n , where n{1,…,N} is the index of the N diagonal elements, thus forming elements of a well-known geometrical series absolutely converging to 1. The elements of the series represent the energy of N source processes, thus their sum is supposed finite.

-The uncorrelated noise part is given by I, where I is the identity matrix and  here is taken as 10 -6 ;

-The structured noise part is given by VkEkVk T , where the Vk matrices are generated as U above, the Ek matrices are generated as Dk above and  is a constant controlling the SNR of the generated points through

      E E T k k T k k k k tr UD U SNR tr V E V I          , ( 18
)
where E denotes expectation. In the ensuing simulations we study relevant outcome parameters as a function of the SNR, which is inversely proportional to noise level as per (18), as well as a function of the size (N) and number (K) of input matrices. We compare the gradient descent algorithm for the geometric mean (section II.A), the naive fixed point algorithm for power means given in [START_REF] Lim | Matrix Power means and the Karcher mean[END_REF] (section II.B) and the MPM algorithm here presented, the latter for several values of p. The gradient descent algorithm is known to converge only up to a certain value of dispersion, i.e., when the SNR is sufficiently high [START_REF] Afsari | On the convergence of gradient descent for finding the Riemannian center of mass[END_REF]. The stopping criterion is chosen identical for all algorithms: the relative error of matrix P with respect to a reference matrix Pref is a dimensionless measure defined as:

22 ref ref FF P P P  . ( 19
)
As stopping criterion, considering two successive iterates P(i-1) and P(i), we use

    2 1 1 i i 1 N F P P I    , ( 20
)
whose magnitude does not depend on the size nor to the norm of the matrices. Figure 2 shows the typical convergence behavior for the geometric mean gradient descent algorithm, the naive algorithm with p=0.5 and the MPM algorithm (p=0.5 and p=0.001), for N=20, K=100 and SNR={100, 10, 1, 0.1}. Figure 3 shows the main effects (bars) and their standard deviation (sd: lines) across 50 simulations of N={10, 25, 50}, K={10, 100, 500} and SNR={100, 1, 0.01} on the number of iterations. "Main effects" means that for each level of N, K and SNR the average and sd are computed across all levels of the other two variables, as in a classical analysis of variance (ANOVA). We see that the number of iterations required by the MPM algorithm is comparable in all situations for p=0.5 (quasi-uniform convergence behavior), whereas for the geometric mean gradient descent it grows rapidly with N and with noise, as it is well known [START_REF] Congedo | Approximate Joint Diagonalization and Geometric Mean of Symmetric Positive Definite Matrices[END_REF]. Even with p=0.001, which is a very small value, the convergence rate of the MPM is superior to the convergence rate of the geometric mean gradient descent flow, which often do not converge at all (the maximum number of iterations allowed was fixed to 50 for all algorithms. See Fig. 2 for examples of divergence). The simulation in Fig. 4 (left) shows the relative error (19) of the power mean with respect to the geometric mean (p=0) for several small values of p, N=10, K=100 and SNR={100, 1, 0.01}. We see that the approximation to the geometric mean in term of log-relative error is linear in log-p, ceiling to 10 -5 with high noise. Figure 4 (right) shows the main effects (across 50 simulations) of N={10, 25, 50}, K={50, 100, 500} and SNR={100, 10, 1} on the error for p=0.001. We see that the error decreases with K and increases with N and with noise.

Notice that an error much lower than this can be achieved by computing the power mean for two small symmetric values around zero and then computing their geometric mean by the closed-form solution for two matrices (see section E.); in practice a relative error < 10 -3 is fully satisfactory, since the noise level in challenging data is much higher than that.

IV. CONCLUSIONS

Power means are generalized means interpolating continuously in the interval p[-1, 1], with p=1 yielding the arithmetic mean, lim p0 yielding the geometric mean and p=-1 yielding the harmonic mean. We have presented a new multiplicative algorithm for estimating power means in the interval p[-1, 1]\{0}. A numerical analysis shows that its convergence rate is very fast and quasi-uniform for values of p close to ½ and -½, while for values of p close to 0 it is still faster as compared to the gradient descent with fixed step-size used for estimating the geometric mean. Furthermore, it converges also in low SNR situations. We can therefore prefer the proposed algorithm also for approximating the geometric mean. In conjunction with expressions ( 5) and ( 6), with the MPM algorithm we can now estimate a number of means sampling along the continuum p=[-1, 1], as in Fig. 5. We name such a sampling a Pythagorean mean field. Applications of mean fields include the possibility to evaluate the most appropriate choice of mean depending on its use and on the data at hand. Mean fields also allow the extension of current Riemannian classifiers, such as in [START_REF] Barachant | Multi-Class Brain Computer Interface Classification by Riemannian Geometry[END_REF][START_REF] Congedo | Riemann Geometry: an Universal BCI Classification Framework[END_REF][START_REF] Fillard | A Riemannian Framework for the Processing of Tensor-Valued Images[END_REF][START_REF] Li | Riemannian Distances for Signal Classification by Power Spectral Density[END_REF], for instance combining current classifiers as applied to all the means of the field. The application of mean fields to real data will be the object of future investigations.

Figure 2 :

 2 Figure2: Typical convergence behavior (on abscissa the number of iterations and on the ordinate the convergence as defined in (20) expressed in dB) for the geometric mean gradient descent, naive fixed point with p=0.5 and the MDM algorithm with p={0.5, 0.001}, for N=20 (dimension of input matrices), K=100 (number of input matrices) and SNR={100, 10, 1, 0.1} (18).

Figure 3 :

 3 Figure 3: main effects average (50 simulations) and sd number of iterations for N={10, 25, 50}, K={10, 100, 500} and SNR={100, 1, 0.01} for the geometric mean gradient descent and the MPM algorithm with p={0.5, 0.4, 0.001}.

Figure 4 :

 4 Figure 4: Relative Error to the geometric mean obtained with the MPM algorithm. Left: N=10, K=100, SNR={100, 1, 0.01} and p={10 -5 ,…, 10 -1 }. Right: Main effects (bars) and sd (lines) obtained across 50 repetitions for p=0.001, N={10, 25, 50}, K={10, 100, 500} and SNR={100, 1, 0.01}.

Figure 5

 5 Figure 5 is a TraDe plot (log-trace vs. log-determinant) for a sampling of power means along continuum p[-1, 1], showing the log-log linear relation of the trace and determinant of power means along the Pythagorean continuum.

Figure 5 :

 5 Figure 5: TraDe plot obtained with N=10, K=10 and SNR=1, for power means corresponding to p=1(Arithmetic), 0.5, 0.1, 0 (Geometric), -0.1, -0.5 and -1 (Harmonic). The relationship between the trace and the determinant of power means is log-log linear.

This research has been partially supported by European project ERC-2012-AdG-320684-CHESS