Thi-Van-Anh Nguyen
email: thi-van-anh.nguyen@unicaen.fr

Laurent D'orazio
email: laurent.dorazio@univ-bpclermont.fr

Sandro Bimonte
email: sandro.bimonte@irstea.fr

Jérôme Darmont
email: jerome.darmont@univ-lyon2.fr

Cost Models for View Materialization in the Cloud

Keywords: H.2.4 [Database Management]: Systems-Distributed databases, Query processing Algorithms, Design, Economics, Performance Cloud computing, Cost models, Materialized views, Performance and cost optimization

In classical databases, query performance is casually achieved through physical data structures such as caches, indexes and materialized views. In this context, many cost models help select a "best set" of such data structures. However, this selection task becomes more complex in the cloud. The criterion to optimize is indeed at least two-dimensional, with the monetary cost of using the cloud balancing query response time. Thus, we define in this paper new cost models that fit into the pay-as-you-go paradigm of cloud computing. These cost models help achieve a multi-criteria optimization of the view materialization vs. CPU power consumption problem, under budget constraints. Finally, we present experimental results that provide a first validation of our contribution and show that cloud view materialization is always desirable.

INTRODUCTION

Cloud computing has recently emerged as a new challenging paradigm. Cloud computing may be defined as a payper-use model for enabling on-demand access to reliable and configurable resources that can be quickly provisioned and released with minimal management. Users only pay for the resources they use. They need not set up the infrastructure, nor buy the software. These features make the design of cloud data management systems desirable. Unfortunately, performance in cloud computing usually relies on scalability, using bigger instances (scale-up) or more nodes (scale-out), making the bill more expensive.

On the other hand, performance optimization in databases has been intensively studied for decades. Reusing methods such as indexing, buffering or view materialization in the cloud could help improving performances without increasing costs. For instance, materialized views, which we particularly focus on in this paper, physically store the results of some relevant, frequent queries. One of the main challenges is then to select what views to materialize. Traditionally, the criteria to consider in the selection process mainly include view storage and maintenance costs [START_REF] Aouiche | Data mining-based materialized view and index selection in data warehouses[END_REF][START_REF] Baril | Selection of materialized views: A cost-based approach[END_REF].

In cloud computing, storage is virtually infinite, enabling to store all possible materialized views. However, each materialized view implies an extra storage cost. Thus, the performance optimization problem becomes finding the best tradeoff between raw scalability (i.e., increasing resources) and materialized views under budget constraints. To illustrate the complexity of selecting materialized views in the cloud, let us introduce a simple, fictitious example. Let the storage cost be $0.10 per GB and per month, and the computing cost be $0.24 per hour. A 500 GB dataset is stored in the cloud for a month. Let Q, the monthly query workload, process in 50 hours. Then, storage cost is $50, computing cost is $12, for a total of $62. If some materialized views are used, let us assume that workload processing time becomes 40 hours. Thus, computing cost becomes $9.6. However, materialized views use up additional storage space, e.g., 50 GB. Thus, storage cost becomes $55, for a total cost of $64.6. Overall, performance has improved by 20%, but cost has also increased by 4%.

The selection of views to materialize depends on several parameters, such as the workload (query execution time, database size, etc.), the cloud pricing model (CPU, storage, ). This choice is a trade-off between response time and cost, and depends on the needs of a particular user. At one end of the spectrum, users under a strong budget constraint may accept long response times, whereas at the other end, users may disregard cost if they need fast response.

In this paper, we address the multi-criteria optimization problem of selecting a set of materialized views while optimizing the global monetary cost of storing and querying a database in the cloud. To achieve this goal, our main contribution is the design of cost models for storing, maintaining and querying data with materialized views over a cloud architecture.

The remainder of this paper is organized as follows. In Section 2, we provide the background information that is used throughout the paper. In Sections 3 and 4, we define cost models for cloud data management and materializing views, respectively. In Section 5, we describe the optimization process that is based on these cost models. In Section 6, we present an experimental evaluation and the first performance analyses of our models. In Section 7, we discuss the state of the art and compare it to our approach. Finally, in Section 8, we conclude this paper and hint at future research directions.

BACKGROUND

We present in this section the background information related to view materialization in the cloud. We first introduce a simple fictitious use case that serves as a running example throughout this paper. Then, we describe a typical pricing model in the cloud, illustrated by some of Amazon Web Services (AWS). Then, we briefly recall the principle of view materialization.

Running example

To illustrate our work, we rely on a simulated dataset storing the sales of an international supply chain. Business users need to analyze the total profit per day, month, and year; and per administrative department, region, and country. Table 1 provides an excerpt of this dataset.

Our full dataset stores 10 years (2000-2010) of sale data. Its size is 500 GB. We run over this dataset a query workload Q that includes such queries as Q1 = "sales per year and country", whose processing time is 0.2 hour. The size of the Q's result is 10 GB. A typical materialized view we may consider to optimize overall response time is V1 = "sales per month and country", whose processing time is 0.1 hour. The whole set of selected materialized views is denoted V . V 's size is 50 GB. Finally, the times to process Q with and without exploiting V are 40 hours and 50 hours, respectively.

Cloud pricing policies

Cloud Service Providers (CSPs) supply a pool of resources, such as hardware (CPU, storage, networks), devel-opment platforms and services. There are many CSPs on the market, such as Amazon, Google and Microsoft. Each CSP offers different services and pricing. This paper relies on a limited, yet representative enough model that includes the main, commonly billed elements, i.e., CPU, storage and bandwidth consumption. In order for the reader to have an overview of such a pricing policy, the following examples present a simplified version of AWS's offer.

In AWS, Elastic Compute Cloud (EC2) [1] provides computing resources. Different instance configurations can be rent (micro, small, large, extra large, etc.) at various prices, as illustrated in Table 2. For example, the cost for a small instance (1.7 GB RAM, 1 EC2 Compute Unit, 160 GB of local storage under Windows) is $0.12 per hour. Let us apply this pricing model onto our use case running on two small instances. Without materialized views, processing time is 50 hours, and thus computing cost is 2 × 0.12 × 50 = $12. With materialized views, processing time falls down to 40 hours and computing cost down to 2 × 0.12 × 40 = $9.6.

Instance configuration Price per hour

Micro

$0.03 Small $0.12 Large $0.48 Extra large $0.96 ... ...

Table 2: EC2 computing prices

Bandwidth consumption is billed with respect to data volume (Table 3). In this model, input data transfers are free, whereas output data transfer cost varies with respect to data volume, with an earned rate when volume increases. When applying this pricing model onto our use case, the cost of bandwidth consumption (query result of 10 GB) is (10 -1) × 0.12 = $1.08.

Data volume

Price per month Input data Any input data Free Output data First 1 GB Free Up to 10 TB $0.12 per GB Next 40 TB $0.09 per GB Next 100 TB $0.07 per GB ... ...

Table 3: Amazon bandwidth prices

Finally, the Amazon Simple Storage Service (S3) [2] supplies storage capabilities. In this model, the price varies with respect to data volume, with an earned rate when volume increases (Table 4). In our running example, monthly storage price when not using materialized views (500 GB dataset) is 0.14 × 500 = $70, and 0.14 × (500 + 50) = $77 when using materialized views (additional 50 GB storage).

Data volume Price per month

Materialized views

In Database Management Systems, a view is a virtual table associated to a query answer. Views help indirectly save complex queries, present the same data in different forms, support logical independence and reinforce security by masking some pieces of data from unauthorized users. Materializing a view, i.e., storing it physically into a table, further helps improve response time by avoiding to recompute the corresponding query each time the view itself is queried. However, materialized views must be refreshed when source data are updated, which induces some maintenance overhead.

In this work, we assume that we have at our disposal a set of candidate materialized views that have already been generated by an existing materialized view selection method (e.g., [START_REF] Baril | Selection of materialized views: A cost-based approach[END_REF]). We aim at choosing the best candidates with respect to the cloud's pay-as-you-go model, taking pricing constraints into account before any view materialization.

DATA MANAGEMENT COST MODELS

This section presents general cost models for data management in the cloud, i.e., without considering the use of materialized views. In cloud computing, customers rent resources to a CSP to run some applications. 1 recalls the costs involved (Section 2.2), i.e., bandwidth consumption for input data transfers and query result retrieval, data storage, and applications' processing time. Let Cc be the sum of computing costs, Cs be the sum of storage costs and Ct be the sum of data transfer costs. Then, the total cost C for cloud data management is:

C = Cc + Cs + Ct.
(1)

Let us define the general parameters that we use to express our cost models (Table 5), as well as two functions. Function s() returns the size in GB of any of these parameters, e.g. s(DS) is the size of the dataset. Function ts() returns the storage time of any dataset, e.g. ts(DS) is the storage time of the dataset in the cloud.

Parameter

Data transfer cost

Data transfer cost depends on several parameters: the size of the initial dataset, the amount of data related to queries and query results, the volume of inserted new data and the pricing model applied by the CSP. The total data transfer cost Ct is the product of the CSP's atomic transfer cost ct by the total size of transferred data:

Ct = ( m i=1 (s(Ri) + s(Qi)) + (s(DS) + s(insertedData)) × ct.
(2) If we adopt Amazon's EC2 like pricing model, which does not charge for input data transfers, we may ignore input queries, the initial data set and inserted data. As a consequence, the total data transfer cost becomes:

Ct = m i=1 s(Ri) × ct.
(3) Moreover, recall that Amazon EC2's cost is variable. It is null for the first GB. Then, it becomes $0.12 up to 10 TB, and so on (Section 2.2).

Example 1: In our running example, with 10 GB (with the first one being free) of bandwidth consumption, data transfer cost is: Ct = s(RQ) × ct = (10 -1) × 0.12 = $1.08.

Computing cost

Let us pose that queries are executed on computing instances {ICj}j=1..n. Each instance may bear different performances (with respect to its number of CPUs, its available RAM, etc.), and thus different costs. The cost for renting instance ICj is denoted c(ICj). It must be paid at each connection to the cloud. Processing time of query Qi on instance ICj is denoted tij. Then, the processing cost of running the set of queries Q = {Q1, ..., Qm} is:

Cc = m i=1 tij × c(ICj). ( 4
)
Example 2: In our running example, let us consider that the query workload Q is processed on two small instances. Then, its processing cost is: Cc = t × c(IC) + t × c(IC) = RoundU p(50) × 0.12 + RoundU p(50) × 0.12 = 50 × 0.12 + 50×0.12 = $12. We must use a function to round processing time up because every started hour is charged.

Storage cost

Storage cost depends on parameters such as the CSP's pricing policy, the size of the data (original dataset and in-serted data) and the storage time. We assume that the storage period in the cloud is devided into intervals. In each interval, the size of the stored data is fixed. The total storage cost is the CSP's storage cost function cs(s(DS)) (where DS is the data storage size, that is to say the size of the original data set and the inserted data) multiplied by the sum of sizes of the initial dataset and inserted data multiplied by their respective storage time during the intervals.

Cs = intervals cs(DS) × (t end -tstart) × s(DS) (5)
Where tstart, t end are start point and end point of an interval.

Example 3: Using Amazon S3 for storage pricing (referred to table 4) and considering that 0.5 TB (512 GB) data has been stored for 12 months. At the beginning of the eighth month, we insert 2 TB (2048 GB) of new data in the cloud. Thus we have two intervals. The storage cost is:

Cs = intervals cs(DS) × (t end -tstart) × s(DS) = (512 × 0, 14 × (7 -0) + (512 + 2048) × 0, 125 × (12 -7) = $2131.76.

COST MODELS FOR MATERIALIZING VIEWS IN THE CLOUD

This section presents cost models for materializing views in the cloud, relying on the costs models developed in Section 3. We assume here that queries are executed on a constant number, nbIC of identical instances IC. In future work, we shall consider the evaluation process on multiple, variable instances.

Let V cand = {V k } k=1..p be a set of candidate materialized views output by any existing materialized view selection technique.

Data transfer cost

We assume that materialized views are selected on the client's side, outputting final set of materialized views V that are materialized in the cloud. This tactic helps save bandwidth and benefit from the computing performance of the cloud. With materialized views created in the cloud, transfer costs of materialized views are null. As a consequence, total transfers cost Ct is not impacted and remains expressed by Formula 3.

Computing cost

Using materialized views implies modifying computing costs, since query processing may exploit materialized views, and views must be materialized and maintained. Computing cost expressed by Formula 4 thus becomes:

Cc = CprocessingQ + CmaintenanceV + C materializationV (6)
where CprocessingQ is the cost of processing queries, CmaintenanceV is the cost of materialized views' maintenance, and C materializationV is the cost of materializing views.

Materialization cost

If a view is materialized, its associated query must be executed, which must be paid for in the cloud. Let the materialization time of view V k be t materialization (V k ). We assume that all views are materialized on a constant number of identical instances, nbIC .

The total materialization time is:

T materializationV = p k=1 t materialization (V k ). (7)
According to Formula 4, materialization cost is:

C materializationV = T materializationV × c(IC) × nbIC . ( 8
)
Example 4: In our running example with two nodes, we assume that V = {V1}, the total materialization time is

T materializationV = p k=1 t materialization (V k ) = t materialization (V1) = 1 hour. Then, materialization cost is: C materialisationV = t materialisation (V1) × c(IC) × nbIC = 1 × 0.12 × 2 = $0.24.

Query processing cost

When using materialized views, query processing time is defined by two main parameters: query workload Q and set of materialized views V . Queries may use the contents of materialized views instead of recomputing their result. Note that we consider that Q is fixed. Since views are usually materialized with respect to a given workload and ours is fixed, then V is also fixed.

Let tiV be the processing time of query Qi when exploiting the set of materialized views V .Thus, the total processing time of Q against V is:

TprocessingQ = m i=1 tiV . (9)
Example 5: In our running example, total processing time when using V is 40 hours.

After taking view materialization into account, Formula 4 becomes:

CprocessingQ = TprocessingQ × c(IC) × nbIC . ( 10
)
Example 6: So, the computing cost becomes: CprocessingQ = m i=1 tiV × c(IC) × nbIC = 40 × 0.12 × 2 = $9.6.

Maintenance cost

The maintenance cost of materialized views is directly proportional to the time required for updating materialized views when they are impacted by modifications of the source dataset. Note that we consider that querying and maintenance do not occur at the same time. For example, queries are posed during day-time and maintenance is performed during night-time. Let the maintenance time of view V k be tmaintenance(V k ). Then, the total maintenance time of V is:

TmaintenanceV = p k=1 tmaintenance(V k ). ( 11
)
Example 7: The maintenance time in the running example is 5 hours.

As a consequence, the cost of maintaining the set of materialized views V is:

CmaintenanceV = TmaintenanceV × c(IC) × nbIC . ( 12
)
Example 8:

Thus, the maintenance cost is:

CmaintenanceV = TmaintenanceV × c(IC) × nbIC = 5 × 0.12 × 2 = $1.2.

Storage cost

Using materialized views does not impact the storage cost as presented by Formula 5. Exploiting materialized views to enhance query performance implies storing them in the cloud and pay the corresponding cost. As a consequence, some data can be duplicated. In that case DS represents the data set, the duplicated data due to materialized views and possible inserted data. Note that original data and materialized views are stored for the whole considered storage period.

Example 9: In our running example, the data set has been stored for a year, the size of duplicated data due to materialized views is 50 GB. In addition, no data are inserted during the considered period. Thus, we have a single interval and storage cost is Cs = ((500 + 50) × (12 -0) × 0.14 = $924

OPTIMIZATION PROCESS

In this section, we investigate how the use of materialized views in the cloud impacts query performance, while trying to minimize overhead storage cost. We propose an algorithm to select the best set of materialized views by exploiting the cost models introduced in Section 4.

Objective functions

Base on the ideas in [START_REF] Kllapi | Schedule optimization for data processing flows on the cloud[END_REF], we distinguish in this section three scenarios, i.e., three objective functions, with respect to the needs and capacity of customers (budget limit, response time limit).

Budget limit

Given a predefined financial budget Bl, our objective in this scenario (denoted M V1) is to select, from the set of candidate materialized views V cand , the set of views to materialize in the cloud V that minimizes query processing time (Figure 2: The chosen solution are colored in red).

M V1 = Minimize TprocessingQ C = Cc + Ct + Cs ≤ Bl ( 13
)
Figure 2: Minimizing processing time under budget constraint

Response time limit

Given a predefined response time limit T l, our objective in this scenario (denoted M V2) is to select, from the set of candidate materialized views V cand , the set of views to materialize in the cloud V that minimizes financial cost (Figure 3).

M V2 = Minimize C = Cc + Ct + Cs TprocessingQ ≤ T l (14)
Figure 3: Minimizing cost under response time constraint

Response time vs. cost tradeoff

In this scenario (denoted M V3), our objective is to select, from the set of candidate materialized views V cand , the set of views to materialize in the cloud V that offers the best tradeoff between query processing time and financial cost (Figure 4). To provide the user with control over this process, we introduce a weight parameter on processing time (α) and cost (1 -α).

M V3 = Minimize (α × TprocessingQ + (1 -α) × C) (15)

Optimization algorithm

In order to select the set of views to materialize, we rely on an existing algorithm such as [START_REF] Baril | Selection of materialized views: A cost-based approach[END_REF] enabling to obtain a set of candidate views V cand . Then we solve the Knapsack 0/1 problem [START_REF]Hard knapsack problems[END_REF] considering this set. To solve this problem, we have opted for a dynamic programming approach.

As a input, we provide elements of our cost models. For example, in the case of MV1 (budget limit), the result is: V = Knapsack(V cand , FCt, FCc, FCs, Bl), where FCt, FCc, FCs are functions which represent the parameters of the models described in the section 4. In particular, FCt involves transfer cost parameters (section 4.1), FCc involves computing cost parameters(section 4.2), and FCs involves storage cost parameters (section 4.3).

EXPERIMENTAL VALIDATION

In this section, we describe the overall setup of our preliminary experimentation effort and the results we have obtained from it. This experiment is realized on both sides: in the cloud and on the client's side. The dataset and materialized views are stored in the cloud, and queries are processed in the cloud. The view selection algorithms are implemented in Java on the client's side. The idea is to select views on the client's side and to materialize them in the cloud.

Experimental setup

The experiments were conducted on a virtual cluster composed of five virtual machines with a 50 GB disk, 2 GB of RAM and 4 vCPU. This virtual cluster is deployed on a physical architecture consisting of two bi-pros 493 GHz with 32 GB of RAM, interconnected by a 1 GB/s network in a 8x2 DRS cluster; and a bi-pro 2,526 GHz with hyperthreading enabled, outside of the DRS cluster, with 16 logical cores with 72 GB of RAM on a 1 GB/s network. All machines feature Hadoop (version 0.20.2) [3] and Pig Latin (version 0.7.0) [START_REF] Olston | Pig latin: a not-so-foreign language for data processing[END_REF]. Queries are written in Pig Latin and are executed by the Pig compiler on the Hadoop cluster within the MapReduce framework [START_REF] Dean | Mapreduce: simplified data processing on large clusters[END_REF]. On the client's side, algorithms are implemented on a PC Dual Core 2 T6660 2.2GHz, with 2 GB of RAM.

We used in our experiments a subset of the dataset from our case study (Section 2.1), whose size is 10 GB. The workload we applied on the dataset was made of 10 queries that calculate the total profit per day, month, year and per country, department, and region such as "per year and per country".

Experimental results

We experimented the three scenarios M V1, M V2 and M V3 described in Section 5.1, with a variable workload made of 3, 5, and 10 queries, with and without exploiting materialized views. Our results are plotted on Figure 5. They clearly show that, in all scenarios, creating materialized views in the cloud is desirable.

In scenario M V1, a significantly better response time is achieved with materialized views, under the same budget constraint than without (Figure 5(a)). The improved performance rates are shown in the Table 6. Similarly, in scenario M V3, materialized views help achieve a tradeoff (chosen by the client) between cost and response time, whether the priority is put on cost (Figure 5(c)) or response time (Figure 5(d)). Improved tradeoff rates are shown in the Table 8.

Number of queries

Number of queries

RELATED WORK

Although there has been a lot of research related to the cloud for a couple of years, relatively few approaches relate to the optimization of data management in the cloud. In this section, we provide a brief overview of cost models, databaselike features and the use of optimization techniques in the cloud.

Cost models have been proposed in the cloud for scheduling of dataflows with regard to monetary cost and/or completion time [START_REF] Kllapi | Schedule optimization for data processing flows on the cloud[END_REF], and cost amortization of data structures to ensure the economic viability of the provider [START_REF] Kantere | Predicting cost amortization for query services[END_REF], particularly for self-tuned caching [START_REF] Dash | An economic model for self-tuned cloud caching[END_REF]. Other cost models have been developed for a real-life astronomy application using the Amazon cloud fee structure [START_REF] Berriman | The application of cloud computing to astronomy: A study of cost and performance[END_REF][START_REF] Deelman | The cost of doing science on the cloud: the montage example[END_REF]. The cost performance tradeoffs of different execution and resource provisioning plans have been simulated, showing that by provisioning the right amount of storage and computing resources, cost can be significantly reduced with no significant impact on application performance [START_REF] Deelman | The cost of doing science on the cloud: the montage example[END_REF]. The performance of three workflow applications with different I/O, memory and CPU requirements has also been compared on Amazon EC2 and a typical high-performance cluster (HPC) to identify what applications achieve the best performance in the cloud at the lowest cost [START_REF] Berriman | The application of cloud computing to astronomy: A study of cost and performance[END_REF].

Recent research takes interest in various aspects of database and decision support technologies in the cloud. For example, different studies investigate the storage and processing of structured data [START_REF] Chatziantoniou | Asset queries: a declarative alternative to mapreduce[END_REF], the optimization of join queries, and how to support analysis operations such as aggregation [START_REF] Condie | Online aggregation and continuous query support in mapreduce[END_REF]. Cloud data warehousing and OLAP systems also raise various problems related to storage and query performance [START_REF] Mahboubi | Enhancing xml data warehouse query performance by fragmentation[END_REF]. Adaptations of these technologies to the cloud's specificities are also addressed, e.g., pay-as-you-go pricing and elasticity [START_REF] Armbrust | A view of cloud computing[END_REF], or the calculation of OLAP cuboids using the MapReduce runtime environment [START_REF] Sergey | Applying map-reduce paradigm for parallel closed cube computation[END_REF]. Finally, although massive scale distributed database systems such as BigTable have begun to incorporate indexes and materialized views [START_REF] Cafarella | Data management projects at google[END_REF], the many techniques developed for maintaining materialized views in traditional, relational databases [START_REF] Agrawal | Asynchronous view maintenance for vlsd databases[END_REF], such as incremental view maintenance [START_REF] Ceri | Deriving production rules for incremental view maintenance[END_REF], deferred maintenance [START_REF] Zhou | Lazy maintenance of materialized views[END_REF] and distributed view maintenance [START_REF] Luo | A comparison of three methods for join view maintenance in parallel rdbms[END_REF], have not been incorporated into the cloud yet. This is thus the main incentive for our own approach.

CONCLUSION

In this paper, we proposed an approach for decreasing the cost of data management in the cloud, by using a classical database performance optimization technique, i.e., view materialization. Our main contributions are novel cost models that complement the existing materialized view cost models with a monetary cost component that is primordial in the cloud. Then, we exploit these cost models into an optimization process that achieves a tradeoff between computing power enhancement and view materialization, under budget constraints. Our first experimental validation hints that view materialization is always desirable.

Many perspectives are opened by this preliminary research. We indeed have many plans on improving our cost models. First, we should include pricing models from several CSPs but Amazon, to render our cost models more versatile. Second, it is well-known that optimization techniques are the most efficient when combined, e.g., when materialized views and index are jointly exploited [START_REF] Aouiche | Data mining-based materialized view and index selection in data warehouses[END_REF]. Thus, we should incorporate indexing, caching and/or fragmentation into our approach. Third, we should exploit the specificities of cloud languages (such as HiveQL [START_REF] Thusoo | Hivea petabyte scale data warehouse using hadoop[END_REF], Jaql [START_REF] Beyer | Jaql: A scripting language for large scale semistructured data analysis[END_REF] and Pig Latin [START_REF] Olston | Pig latin: a not-so-foreign language for data processing[END_REF]) with respect to semantics and parallelism to further improve performance. Finally, in this paper, we worked on one single cloud instance, thus minimizing the primordial elasticity characteristic of the cloud. Thus, we should expand our cost models on variable resources.

Another line of future work relates to the experimental validation of our proposals. Here, due to limited resources, we experimented our cost models on a small cluster and a toy dataset. Thus, we plan to design a wider-scale experimentation plan that would exploit a full-fledged database or data warehouse benchmark, such as TPC-E [4] or the Star Schema Benchmark [START_REF] O'neil | The star schema benchmark and augmented fact table indexing[END_REF], onto AWS or other infrastructure providers.

Figure 1 :

 1 Figure 1: Costs involved in cloud data management

  Description DS Dataset Q = {Qi}i=1..m Query workload R = {Ri}i=1..m Set of query answers

Figure 4 :

 4 Figure 4: Optimizing response time and cost

Figure 5 :

 5 Figure 5: Experimental results of application of materialized views

Table 1 :

 1 Sales dataset excerpt network, etc.

	Year Month Day Country	Region	Department Profit
	2000	12	31	France	Auvergne Puy-de-Dôme $35.000
	2000	1	1	France	Auvergne Puy-de-Dôme $40.000
	2000	12	31	Italy	Campanie	Naples	$23.000
	1999	1	1	Italy	Campanie	Naples	$50.000

Table 4 :

 4 Amazon storage prices

	First 1 TB	$0.14 per GB
	Next 49 TB	$0.125 per GB
	Next 450 TB	$0.11 per GB
	...	...

Table 5 :

 5 General parameters

Table 6 :

 6 Improved performance rates between with and without materialized views under the same budget limit in the scenario M V1In scenario M V2, a significantly lower global cost is achieved with materialized views, under the processing time constraint than without (Figure5(b)). The improved cost rates are shown in the Table7.

		Budget limit IP Rate
	3	$0.8	25%
	5	$1.2	36%
	10	$2.4	60%

Table 7 :

 7 Improved cost rates between with and without materialized views under the same time limit in the scenario M V2

Table 8 :

 8 Improved tradeoff rates between with and without materialized views in the scenario M V3

		Rate (α = 0.3) Rate (α = 0.7)
	3	55%	32%
	5	50%	35%
	10	68%	45%

Acknowledgment

This work is sponsored by the ANR under grant SYSEO ANR-10-TECSAN-005-01. We would like to sincerely thank all the colleagues in IRSTEA, LIMOS and ERIC laboratories for the interesting discussions and the reviewers for their extremely insightful remarks for improving this paper.


Oriol Pont
email: oriol.pont@inria.fr

Antonio Turiel

Hussein Yahia

Singularity analysis of digital signals through the evaluation of their Unpredictable Point Manifold

Keywords: multiscale measures, reconstructible systems, MSC-94A08, MSC-43A32, MSC-46T12, CCS-G.1

The local singularity exponents of a signal are directly related to the distribution of information in it. This fact implies that accurate evaluation of such exponents opens the door to signal reconstruction and characterisation of the dynamical parameters of the process originating the signal. Many practical implications arise in a context of digital signal processing, since the information on singularity exponents is usable for compact encoding, reconstruction and inference. Since singularity exponents are conceptually associated to differential calculus, its evaluation in a digital context is not straightforward and it requires the calculation of the Unpredictable Point Manifold of the signal. In this paper, we present an algorithm for estimating the values of singularity exponents at every point of a digital signal of any dimension. We show that the key ingredient for robust and accurate reconstructibility performance lies on the definition of multiscale measures in the sense that they encode the degree of singularity and the local predictability at the same time.

Introduction

In digital signals of all kinds, the estimation of their local singularity exponents provides a valuable information that can be used to compactly encode the signal values [START_REF] Mallat | Singularity detection and processing with wavelets[END_REF][START_REF] Mallat | Wavelet transform maxima and multiscale edges[END_REF]. Furthermore, the singularity exponents can be accurately estimated at the top points of a Wavelet Transform Modulus Maxima (WTMM) line, for the case of turbulent flows or analogous scale-invariant systems [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Mallat | Characterization of signals from multiscale edges[END_REF][START_REF] Muzy | Wavelets and multifractal formalism for singular signals: Application to turbulence data[END_REF]. However, the generalization of this methodology from the WTMM points to any point of the signal (at its finest resolution level) is not straightforward: such generalizations become convoluted and they provide rather imprecise estimates of singularity exponents outside the WTMM lines [START_REF] Struzik | Determining local singularity strengths and their spectra with the wavelet transform[END_REF] and poor description of their distributions [START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data: a comparative study[END_REF].

There is a completely different approach based on local gradient modulus measures that are numerically estimated by convolution with wavelet or other positive multiscale kernels. This method provides more stable and accurate results for discretised signals even at maximum resolution [START_REF] Turiel | Relevance of multifractal textures in static images[END_REF][START_REF] Turiel | The multi-fractal structure of contrast changes innatural images: from sharp edges to textures[END_REF][START_REF] Turiel | The self-similarity properties of natural images resemble those of turbulent flows[END_REF]. This approach provides the basic bricks of the so-called Microcanonical Multiscale Formalism (MMF) [START_REF] Pont | On optimal wavelet bases for the realization of microcanonical cascade processes[END_REF][START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF], which gives access to unveil key physical descriptors for turbulence and turbulent-like systems [START_REF] Isern-Fontanet | Microcanonical multifractal formalism: application to the estimation of ocean surface velocities[END_REF][START_REF] Pont | Empirical evidences of a common multifractal signature in economic, biological and physical systems[END_REF][START_REF] Pont | Application of the microcanonical multifractal formalism to monofractal systems[END_REF][START_REF] Pont | Description, modeling and forecasting of data with optimal wavelets[END_REF][START_REF] Turiel | Role of multifractal sources in the analysis of stock market time series[END_REF][START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF][START_REF] Turiel | Multifractal method for the instantaneous evaluation of the stream function in geophysical flows[END_REF][START_REF] Turiel | Tracking oceanic currents by singularity analysis of micro-wave sea surface temperature images[END_REF]. This way, estimating the singularity exponents of a signal also provides information about the physical process behind it.

At the local level, the accuracy of MMF-based singularity-analysis estimates is significantly sensitive to the shape of the multiscale measure used to obtain them [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF], even if the estimated distributions are not [START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data: a comparative study[END_REF]. Therefore, if we want to characterize only the statistical properties of the singularity exponents we can project the signal over any well-behaved kernel, such as positive multiscale functions with fast enough decay [START_REF] Turiel | Relevance of multifractal textures in static images[END_REF], but if we need a precise geometrical identification of the singularity structures and their underlying patterns from edges to textures, then we require a fine-tuned multiscale measure.

The work we present in the following describes a method that allows high-performance determination of the singularity exponents, with an algorithm that is fast, computationally cheap, stable, accurate and provides fine resolution. Its definition relies in the connection of singularity exponents with the concept of signal reconstruction from the Unpredictable Point Manifold (UPM). Singularity analysis establishes a link between singularity exponent and reconstructibility degree of a point, in a sense that the UPM can be estimated as the Most Singular Component (MSC) of the signal [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF].

To validate the detection of the information content provided by the singularity exponents as calculated from UPM-measures, we illustrate a comparison of them with the unrelated method of multiscale entropy [START_REF] Costa | Multiscale entropy analysis of complex physiologic time series[END_REF][START_REF] Starck | Multiscale entropy filtering[END_REF]. We observe that all the features and their informative degree coincide and furthermore the singularity exponents provide a finer resolution and so a better localization of the information.

This paper is presented as an extended version of our previous findings that we have presented at the 14th International Workshop on Combinatorial Image Analysis, IWCIA 2011 held in Madrid in May 2011 [START_REF] Pont | An optimized algorithm for the evaluation of local singularity exponents in digital signals[END_REF]. Here we summarise our previous findings, while at the same time we expand them accordingly: we give a deeper description of the analysis performed and the methods used, and we present a more exhaustive evaluation (which includes both signal reconstructibility and illustrating a comparison of information content with multiscale entropy). Namely, since analysis of multidimensional signals may also be of interest to the reader, we present here a generalisation of our technique to nD and give a wider insight to it by illustrating the singularity analysis and reconstruction performance for the cases of empirical signals of 1D (stock-market series), 2D (Meteosat IR map) and 3D (ERA-40 global climate re-analysis).

The paper is structured as follows: in the next Section we introduce the basics of singularity analysis. Section 3 explains the concept of reconstruction from the MSC. Section 4 is the core of the paper: first the conditions defining UPM-measures are discussed, then Subsection 4.1 explains the definitions for calculus on reduced neighbourhoods through the cross-gradient and cross-reconstruction operators, and finally in Subsection 4.2 our algorithm to retrieve the UPM and thus the singularity exponents from digital signals is presented. Some analyses on empirical signals of 1D, 2D and 3D are shown in Section 5. The last section, Section 6, presents the conclusions of our work.

Evaluating singularity exponents

In mathematical analysis, the concept of singularity exponents refers to the analysis of singularities of differentiable functions in a broad sense and different non-equivalent defintions are possible; in the present work we focus on its usual meaning in the theory of complex systems. Using singularity analysis we intend to characterise the local behaviour of a R m -valued function f (x) defined on R d around each one of its domain points x according to their so-called singularity exponents, Hölder exponents [START_REF] Jaffard | Multifractal formalism for functions. I. Results valid for all functions[END_REF] or Hurst exponents [START_REF] Jones | Wavelet packet computation of the hurst exponent[END_REF][START_REF] Simonsen | Determination of the hurst exponent by use of wavelet transforms[END_REF]. If the signal behaves at point x according to:

f (x + r) -f (x) = α(x)r h(x) + o r h(x) (r → 0) (1)
then h(x) is the singularity exponent at x: small perturbations around x lead to changes in the function value as a power of the displacement modulus r = r . A strictly n-times derivable function naturally leads to an integer Hölder exponent h(x) = n. As a result, this formulation allows to generalise the concept of integer differentiability to fractional differentiability, but to complete that transposition, a slightly more exigent formulation of eq. ( 1) is required: we need to assume that there exists a (1, 1) continuous tensor from R d to R m , α(x), such that

f (x + r) -f (x) = α(x)|r r h(x)-1 + o(r h(x) ) (r → 0) (2)
where α(x)|r is the standard duality bracket for (1, 1) tensors. If such representation is possible, the exponent h(x) is called the Hurst exponent of the function, but if only eq. ( 1) applies then we use the term Hölder exponent.

In real-world empirical signals, Hölder exponents often appear masked by long-range correlations and fluctuations due to noise and discretisation effects that obstruct their direct estimation [START_REF] Turiel | The multi-fractal structure of contrast changes innatural images: from sharp edges to textures[END_REF][START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF]. A generalization that avoids such problems consists of defining singularity exponents from gradient measures [START_REF] Turiel | The multi-fractal structure of contrast changes innatural images: from sharp edges to textures[END_REF]. Let s(x) be a signal defined in R d with values in R, we define the gradient-based measure µ through its density dµ(x):

dµ(x) = ∇s (x) dx (3)
which is defined as absolutely continuous with respect to the Lebesgue measure. As a consequence, the measure of any Borelian A is [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF]:

µ(A) = A dx ∇s (x) (4)
The main advantage of gradient measures is that they characterise the local degree of singularity of the signal [START_REF] Turiel | The multi-fractal structure of contrast changes innatural images: from sharp edges to textures[END_REF]. From eq. ( 1), let us now consider a function f (x) that has h(x) + 1 as Hölder exponents (whith a +1 shift introduced for later convenience). Let Br(x) be a ball of radius r (for some arbitrary norm in R d ) centred around x. This measure behaves as:

µ (Br(x)) = α(x) r d+h(x) + o r d+h(x) (r → 0) ( 5
)
where d is the dimension of the signal domain. For discrete signals, we can interpolate this measure by means of wavelet projections [START_REF] Lee | Signal smoothness estimation in hölder spaces[END_REF][START_REF] Lee | Noisy image enhancement on smooth function spaces[END_REF]. Let Ψ be a wavelet function, T Ψ µ(x, r) is the projection of the measure µ at the point x and scale r [START_REF] Daubechies | Ten Lectures on Wavelets, Cbms-Nsf Regional Conference Series in Applied Mathematics[END_REF][START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]:

T Ψ µ(x, r) = R d dµ(x ′ ) 1 r d Ψ x -x ′ r (6)
This way, T Ψ is an operator that maps from a set M of σ-finite measures on R d to the set of

R d × R + → R functions.
The singularity exponents of the signal as per eq. ( 5) are directly expressed on the wavelet projection too [START_REF] Daubechies | Ten Lectures on Wavelets, Cbms-Nsf Regional Conference Series in Applied Mathematics[END_REF][START_REF] Turiel | The multi-fractal structure of contrast changes innatural images: from sharp edges to textures[END_REF]:

T Ψ µ(x, r) = α Ψ (x) r h(x) + o r h(x) (r → 0) (7)
with the same singularity-exponent values h(x).

Measures can then be used on discretised data with any multiscale function Ψ in eq. ( 7), without any requirement of the admissibility property [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF]. As the attainable resolution level depends on the number of zero-crossings of the wavelet used [START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data: a comparative study[END_REF], hence we can even use positive kernels that maximise the resolution at which we estimate the singularity exponents [START_REF] Turiel | Relevance of multifractal textures in static images[END_REF].

A reconstruction formula from the most singular set

Signal reconstructibility from a partial set of its gradient is the key concept to define a wavelet that precisely identifies the singularity exponents. The theoretical development and its practical implementation give raise to a reconstruction algorithm [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF] relating reconstructibility to singularity.

We require that every point of the signal has a singularity exponent defined [START_REF] Frisch | Turbulence: The legacy of A.N. Kolmogorov[END_REF][START_REF] Parisi | On the singularity structure of fully developedturbulence[END_REF] and that these exponents organise in sets that follow the signal multiscale structures with a geometrical hierarchy [START_REF] Falconer | Fractal Geometry: Mathematical Foundations andApplications[END_REF]. As a consequence, such structures follow cascading relations whose distribution matches those of the macroscopic descriptors of the system under study [START_REF] She | Universal scaling laws in fully developed turbulence[END_REF]. Since it is difficult for classical methods to consistently assess accurate values of singularity exponents on empirical signals, these structures have been often characterised in a statistical way only. However, it is possible to have a direct access to this geometrical hierarchy through the reconstruction of the full structure from its vertex [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF]. In multifractal signals, the set associated to the vertex is formally well-known and it corresponds to the Most Singular Component (MSC), i.e., the set of points having the most singular values of singularity exponent [START_REF] Turiel | The multi-fractal structure of contrast changes innatural images: from sharp edges to textures[END_REF][START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF].

For the gradient of a signal over the MSC to fully reconstruct the signal, it needs to concentrate all its information [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF]. We look for a universal operator able to reconstruct a signal that follows the known statistical symmetries of turbulent-like and multiscale signals [START_REF] Frisch | Turbulence: The legacy of A.N. Kolmogorov[END_REF]. Such operator, if it exists, it is uniquely defined under the requirements of being deterministic, linear, translational invariant, isotropic and leading to the known power-spectrum shape.

Given a multiscale signal s, F∞ is its MSC: F∞ is the set of points x such that h(x) ∈ ] h∞ -∆, h∞ + ∆ [ with h∞ being the minimum value of h(x) over the finite domain of the discrete signal and ∆ is an ideally infinitessimal threshold parameter. Let ∇ F ∞ s be the essential gradient of s defined as:

∇ F ∞ s(x) = ∇s(x) δ F ∞ (x) (8)
where δ F ∞ consists of a delta distribution on the continuum of the F∞: a distribution that is homogeneous in Hausdorff topological dimension to a repartition in between dimensions d -1 and d, so that it assigns a uniform weight inside the MSC points F∞ while it vanishes outside.

If we denote the gradient reconstructing functional by G, a deterministic representation of the gradient of the signal means:

∇s(x) = G ∇ F ∞ s(x) (9)
Considering now the fact of G being linear, an integral representation is given by:

∇s(x) = F∞ dx ′ G(x, x ′ ) ∇s(x ′ ) (10)
The property of translational invariance implies:

∇s(x) = F∞ dx ′ G(x -x ′ ) ∇s(x ′ ) ( 11
)
where G is the R d × R d density of the linear operator G. By considering the differential forms involved, since ∇s(x) is a gradient in the distributional sense, G can also be represented as a gradient of some functional-vector g:

s(x) = F∞ dx ′ g(x -x ′ ) • ∇s(x ′ ) ( 12
)
which actually is a convolution between g and ∇ F ∞ s, eq. ( 8), expressed in the Fourier domain as:

ŝ(k) = ˆ g(k) • ∇ F ∞ s(k) (13)
which is the reconstruction formula [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF].

The vector field g is the universal reconstruction kernel. As mentioned, we require that it is isotropic and does not alter the power spectrum (e.g., a f -2 signal has a white gradient). Isotropy means it has the same direction as the unit vector of the frequency, namely: ˆ g(k) = ĝ(k) u k where k is the frequency modulus: k = k u k . Then the modulus of the reconstruction kernel, ĝ(k), is:

ˆ g(k) = i u k k (14)
with i = √ -1 as the imaginary unit, as we require it to work also in the case where the MSC is the whole domain. This definition implies that the reconstruction kernel g is actually an inverse gradient operator. The reconstruction is always perfect if we consider the whole domain, i.e., if F∞ = R d then eq. ( 13) reduces to a trivial identity: ∇ F ∞ s = ∇s. It is not evident at all if there exists a smaller set F∞ ⊂ R d for which eq. ( 13) returns a perfect reconstruction, at least without assuming any further conditions. Nevertheless, any set F that leads to a perfect reconstruction with the reconstruction formula verifies that [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF]:

div ∇ F c s = 0 ( 15
)
where F c is the complementary set of F . As the gradient and divergence operators are local and linear, the formula above implies that whether a point needs to be included into F or not is something that can be evaluated locally at the neighbourhood of that point. In other words, the points that must be in F to get a perfect reconstruction or unpredictable points are those whose value cannot be predicted from their neighbours. Moreover, signal predictability is often considered a core concept for the analysis of complex systems [START_REF] Aurell | Predictability in the large: an extension of the concept of lyapunov exponent[END_REF][START_REF] Boffetta | Predictability: a way to characterize complexity[END_REF], by often considering measures of information growth such as e.g., Lyapunov exponents and Kolmogorov-Sinai entropy. We call the Fu set of unpredictable points the Unpredictable Points Manifold (UPM) of the signal, which is by definition the smallest set possible for having perfect eq. ( 13) reconstruction. Also, we are considering the hypothesis that Fu = F∞ for the signals under study, which is the basic conjecture giving rise to the framework of reconstructible systems [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF]. Most practical experiences evidence that properly calculated MSCs are able to reconstruct the signals [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF].

Unpredictable point manifold (UPM) and its associated measures

The concept of UPM-measure generalizes that of gradient measures introduced in the previous section. The first requirements for any singular positive UPM-measure µ are:

i) It is concerned with the local singularity structure of functions.

ii) It leads to a MSC as close to the UPM as possible.

In essence, we can think of a UPM-measure as a gradient measure that is sensitive to the local predictability degree, eq. ( 15). They represent the gradient form as finite differences in discrete signals and in a way that highlights unpredictability. Vector-valued wavelet projections are able to combine both features at the same time when carefully tuned.

Some standard singularity-analysis methods perform multiple wavelet projections of the measure to estimate the singularity exponents through log-log regressions on eq. ( 7). However, projecting at several scales penalizes resolution -especially for the most singular structures [START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data: a comparative study[END_REF]-and it increases computational complexity; that is why we concentrate on punctual estimates of the singularity exponents from the UPM measure [START_REF] Pont | Application of the microcanonical multifractal formalism to monofractal systems[END_REF][START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data: a comparative study[END_REF], namely:

h(x) = log (T Ψ µ(x, r 0 )/ T Ψ µ(•, r 0 ) ) log r 0 + o 1 log r 0 ( 16
)
The denominator normalising the main logarithm T Ψ µ(•, r 0 ) is the expectation of the wavelet projection (in practice, its average over the whole signal or as wide as stationarity can be assumed) and lessens the relative amplitude of the fluctuations due to the o 1 log r 0 correction term. To effectively apply eq. ( 16) we require that r 0 is small enough so that the leading order dominates and we can neglect this correction term. The scale r 0 is thus the smallest accessible one, i.e., the pixel scale in practical terms. Without loss of generality, we conventionally assign a Lebesgue measure of 1 to the whole space domain, so for a

N 1 × N 2 × • • • × N d d-dimensional signal the value of r 0 is fixed to r 0 = 1 d √ N 1 •••N d
. This condition means that we require that the signal is large enough to make the first term in the right hand side of eq. ( 16) the dominant contribution to the singularity exponent. For practical purposes, any linear resolution r 0 finer than a hundredth is generally good enough.

The Cross Fourier Transforms

To maximise the attainable resolution of the singularity exponents function, we reduce the reconstruction formula, eq. (13, to the minimal meaningful neighbourhood of each point, i.e., its 2d nearest neighbours in 2d connexity neighbourhoods. These are expressed as the 4 neighbours at the ends of a cross for the case of 2D signals (d = 2) and that is why we talk about cross-gradient and gross-reconstruction operators.

In the general case and for any quantity p(x), we represent the neighbourhood of a given point x 0 by a (2d + 1)-component vector comprising this point and its nearest neighbours and following the indexing convention shown in Figure 1. The result is the vector

(p 0 , • • • , p 2d ).
The notion of local predictability directly extends to any dimension, independently of the number of neighbouring directions, with vectors of (2d + 1) components linearly growing with the dimension d.

Let us remark the fact that the nearest connexity in such cross-like structuring remains inherent to our processing of the (p 0 , • • • , p 2d ) vectors. We cannot just apply the standard harmonics of the Fourier Transform, i.e. (e 2ikπ/(2d+1) ) k , as if they were standard linear vectors. Instead, since the centre of the cross and the position of all the other components are related to ±1 pixel displacements, the basic Nyquist frequency in each direction is 2π/3. So that, the equivalent Fourier transform adapted to this cross formation has the following components:

j = e 2πi/3 = cos(2π/3) + i sin(2π/3) = - 1 2 + i √ 3 2 , j = j 2 (17)
With them, we define the direct Cross Fourier Transform of any (2d + 1)-vector p = (p 0 , • • • , p 2d ) as the complex (2d + 1)-vector ˆ p = (p 0 , • • • , p2d ) obtained according to the following formula: Figure 1. Schematic representation of the indexing of the points in the 2D cross (left) and the 3D case (right). The central point will be assigned the index 0, the next point in the positive direction of the first coordinate axis will be indexed 1, the opposite neighbour of it in the negative direction is indexed 2, the next point in the second coordinate axis, positive direction is indexed as 3, its opposite is indexed as 4, and so on we proceed equivalently for any other coordinate axes.

ˆ p = F p ( 18
)
where F is the following (2d + 1) × (2d + 1) complex matrix:

F = 1 3             1 1 1 1 • • • 1 1 1 1 j j 1 1 1 1 1 j j 1 1 1 1 1 1 1 j 1 1 1 . . . . . . . . . 1 1 1 1 j 1 1 1 1 1 1 1 j j 1 1 1 1 • • • 1 j j             (19)
which represents the linear combination of the harmonics associated to the nearest cross-like displacements. The inverse of this matrix is given by:

F -1 =             1 -2d 1 1 1 • • • 1 1 1 1 j j 0 0 0 0 1 j j 0 0 0 0 1 0 0 j 0 0 0 . . . . . . . . . 1 0 0 0 j 0 0 1 0 0 0 0 j j 1 0 0 0 • • • 0 j j             (20)
We now introduce surrogates of the gradient and reconstruction operators restricted to the cross-like nearest neighbourhood. The goal is to use these to evaluate the local degree of predictability of the central point, and that is why they are based on the Cross Fourier Transform that we have just introudced.

The Cross Gradient Operator is the operator

(∂ k1 , • • • , ∂ kd ) = F -1 • ( ∂ k1 , • • • , ∂ kd ) • F.
When expressed in Fourier space, the operator is just a component-by-component multiplication of a given function by the ∂ k1 • • • ∂ kd factors, defined as:

∂ kn = (0, • • • , 0, 0, i √ 3, -i √ 3, 2n-1,2n 0, 0, • • • , 0, 0) (21)
where the index convention for the vectors runs from 0 to 2d as before, and the index convention for the coordinate axes runs 1 to d (as for n).

The Cross Reconstruction Operator is one of the inverses of the Cross Gradient Operator (the gradient cancels any constant value so that the reconstruction is defined up to a constant). We canonically define it such that the resulting reconstructed vector is of zero mean along its (2d + 1) components: 2d i=0 s i = 0. Therefore, the algorithm subtracts the mean of the vector as a first step (see the next subsection).

The Cross Reconstruction is the operator R = F -1 • R • F which in Fourier has 2d functional components, namely R = ( Rk1 , • • • , Rkd ). The reconstruction operator acts as the sum of all the products of each component by the corresponding components of the inputted gradient vector on which it operates:

Rkn = (0, • • • , 0, 0, -i/ √ 3, i/ √ 3, 2n-1,2n 0, 0, • • • , 0, 0) (22)
analogously to the gradient case. Also here their operation can be considered analogous to diagonal matrices acting on the (2d + 1)-dimensional ˆ p vectors. The reconstructed vector is obtained as the sum of its component in each coordinate axis, from k1 to kd.

The Cross Gradient and the Cross Reconstruction are the two basic algorithms for the design of the UPM-measures. Let us summarise again the algebraic structure we are considering: For every point in a d-dimensional, scalar, digital signal we have a (2d + 1)-dimensional vector p made of the value at that point and its nearest neighbours. For every such vector we have a total of d gradient components (which are cross-like structured (2d + 1)-vectors themselves), one for each dimension. For each one of these d gradient components there is (2d + 1)-dimensional reconstruction operators and the sum of the resulting d components results in a reconstructed (2d + 1)-vector. Cross-derivation followed by cross-reconstruction is equivalent to projection to a (2d + 1)-identity matrix (up to a normalisation constant in the 0-frequency in the cross-Fourier domain). The same procedure applies when the considered digital signal is vector-valued itself; then we just independently treat every component separate.

A Local Correlation Singularity Measure

One of the simplest and most generic ways of measuring the local unpredictability consists of just subtracting the signal value at a given point from the value inferred from its neighbour points. These values are previously detrended to cancel any global offset influence. The result then measures local correlation and so it is called the Local Correlation Singularity Measure, which is calculated with the following algorithm: Goal: To evaluate T Ψ lcsm µ(x 0 , r 0 ) at a given point x 0 . Algorithm:

(1) The nearest neighbourhood of x 0 is expressed as a (2d + 1)-vector s = (s 0 , • • • , s 2d ) according to the cross scheme that we have presented in Figure 1. (2) This vector is conveniently detrended: let S = 1 2d-1 2d i=0 s i , we define the detrended vector p = (p 0 , • • • , p 2d ) as:

p 0 = s 0 + S ; p i = s i -S , i = 1, • • • , 2d (3)
We apply the Cross Gradient Operator to p; the result is the gradient vectors g kn . (4) The first component of these gradient vectors are kept for later use, An = g kn,0 (the zeroth/central component for the kn axis direction), n = 1, • • • , d. (5) We then set these components to zero, g k1,0 = • • • = g kd,0 = 0. (6) We apply the Cross Reconstruction Operator to the resulting vectors g kn and we sum them to obtain the reconstructed vector r. [START_REF] Falconer | Fractal Geometry: Mathematical Foundations andApplications[END_REF] We apply once more the Cross Gradient Operator on the reconstructed vector r to obtain ρ kn . [START_REF] Frisch | Turbulence: The legacy of A.N. Kolmogorov[END_REF] We define the Local Correlation Singularity Measure as the modulus of the difference of the cross gradients at the centre of the cross, namely:

T Ψ lcsm µ(x 0 , r 0 ) = (An -ρ kn,0 ) 2
Which in fact is equivalent of keeping the modulus of a certain vector-valued wavelet projection. ( 9) Finally,we obtain the singularity exponent h(x 0 ) directly through eq. ( 16).

Most singular features in digital signals

In the following, we show the application of the singularity analysis based on the Local Correlation Singularity Measure to empirical real-world signals of different dimensions.

We begin with a time series, i.e., a 1D case. It consists of the evolution of a stock market index: the Spanish IBEX-35 on a daily basis for 10 years. Singularity analysis has many applications in this case because the singularity exponents encode key aspects of the multiscale structure of the stock dynamics and identification of its MSC can be used not only in reconstruction but also in inference and stochastic forecast of future volatilities [START_REF] Pont | Description, modeling and forecasting of data with optimal wavelets[END_REF]. In Figure 2 we show the time series, its MSC and its reconstruction. Centre: MSC points of this time series. When projected to the same discrete resolution of the original signal it covers somewhat less than 30% of the points. Right: Reconstructed time series from the gradient of the original signal on the MSC, using the reconstruction formula. The big part of the dynamics is recovered including most of its relevant aspects, resulting in a quality reconstruction of 27.9 dB as PSNR.

The Local Correlation Singularity Measure provides equally good results on 2D signals.

We show examples of two different cases with different but related applications. The singularity exponents on natural images describe some kind of multiscale edges that provide a compact coding [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF]. When applied on geophysical maps the objective consists in giving the descriptors of the turbulent dynamics as singularity exponents define eddies and filaments as well as they encode the interscale energy injection known as turbulent cascade process [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF].

In Figure 3 we show an example of a natural scene and its reconstruction from its MSC. The singularity exponents are directly related to the local informativity degree and distribution of information in the image, as we can see when comparing them to the methodologically unrelated notion of multiscale entropy [START_REF] Costa | Multiscale entropy analysis of complex physiologic time series[END_REF][START_REF] Starck | Multiscale entropy filtering[END_REF] but with the advantage that they are sensibly finer and they are able to reach the resolution of 1 pixel, i.e., the same resolution of the original signal, Figure 4. In Figure 5 we show a case of a geophysical turbulence map: a Meteosat image in the infrared channel. Meteosat is a satellite specialised in the observation of clouds. Its infrared radiometry maps have two important advantages over visible: they are acquired 24 hours a day and they indicate the temperature of the clouds, thus permitting the identification of fronts and rain. Since clouds at a large scale are principally advected by the atmospheric flow, the structure of their global distribution describe the dynamics of atmospheric turbulence [START_REF] Chigirinskaya | Unified multifractal atmospheric dynamics tested in the tropics part 1: horizontal scaling and self organized criticality[END_REF]. Our results in this case show that the MSC captures the most of the signal dynamics, as in the previous cases.

To illustrate a 3D case, we rather consider a 2D+1 spatio-temporal signal, namely the time of evolution of 2D maps. Video constitutes a paradigm in that case and space-time correlations and transitions play a key role in the majority of their compression algorithms, but the spatial and temporal axis are not trivially comparable and isotropic methods need to be modified. Instead, we show a case of geophysical turbulence where typical velocities allow relating the spatial structure of turbulence to its temporal dynamics [START_REF] Chaté | Transition to turbulence via spatio-temporal intermittency[END_REF][START_REF] Frisch | Turbulence: The legacy of A.N. Kolmogorov[END_REF][START_REF] Houra | Spatio-temporal turbulent structures of thermal boundary layer subjected to non-equilibrium adverse pressure gradient[END_REF]. In Figure 6 we show several slices of a sequence showing the time evolution of global maps of air temperature at the atmospheric boundary layer. We also show the corresponding MSCs where singularity exponents have been evaluated on a space-time domain. This signal corresponds to the ERA-40 re-analysis product where global high resolution remote sensing variables are assimilated into an atmospheric level to estimate atmospheric circulation in a global and continuous way. We can see in the detail Figure 7 the comparison between the 2D singularity exponents, which take into account only the spatial structure, and the spatio-temporal singularity exponents, which appear more informative.

Conclusions

Having an accurate and robust estimation of the singularity exponents of a signal becomes essential in order to retrieve its key features and localise its information content, particularly in nonlinear and multiscale settings. Furthermore, the singularity degree as encoded by the local singularity exponent of a signal point is directly related to the exceptionality of that point in the signal or, in other words, how unpredictable it is or how difficult to reconstruct from the rest of the signal. In contrast, the regular, predictable points can be easily inferred from the singular ones, which means that the Most Singular Component (MSC) of the signal can be used for compact encoding, reconstruction or inference. In this work, we have introduced an algorithm to accurately retrieve the singularity exponent at each point of a given signal. This algorithm calculates singularity exponents through a discrete mask, which is constructed as a measure of both the singularity degree and the local reconstructibility of that point. This mask actually represents a discretised numerical equivalent of a wavelet basis. With this basis, we introduce the notion of the Unpredictable Point Manifold (UPM) of a signal, which provides a rigorous definition for an intuitive concept: the degree of singularity of a given point in a signal is directly related to the degree of how exceptional is the value of the signal on that point compared to its direct neighbourhood.

Additionally, we have detailed a singularity-analysis algorithm based on the local correlations of the signal [START_REF] Turiel | Method and system for the singularity analysis of digital signals[END_REF]. The algorithm introduces the so-called "Local Correlation Singularity Measure" as the basis to retrieve estimates of the singularity exponents with optimal accuracy and resolution. As a consequence, the signal can be reconstructed from its MSC points at high quality, which can be used in image compression [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF], also in other related applications such as assessment of streamlines in turbulent flows [START_REF] Turiel | Tracking oceanic currents by singularity analysis of micro-wave sea surface temperature images[END_REF][START_REF] Turiel | The multifractal structure of satellite temperature images can be used to obtain global maps of ocean currents[END_REF], detection of meteorological convection [START_REF] Turiel | Multiscale techniques for the detection of precipitation using thermal IR satellite images[END_REF] and identification of stock-market investment cycles [START_REF] Turiel | Role of multifractal sources in the analysis of stock market time series[END_REF].

We have shown the performance of the Local Correlation Singularity Measure algorithm on empirical signals of different types. We have illustrated a 1D case with stock-market series, where estimation of the UPM allows inference and forecasting methods [START_REF] Pont | Description, modeling and forecasting of data with optimal wavelets[END_REF], the 2D cases of natural scenes and the atmospheric dynamics as seen by Meteosat IR maps, and a 3D case of the type 2D-space+1D-time, being careful with the space-time correspondence, illustrated with a sequence of global air temperature maps from the ERA-40 global climate re-analysis, where we have observed how the singularity exponents taken in the spatio-temporal domain give more details in that case than the 2D singularity exponents evaluated independently.

Figure 2 .

 2 Figure 2. Left: Spanish blue chip index IBEX-35 opening daily values from 1990 to 2000.Centre: MSC points of this time series. When projected to the same discrete resolution of the original signal it covers somewhat less than 30% of the points. Right: Reconstructed time series from the gradient of the original signal on the MSC, using the reconstruction formula. The big part of the dynamics is recovered including most of its relevant aspects, resulting in a quality reconstruction of 27.9 dB as PSNR.

Figure 3 .

 3 Figure 3. Left: Original Lena image. Centre: MSC, defined as {h < -0.5}; it comprises 30% of the points of the image. Right: Reconstruction. Though some of the finest details are missing due to the lack of capability of the method to capture every UPM point, the obtained reconstruction exhibits a fine quality of 24.5 dB PSNR.

Figure 4 .

 4 Figure 4. Left: Singularity exponents of the Lena image estimated using the Local Correlation Singularity Measure; they are represented using a inverse grayscale palette (the brightest the smaller, so more singular). Right: Multiscale entropy of the Lena image, in the same palette. All details of the detected features persist, confirming the correspondence between singularity exponents and degree of informativity. Resolution is nevertheless coarser for the entropy estimation.

Figure 5 .

 5 Figure 5. Top: Infrared radiance captured by one of the Meteosat satellites on 31st Jul 1998 on the area 20 • E to 20 • W and 0 • to 20 • N. Middle: MSC, defined as {h < -0.5}; it comprises 20% of the points of the image. Bottom: Reconstruction formula applied on the MSC points. Key aspects of atmospheric turbulence as the structure-rich subtropical winds are detected in the MSC and the reconstructed image reproduces them.

Figure 6 .

 6 Figure 6. Left: ERA-40 re-analysis showing a sequence of global maps of the air temperature at the boundary layer of the atmosphere. Resolution is 2 • 5 in space and 6h in time. Temperature is expressed in Kelvins. Centre: MSC estimated from the singularity exponents in the spatio-temporal domain. Right: Reconstructions from the MSCs, capturing most of atmospheric circulation. Top to bottom: sequence slices at the 876576th, 876600th, 876624th and 876648th hour past midnight Jan 1st 1900.

Figure 7 .

 7 Figure 7. Top: 876576th hour slice of ERA-40 artificially rescaled 4x with bicubic interpolation for the purpose of clarity of illustration. Middle: singularity exponents calculated only in the space. Some rescaling artefacts visibly appear but without significant disturbance of the fine structure details. Bottom: singularity exponents calculated in the space-time domain. Notice the increased degree of detail when the temporal information is taken into account.

Acknowledgements

The authors acknowledge financial support from Fonds LEDUCQ, the ARC FIBAUR project and project MIDAS-6 AYA2010-22062-C05.