Pierre Matri
email: pmatri@fi.upm.es

Alexandru Costan
email: alexandru.costan@irisa.fr

Gabriel Antoniu
email: gabriel.antoniu@inria.fr

Jesús Montes
email: jmontes@fi.upm.es

María S Pérez
email: mperez@fi.upm.es

Towards Efficient Location and Placement of Dynamic Replicas for Geo-Distributed Data Stores

Keywords: Storage Networks, Geo-Replication, Wide-area replication, Content Distribution Network, Data Warehousing, Metadata, Data Consistency, Cloud, Availability

Large-scale scientific experiments increasingly rely on geodistributed clouds to serve relevant data to scientists worldwide with minimal latency. State-of-the-art caching systems often require the client to access the data through a caching proxy, or to contact a metadata server to locate the closest available copy of the desired data. Also, such caching systems are inconsistent with the design of distributed hashtable databases such as Dynamo, which focus on allowing clients to locate data independently. We argue there is a gap between existing state-of-the-art solutions and the needs of geographically distributed applications, which require fast access to popular objects while not degrading access latency for the rest of the data. In this paper, we introduce a probabilistic algorithm allowing the user to locate the closest copy of the data efficiently and independently with minimal overhead, allowing low-latency access to non-cached data. Also, we propose a network-efficient technique to identify the most popular data objects in the cluster and trigger their replication close to the clients. Experiments with a real-world data set show that these principles allow clients to locate the closest available copy of data with small memory footprint and low error-rate, thus improving read-latency for non-cached data and allowing hot data to be read locally.

INTRODUCTION

The rapid development of cloud computing solutions such as Amazon Web Services [2], Microsoft Azure [START_REF]Microsoft Azure[END_REF] or Google Cloud Platform [START_REF]Google Cloud[END_REF] has lead a number of applications to move from dedicated hardware to public clouds. This commonly allows the user to reserve and to use resources in multiple datacenter locations. Among these applications, scientific experiments tend to become geo-distributed as well. Such geo-distribution provides low-latency data access for scientists worldwide while minimizing bandwidth utilization and improving fault-tolerance as well as disaster-recovery. For instance, the MonALISA monitoring backend [START_REF] Legrand | MonALISA: An agent based, dynamic service system to monitor, control and optimize distributed systems[END_REF] of AL-ICE (A Large Ion Collider Experiment) [START_REF] Aamodt | The ALICE experiment at the CERN LHC[END_REF] is distributed over 300 sites around the world. This raises the question of data locality, i.e. the location of the data relative to the computational resources.

Content Distribution Networks [START_REF] Peng | CDN: content distribution network[END_REF], or CDNs, located between the end-user and the origin data location, help reducing content access latency by caching it as close as possible to the end-user. However, being loosely coupled with the underlying data storage, they are usually only suitable for content that changes rarely, such as media files or static resources. Different replication strategies are needed to provide the same data locality properties to frequently updated content such as database objects.

Replicating all content at each site is a simple way to increase read performance. However, it can lead to significantly degraded read performance due to the additional synchronization needed. It also results in poor network and storage resource usage by unnecessarily replicating data to sites from which it is only rarely queried. To alleviate from this, an option is to statically choose one (or several) sites for replication when the data is first saved into the system. Although resource overhead will be significantly lower, and assuming a pertinent site choice, this still ignores the fact that data popularity can change over time, potentially resulting in a non-optimal data locality.

Using real-time metrics to dynamically decide when to create (or remove) additional replicas while retaining strong data consistency can be challenging. Many previous attempts [START_REF] Sivakumar | D-tunes: self tuning datastores for geo-distributed interactive applications[END_REF] rely on one or multiple centralized metadata servers that indicate to the clients if and where a specific piece of data can be read locally. These central servers can potentially become hot spots under highly concurrent reads, being placed on the critical path of any read or write. The added synchronization needed between metadata servers will further decrease the write performance of the system. This centralized approach is hardly compatible with the largely distributed DHT-based design of storage systems such as Dynamo [START_REF] Decandia | Dynamo: Amazon's highly available key-value store[END_REF], Cassandra [START_REF] Lakshman | Cassandra: A decentralized structured storage system[END_REF] or Týr [START_REF] Matri | Týr: Efficient Transactional Storage for Data-Intensive Applications[END_REF], making poor usage of their communication protocols between nodes. Since all of these systems rely on a gossip-based [START_REF] Demers | Epidemic algorithms for replicated database maintenance[END_REF] failure detection algorithm and periodically exchange messages between nodes, we argue that it is desirable for caching algorithms to use these messages in order to disseminate live object usage information, consequently reducing the network overhead.

In this paper, we are introducing a set of algorithms and design principles that allow to dynamically replicate popular data objects as close as possible to the end-users and to consequently guarantee the lowest-possible request latency. Our approach leverages completely decentralized metadata in order to avoid traversing unnecessary network hops and to guarantee efficient location of all data objects. These features make our proposal particularly well-suited for looselycoupled DHT-and gossip-based storage systems.

To the best of our knowledge, we are proposing the first decentralized, dynamic replication system that allows the clients to locate the closest copy of the data independently of any metadata server. Our contributions are these:

• A probabilistic usage-aware dynamic replication module running on each node of the cluster, which efficiently collects and disseminates usage metrics (Section 2.1). It independently orders additional copies of the data to be created or deleted based on this information. The underlying algorithms rely on probabilistic principles to keep the metadata overhead and performance cost as small as possible.

• A probabilistic decentralized data-location algorithm embedded in the cluster nodes and in client libraries, allowing users to efficiently locate additional data copies without having to contact any centralized metadata server (Section 2.2). This ultimately allows for low-latency queries, avoiding hot spots in the cluster and improving its horizontal scalability.

• An experimental study of a proof-of-concept implementation leveraging the above principles and showing that the dynamic data replication set forth by our proposal does not come at the cost of reduced performance or increased request latency (Section 3).

We briefly discuss the obtained results (Section 4) and review the related work (Section 5). Finally, we conclude by outlining future work on our proposal (Section 6).

DESIGN OF A SELF-ADAPTIVE DATA STORAGE SYSTEM

In this paper, we assume one cluster whose nodes are located in different geographical locations, or sites. Each node can talk directly to any other on using the network, and all nodes are equal in terms of responsibilities in the cluster.

State-of-the-art distributed hash table based storage systems such as Dynamo introduce the concept of smart clients, i.e. clients that use cluster state information to route requests to the appropriate node directly. We propose to extend these capacities with dynamic replication awareness. Instead of using a metadata server to locate a piece of data it wants to read, a client is able to independently determine the closest available location of the data. Clients supporting this level of functionality are called smart clients. In order to preserve backward compatibility with other clients, referred to as naive clients, they will continue working properly but will not benefit from an optimal location of the data objects.

Collection of usage metrics

Preferred site indication

We assume that each client maintains a list of every site in the cluster ordered by preference (such as average network latency, geographical proximity or any other relevant metric). When submitting a read request to any node in the cluster, smart clients include an additional piece of information: the preferred site. The preferred site is the site the client would like the piece of data being read to be available in. This site is usually the closest site from the client, or the site with the lowest read latency.This preferred site is embedded in the request even if it has been sent to a node belonging to it.

Per-node popularity metrics collection

In order to be able to suggest additional copies of objects to be created in the cluster, each node keeps track of the most popular couples (object, pref erredsite) from the requests it receives from clients during a read.

The identification of these couples is performed by adapting a lightweight streaming algorithm designed to find the elements with the highest frequency in a stream of data: Space-Saving [START_REF] Metwally | Efficient computation of frequent and top-k elements in data streams[END_REF]. We recall here a few of its key features that will be leveraged later. Its probabilistic nature allows to identify the most frequent elements while keeping tight error guarantees and using a small and predictable amount of memory. The output of the algorithm is a list of candidate frequent elements, with multiple associated counters:

• ce is the counter structure for an element e,

• ce. ĥ the hit counter for this element, i.e. the total estimated number of hits of the element e,

• ce.ε is the error counter for the element e, i.e. the error margin of the hit counter for this element.

The basic Space-Saving algorithm is detailed in the pseudocode of Algorithm 1.

The Space-Saving algorithm identifies the k most frequent elements in a stream S. We are more interested in the most recent frequent items than in the historic frequent items. We then need to reinitialize the Space-Saving summary frequently, using only one configurable time window. In order to avoid starting each of these periods with an empty summary and because the Space-Saving structure has been demonstrated to be mergeable [START_REF] Agarwal | Mergeable summaries[END_REF], we merge the Space-Saving summary from the current time window t, Ct, with the one from the previous time window Ct-1. As such, the frequent element summary for the time window t, Ft, is given by Ft = MergeSpaceSaving(Ct, Ct-1).

Algorithm 1 Basic Space-Saving algorithm Input: k > 0 counters, element stream S Output: summary containing most frequent elements and their estimated frequency, C procedure SpaceSaving(S, k) C ← InitializeCounters(k) for all elements e in S do if e is monitored then let ce be the counter structure for e ce. ĥ ← ce. ĥ + 1 else let em be the element with least hits, ce m . ĥ replace em with e ce.ε ← ce m . ĥ ce. ĥ ← ce m . ĥ + 1 end if end for end procedure

Cluster-wide popularity metrics dissemination

Our algorithm builds on a gossip infection-style, weaklyconsistent protocol [START_REF] Demers | Epidemic algorithms for replicated database maintenance[END_REF][START_REF] Das | SWIM: scalable weakly-consistent infection-style process group membership protocol[END_REF] to disseminate object popularity information across the cluster with the lowest possible network overhead. Each node periodically sends to other randomly selected nodes in the cluster its own Space-Saving summary for the current time window, tagging it with the local time. In order to quickly relay the information clusterwide, the nodes piggyback these messages with recentlyreceived information from other nodes, including their associated time. Each node keeps a view of the most popular items for each node in the cluster, updating its state from another node only if the received one is more recent (based on its timestamp value). As such, for a n-node cluster and a configured value of k counters per summary, the size of the complete structure is predictable and in the order of O(kn).

Additional replica placement

Our proposal for replica placement relies on consistent hashing [START_REF] Karger | Consistent hashing and random trees: Distributed caching protocols for relieving hot spots on the world wide web[END_REF], as implemented by various decentralized storage systems such as Chord [START_REF] Stoica | Chord: A scalable peer-to-peer lookup service for internet applications[END_REF], Dynamo [START_REF] Decandia | Dynamo: Amazon's highly available key-value store[END_REF] or Týr [START_REF] Matri | Týr: Efficient Transactional Storage for Data-Intensive Applications[END_REF] for data placement. Given a hash function h(x), the output range [hmin, hmax] of the function is treated as a circular space (hmin sticking around to hmax). Every node is assigned a different random value within this range, which represents its position on the ring. The node responsible for an object obj with key k is determined by the result of h(k), giving a unique position h k on the ring. The first node encountered while walking the ring past this position is called the coordinator node of obj. The nodes to store additional replicas of an element are obtained by continuing walking the ring passed the coordinator node until an appropriate number of nodes are found. All these nodes are called natural storage nodes of obj.

Using the same ring, the node that will be chosen to hold an additional replica of the object obj in a site s, if any, is the first node in the site s encountered while continuing walking the ring past the last natural storage node for obj. This node being chosen deterministically, any additional replica for a given object on a specific site will always be stored on the same node.

Dynamic replica creation

Merging all the collected summaries received from every node gives each node an overview of the most frequent couples (object, site) cluster-wide. This merge is performed using the same Space-Saving merging function used in Section 2.1.2. Each node checks periodically whether such couples for objects it is the coordinator node for appear in this cluster-wide summary. If so, for each site at which an object obj is popular, a new replica is created at that site, if none exists already. The specific node on which such replica is created is determined as explained in Section 2. 1.4. Replica creations are advertised to all other nodes in the cluster using the gossip protocol defined in Section 2.1.3. Each node keeps a list of all objects for which additional replicas have been created, as well as the sites they have been replicated on.

Dynamic replica removal

Whenever the popularity of an object drops, additional replicas dynamically created are no longer needed and must be deleted. The node holding an additional replica for an object obj is responsible of independently deleting that copy in order to make room for other replicas. This deletion happens after an object previously present in the cluster-wide summary disappears from it.

In order to prevent replicas being repeatedly created and deleted for objects with popularity fluctuations near the Space-Saving threshold, this deletion may be delayed after a grace period greater than the summary window length defined in Section 2.1.2. Similarly to replica creation, this deletion is advertised using the same gossip protocol.

Accessing dynamic replicas

Decentralized data location

When a client reads an object, it tries to locate the closest copy of a desired data object within the cluster. To do this, it needs sufficient information in a memory-and networkefficient fashion. Bloom filters [START_REF] Bloom | Space/time trade-offs in hash coding with allowable errors[END_REF] provide an ideal base for this. Their compact nature makes them easy to forward over the network and lightweight to store in the clients memory, while providing fast lookup and tight guarantees on false positive rates.

Since replica creation or removal is advertised cluster-wide using gossip, each node is able to provide this information to any client. As such, we allow smart clients to address a request to any node on any site and receive in response a Bloom filter summarizing the current additionally replicated objects throughout the cluster, identified by the couple (object, site).

To read an object obj, a client needs the Bloom filter obtained from any node of the cluster as well the ordered list S of preferred sites. The site s to send a request to is either the first one in the S list or the one for which the Bloom filter search for (obj, s) returns a positive result. The algorithm is detailed in the pseudocode of Algorithm 2.

At the chosen site s, the request is addressed to the node that would hold a dynamic replica, using the same deterministic algorithm described in Section 2.1.4.

Server read request handling

When a node receives a request for an object it holds, it responds to the client normally. In some cases however, our algorithm can cause the client to contact a node not holding the desired data. This can happen either because the client used an outdated Bloom filter (i.e. the filter returned a false positive because the cluster had not yet converged on the current replication status when the Bloom filter was obtained) or because the client was of naive type. Should the client contact the wrong node in such case, that node will directly forward the request to the closest node holding the data according to its knowledge of the replicated objects in the cluster, or using the DHT information.

Relaying the request instead of responding to the client incurs fewer messages, making a false guess from the client result in only a single message overhead compared to a successful guess, furthermore addressed to a close site.

Bloom filter refresh

As replicated objects may change depending on the performed user requests, the Bloom filters need to be periodically refreshed to keep the false positive rate of our location algorithm as low as possible. As such, after a configurable period of time, clients will contact a random node of the cluster to request an up-to-date version of the summary.

To cope with cases where the configured summary validity period is too large for some usage patterns, the client keeps a count of false positives or false negatives. Should the client address the request to a server not holding the desired piece of data, the server relays the query as usual to the correct node but flags the answer as erroneous. This may happen for instance because it has been deleted after the client bloom filter has been last updated. We call this case a false positive. Inversely, if the client accesses a remote object replica and a closer one existed, which we call a false negative, the server responding to the client flags the response message as well. A counter of wrong assumptions is maintained by each client, and set to 0 every time the summary is updated. For each false positive or negative, the client increments its counter. It updates its bloom filter as soon as the counter reaches a configurable threshold.

Additional optimizations

Since the Bloom filter summary of the additional replicas can be requested from any node in the cluster, clients can piggyback update requests to any request to the cluster. In turn, nodes are able to piggyback the latest summary to the response. This results in a sensible reduction of the network utilization.

Additionally, a desirable behavior of the location algorithm is to allow nodes to use newly-created replicas of popular data as soon as possible after they have been created, without waiting for a summary refresh by clients. As such, whenever a cluster node receives a read request from a client, it can check if another copy of the data is available closer to the client. In that case, the node piggybacks this information to the response, along with an updated summary that will allow the client to read from the closest available copy in the future.

EVALUATION

Experimental setup

We evaluated our design using a synthetic simulation on 80 processes running on a single machine, replicating the behavior of multiple nodes over multiple sites -5 sites composed of 16 machines each. Gossip between processes is based on the SWIM protocol [START_REF] Das | SWIM: scalable weakly-consistent infection-style process group membership protocol[END_REF], implemented by Serf [START_REF] Serf | [END_REF]. Latency between nodes is simulated using the native Linux Traffic Control [START_REF]Linux Traffic Control[END_REF] utility. Latency settings have been set to realistic values measured by Verizon [START_REF]Verizon Enterprise Solutions IP Latency Statistics[END_REF] on their own network, and are detailed in Table 1. The simulation itself was implemented using approximately 700 lines of Python code, not counting existing open-source implementations of Bloom filters and Space-Saving algorithm.

Object requests were simulated using an open data set composed of real usage data collected at the University of Indiana [START_REF] Meiss | Ranking web sites with real user traffic[END_REF]. It is composed of the URLs of more than 25 billion user requests collected over two years. Usual web access data have been demonstrated to follow a Zipfian distribution [START_REF] Powers | Applications and explanations of zipf's law[END_REF][START_REF] Adamic | Zipf's law and the internet[END_REF], causing hot spots on storage clusters. In our experiments, we consider URLs to be a data object and simulate read requests originating from all 5 different, geographically disperse sites.

For our experiments, the Space-Saving time window length defined in Section 2.1.2 was set to 30 seconds, the number of counters to 128, and the removal grace period defined in 0 100 200 300 400 500 600 700 800 900 5 Section 2.1.6 to 90 seconds. This experiment was ran for a total duration of 15 minutes.

Object replication

We argue that our technique is able to efficiently identify the most popular data objects throughout the cluster and replicate them as close as possible to the interested clients. Improper identification of these popular objects would result in no significant decrease in the average read latency. In that case, non-frequent data objects would be replicated, and these would be unlikely to cause any noticeable difference in the obtained results.

In our results, plotted in Figure 1, we can observe a clear decrease of the mean request latency with dynamic replication enabled compared to the results without replication. We notice at the beginning on the curve the initial period of time in which the popular data objects are being identified and progressively replicated. High variability in the results is caused by the inner characteristics of the data set, where frequent items tend sometimes to be grouped With our dynamic replication technique, after result stabilization, the achieved average read latency dropped by 38%, from 16 ms to 10 ms. It is also interesting to observe that that the frequent elements were identified quickly, in under three times the Space-Saving window length. The number of replicated objects stabilized at an average of 147 replicated objects in the cluster, slightly more than the configured number of counters in the Space-Saving algorithm, due to the grace period before dynamically created replicas of previously-popular objects are deleted.

Object location

Our client-side algorithm is able to locate the additional replicas in the cluster with minimal error rate, thanks to the fast convergence of the gossip information dissemination and to the guaranteed error rate of the Bloom filter summary sent to the clients. We want to evaluate the accuracy of our proposal, by measuring and plotting the error rates of this algorithm. We separate the false positives -when a node is wrongly believed to hold a replica of the object, from the false negatives -when the request is sent to a non-optimal replica of the object, i.e. when a closer replica existed. A false positive indicates that a client requested an object replica at a site where it was not available. This can happen if the bloom filter of the client is outdated and the replica has been deleted. A false negative indicates that the client has requested an origin copy of the data while a local replica existed. This may occur if the bloom filter of the client was last updated before this replica was created.

The results are depicted in Figure 2. We observe that our algorithm correctly identifies the most-optimal replica of the desired data in more than 98% of the cases, independently and without requiring the assistance of a metadata server.

The initial high false negative rate is caused by the first replicas being created, and is explained by the short period of time before the clients are made aware of these copies. The error counter described in Section 2.2.3 helps to quickly reduce this error rate.

Not surprisingly, the false positive rate increases around 120 seconds after the start of the experiment. This is caused by the first additional replicas being deleted shortly after the configured 90-second grace period plus 30 seconds Bloom filter window length.

DISCUSSION

Impact of the workload type

Our work targets applications in which reads dominate over writes. For applications with a lower read-to-write ratio, the potential gains of dynamic replication might not be as clear. Specifically, should an object be updated, dynamic replicas need to be kept synchronized with the original data. As the number of dynamic replicas of any given object grows, the overhead of this synchronization grows as well. This will likely result in a decreased write performance for this object, which might hinder the performance of the application.

Future work will focus on experimenting our algorithm on top of a production system to evaluate the actual impact of our proposal on the performance of a real-world application with different workloads.

Decentralized message dissemination

Our proposal relies heavily on a gossip protocol being used by the underlying storage system. For storage systems not based on a gossip protocol, such as HDFS [START_REF] Shafer | The Hadoop distributed filesystem: Balancing portability and performance[END_REF], the benefits of adding such a communication layer instead of relying on the existing centralized architecture still has to be evaluated.

RELATED WORK

Self-adaptive storage. These systems emerged in the context of autonomic computing, a paradigm inspired by biology (i.e. by how the human nervous system reacts to external changes in an autonomous manner, through unconscious reflexes in order to adapt our body to its needs and to the environment without requiring our attention). IBM has introduced a reference model for autonomic control loops [START_REF]An architectural blueprint for autonomic computing[END_REF], based on monitoring, analysis and reactions, currently used by most autonomic storage systems. In [START_REF] Carpen-Amarie | Bringing Introspection into BlobSeer: Towards a Self-Adaptive Distributed Data Management System[END_REF] Carpen et al. make a first step towards enhancing with self-adaptive support a massively distributed storage system. They add an introspection layer that collects data about the usage of storage resources and data access patterns relying on the MonALISA [START_REF] Legrand | MonALISA: An agent based, dynamic service system to monitor, control and optimize distributed systems[END_REF] monitoring framework. However, in contrast to our solution, this approach lacks the reaction phase that closes the control loop by adjusting the storage system according to the tracked data. More recent solutions turn to adapting some specific parameters, e.g. the replication factor [START_REF] Knežević | DHT-Based Self-adapting Replication Protocol for Achieving High Data Availability[END_REF]. Here, a DHT-based self-adapting replication protocol is introduced to determine the locations of replicas and then to autonomously adjust their number to deliver a configured data availability guarantee. However, these works focus on the relationship between the number of replicas and performance, completely ignoring the impact of replica location on query efficiency. D-Tunes [START_REF] Sivakumar | D-tunes: self tuning datastores for geo-distributed interactive applications[END_REF], the related work most similar to ours, tries to fill this gap. It is based on some selftuning algorithms that can adapt to workload changes over short time-scales by automatically configuring parameters like replication factors, consistency levels and readjusting read/write priorities while judiciously recommending data placements over longer time-scales. However, D-Tunes does change replication settings on a global-scale, making it difficult to perform fine-grained replication optimization on a per-object basis.

Dynamic replication. A vast number of works acknowledge that static, manual replication introduces a serious overhead and can drastically affect the storage performance. This is the case with geo-replicated datastores such as Spanner [START_REF] Corbett | Spanner: Google's globally distributed database[END_REF] and Cassandra [START_REF] Lakshman | Cassandra: A decentralized structured storage system[END_REF], in which replication strategies are manually configured by the application developers. An alternative consists of dynamically distributing workloads by replicating and migrating replicas among storage nodes. Since dynamic multiple-location replication is NP-complete in nature, most of the existing strategies for dynamic replication are typically based on so-called single-location algorithms for identifying a single site for data replication. Dong et al. [START_REF] Dong | On dynamic replication strategies in data service grids[END_REF] transform the multiple-location problem into several classical mathematical problems with different parameter settings, for which efficient approximation algorithms exist. However, they don't consider the impact of replication granularity on performance and scalability. Wei et al. [START_REF] Wei | Dynamic replication management for object-based storage system[END_REF] address this issue by trying to answer the question: how many replicas the system should keep at least to satisfy availability requirements? To this end, a model is developed to express availability as a function of replica number. This approach, however, only works within a single site, as it assumes uniform bandwidth and latency, which is not the case with the geo-distributed workloads that we target. Inspired by the P2P systems, [START_REF] Shen | Ead: An efficient and adaptive decentralized file replication algorithm in p2p file sharing systems[END_REF] proposes an adaptive decentralized file replication algorithm that achieves high query efficiency and high replica utilization at a significantly low cost. The idea is to select query traffic hubs and frequent requesters as replica nodes, and to dynamically adapt to nonuniform and time-varying file popularity and node interest. Unlike current methods and similarly to our approach, they create and delete replicas in a decentralized self-adaptive manner guaranteeing high replica utilisation. While this solution works gracefully in a P2P environment, its assumption of a known underlying network topology (which does not hold in a cloud environment) and the use of additional intermediate nodes (hubs) along the path between the client and the servers make it unusable in modern storage systems.

CONCLUSIONS

In this paper, introduce techniques and algorithms designed to help reduce read latency in a geographically distributed storage cluster. To do so, we propose a method to collect live request metrics using probabilistic algorithms in order to identify the most popular data objects. These metrics are disseminated throughout the cluster using gossiping. We further propose another probabilistic, decentralized algorithm for efficiently locating the closest available copy of the data. It does not require the use of a centralized metadata center, thus improving the fault-tolerance and horizontalscalability of the cluster, while at the same time reducing the client read latency.

Our experiments over a large real-world data set show that these techniques can efficiently identify and duplicate the most frequently requested data objects in a cluster. Additionally, they allow the client to independently and probabilistically locate the additional copies of the data with a low error rate.

In order to demonstrate the performance of our approach, future work will focus on implementing these techniques on top of a real storage system such as Cassandra. We also plan to evaluate its performance with additional real-world use cases and data sets. Finally, we are investigating the use of machine-learning techniques to be able to identify with higher accuracy the most important data objects to be replicated.

Figure 1 :

 1 Figure 1: Achieved read latency with and without dynamic replication enabled.

Figure 2 :

 2 Figure 2: Measured false positive and false negative rates.

Table 1 :

 1 Simulated round-trip latencies between sites (ms)

		US East US West Europe Asia Pacific
	US East	0.25	35	70	145	185
	US West 35	0.25	105	110	150
	Europe	70	105	0.25	120	310
	Asia	145	110	120	0.25	140
	Pacific	185	150	310	140	0.25
	Algorithm 2 Object location algorithm
	Input: object to read obj, list S of all sites ordered by pref-
	erence.
	ss ← ∅
	function LocateClosestReplica(obj)
	if ss = ∅ then
	let sr be a random node in the cluster
	ss = RequestBloomFilter(sr)
	end if
	let N be the list of natural sites for obj
	for all site s in S do
	if s in N or (obj, s) in ss then
	return s
	end if
	end for
	end function

ACKNOWLEDGMENTS

This work is part of the BigStorage project, funded by the European Union under the Marie Sk lodowska-Curie Actions (H2020-MSCA-ITN-2014-642963). The authors would like to thank the reviewers who provided valuable counsel and expertise that greatly helped this research.


Paul Escande
email: paul.escande@gmail.com.

Pierre Weiss
email: pierre.armand.weiss@gmail.com

Approximation of integral operators using product-convolution expansions

Keywords: Integral operators, wavelet, spline, structured low rank decomposition, numerical complexity, approximation rate, fast Fourier transform AMS classifications: 47A58, 41A46, 41A35, 65D07, 65T60, 65T50, 65R32, 94A12

HAL is

Introduction

We are interested in the compact representation and fast evaluation of a class of space or time varying linear integral operators with regular variations. Such operators appear in a large number of applications ranging from wireless communications [START_REF] Chul | An efficient feedback method for MIMO systems with slowly time-varying channels[END_REF][START_REF] Hrycak | Low complexity equalization for doubly selective channels modeled by a basis expansion[END_REF] to seismic data analysis [START_REF] Griffiths | Adaptive deconvolution: A new technique for processing time-varying seismic data[END_REF], biology [START_REF] Gilad | A fast algorithm for convolution integrals with space and time variant kernels[END_REF] and image processing [START_REF] Sawchuk | Space-variant image motion degradation and restoration[END_REF].

In all these applications, a key numerical problem is to efficiently evaluate the action of the operator and its adjoint on given functions. This is necessary -for instance -to design fast inverse problems solvers. The main objective of this paper is to analyze the complexity of a set of approximation techniques coined product-convolution expansions.

We are interested in bounded linear integral operator H : L 2 (Ω) → L 2 (Ω) defined from a kernel K by: Hu(x) = Ω K(x, y)u(y) dy.

(1.1) for all u ∈ L 2 (Ω), where Ω ⊂ R d . Evaluating integrals of type (1.1) is a major challenge in numerical analysis and many methods have been developed in the literature. Nearly all methods share the same basic principle: decompose the operator kernel as a sum of low rank matrices with a multi-scale structure. This is the case in panel clustering methods [START_REF] Hackbusch | On the fast matrix multiplication in the boundary element method by panel clustering[END_REF], hierarchical matrices [START_REF] Börm | Introduction to hierarchical matrices with applications[END_REF], cross approximations [START_REF] Oseledets | TT-cross approximation for multidimensional arrays[END_REF] or wavelet expansions [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF]. The method proposed in this paper basically shares the same idea, except that the time varying impulse response T of the operator is decomposed instead of the kernel K. The time varying impulse response (TVIR) T of H is defined by: T (x, y) = K(x + y, y).

(1.

2)

The TVIR representation of H allows formalizing the notion of regularly varying integral operator: the functions T (x, •) should be "smooth" for all x ∈ Ω. Intuitively, the smoothness assumption means that two neighboring impulse responses should only differ slightly. Under this assumption, it is tempting to approximate H locally by a convolution. Two different approaches have been proposed in the literature to achieve this. The first one is called convolution-product expansion of order m and consists of approximating H by an operator H m of type:

H m u = m k=1 w k (h k u), (1.3)
where denotes the standard multiplication for functions and the Hadamard product for vectors, and denotes the convolution operator. The second one, called product-convolution expansion of order m, is at the core of this paper and consists of using an expansion of type:

H m u = m k=1 h k (w k u). (1.4)
These two types of approximations have been used for a long time in the field of imaging (and to a lesser extent mobile communications and biology) and progressively became more and more refined [START_REF] Trussell | Identification and restoration of spatially variant motion blurs in sequential images[END_REF][START_REF] Nagy | Restoring images degraded by spatially variant blur[END_REF][START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF][START_REF] Gilad | A fast algorithm for convolution integrals with space and time variant kernels[END_REF][START_REF] Hrycak | Low complexity equalization for doubly selective channels modeled by a basis expansion[END_REF][START_REF] Muthuvel Arigovindan | A parallel product-convolution approach for representing the depth varying point spread functions in 3d widefield microscopy based on principal component analysis[END_REF][START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF][START_REF] Miraut | Efficient shift-variant image restoration using deformable filtering (Part I)[END_REF][START_REF] Denis | Fast approximations of shift-variant blur[END_REF]. In particular, the recent work [START_REF] Denis | Fast approximations of shift-variant blur[END_REF] provides a nice overview of existing choices for the functions h k and w k as well as new ideas leading to significant improvements. Many different names have been used in the literature to describe expansions of type (1.3) and (1.4) depending on the communities: sectional methods, overlap-add and overlap-save methods, piecewise convolutions, anisoplanatic convolutions, parallel productconvolution, filter flow, windowed-convolutions... The term product-convolution comes from the field of mathematics [START_REF] Busby | Product-convolution operators and mixed-norm spaces[END_REF] 1 . We believe that it precisely describes the set of expansions of type (1.3) and therefore chose this naming. It was already used in the field of imaging by [START_REF] Muthuvel Arigovindan | A parallel product-convolution approach for representing the depth varying point spread functions in 3d widefield microscopy based on principal component analysis[END_REF]. Now that product-convolution expansions have been described, natural questions arise:

i) How to choose the functions h k and w k ?

ii) What is the numerical complexity of evaluating products of type H m u?

iii) What is the resulting approximation error H m -H , where • is a norm over the space of operators? iv) How many operations are needed in order to obtain an approximation H m such that H m -H ≤ ?

Elements i) and ii) have been studied thoroughly and improved over the years in the mentioned papers. The main questions addressed herein are points iii) and iv). To the best of our knowledge, they have been ignored until now. They are however necessary in order to evaluate the theoretical performance of different product-convolution expansions and to compare their respective advantages precisely.

The main outcome of this paper is the following: under smoothness assumptions of type T (x, •) ∈ H s (Ω) for all x ∈ Ω (the Hilbert space of functions in L 2 (Ω) with s derivatives in L 2 (Ω)), most methods proposed in the literature -if implemented correctly -ensure a decay of type H m -H HS = O(m -s ), where • HS is the Hilbert-Schmidt norm. Moreover, this bound cannot be improved uniformly on the considered smoothness class. By adding a support condition of type supp(T (x, •)) ⊆ [-κ/2, κ/2], the bound becomes H m -H HS = O( √ κm -s ). More importantly, bounded supports allow reducing the computational burden. After discretization on n time points, we show that the number of operations required to satisfy H m -H HS ≤ vary from O κ

1 2s n log 2 (n) -1/s to O κ 2s+1 2s n log 2 (κn) -1/s
depending on the choices of w k and h k . We also show that the compressed operator representations of Meyer [START_REF] Meyer | Wavelets and operators[END_REF] can be used under additional regularity assumptions.

An important difference of product-convolution expansions compared to most methods in the literature [START_REF] Hackbusch | On the fast matrix multiplication in the boundary element method by panel clustering[END_REF][START_REF] Börm | Introduction to hierarchical matrices with applications[END_REF][START_REF] Oseledets | TT-cross approximation for multidimensional arrays[END_REF][START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF][START_REF] Escande | Sparse wavelet representations of spatially varying blurring operators[END_REF] is that they are insensitive to the smoothness of T (•, y). The smoothness in the x direction is a useful property to control the discretization error, but not the approximation rate. The proposed methodology might therefore be particularly competitive in applications with irregular impulse responses.

The paper is organized as follows. In section 2, we describe the notation and introduce a few standard results of approximation theory. In section 3, we precisely describe the class of operators studied in this paper, show how to discretize them and provide the numerical complexity of evaluating product-convolution expansions of type (1.4). Sections 4 and 5 contain the full approximation analysis for two different kinds of approaches called linear or adaptive methods. Section 6 contains a summary and a few additional comments.

Notation

Let a and b denote functions depending on some parameters. The relationship a b means that a and b are equivalent, i.e. that there exists 0 < c 1 ≤ c 2 such that c 1 a ≤ b ≤ c 2 a. Constants appearing in inequalities will be denoted by C and may vary at each occurrence. If a dependence on a parameter exists (e.g. ), we will use the notation C( ).

In most of the paper, we work on the unit circle Ω = R\Z sometimes identified with the interval -1 2 , 1 2 . This choice is driven by simplicity of exposition and the results can be extended to bounded domains such as Ω = [0, 1] d (see section 6.2). Let L 2 (Ω) denote the space of square integrable functions on Ω. The Sobolev space H s (Ω) is defined as the set of functions in L 2 (Ω) with weak derivatives up to order s in L 2 (Ω). The k-th weak derivative of u ∈ H s (Ω) is denoted u (s) . The norm and semi-norm of u ∈ H s (Ω) are defined by:

u H s (Ω) = s k=0 u (k) L 2 (Ω) and |u| H s (Ω) = u (s) L 2 (Ω) . (2.1)
The sequence of functions (e k ) k∈Z where e k : x → exp(-2iπkx) is a Hilbert basis of L 2 (Ω) (see e.g. [START_REF] Katznelson | An introduction to harmonic analysis[END_REF]).

Definition 2.1. Let u ∈ L 2 (Ω) and e k : x → exp(-2iπkx) denote the k-th Fourier atom. The Fourier series coefficients û[k] of u are defined for all k ∈ Z by:

û[k] = Ω u(x)e k (x) dx. (2.2)
The space H s (Ω) can be characterized through Fourier series.

Lemma 2.1 (Fourier characterization of Sobolev norms).

u 2 H s (Ω) k∈Z |û[k]| 2 (1 + |k| 2 ) s . (2.3)
Definition 2.2 (B-spline of order α). Let α ∈ N and m ≥ α + 2 be two integers. The B-spline of order 0 is defined by

B 0,m = 1 [-1/(2m),1/(2m)] . (2.4)
The B-spline of order α ∈ N * is defined by recurrence by:

B α,m = mB 0,m B α-1,m = m α B 0,m . . . B 0,m α times . (2.5)
The set of cardinal B-splines of order α is denoted B α,m and defined by:

B α,m = f (•) = m-1 k=0 c k B α,m (• -k/m), c k ∈ R, 0 ≤ k ≤ m -1 . (2.6)
In this work, we use the Daubechies wavelet basis for L 2 (R) [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF]. We let φ and ψ denote the scaling and mother wavelets and assume that the mother wavelet ψ has α vanishing moments, i.e.

∀0 ≤ m < α, [0,1] t m ψ(t)dt = 0.
(2.7)

Daubechies wavelets satisfy supp(ψ) = [-α + 1, α], see [28, Theorem 7.9, p. 294]. Translated and dilated versions of the wavelets are defined, for all j > 0 by

ψ j,l (x) = 2 j/2 ψ 2 j x -l . (2.8)
The set of functions (ψ j,l ) j∈N,l∈Z , is an orthonormal basis of L 2 (R) with the convention ψ 0,l = φ(x -l). There are different ways to construct a wavelet basis on the interval [-1/2, 1/2] from a wavelet basis on L 2 (R). Here, we use boundary wavelets defined in [START_REF] Cohen | Wavelets on the interval and fast wavelet transforms[END_REF]. We refer to [START_REF] Daubechies | Ten lectures on wavelets[END_REF][START_REF] Mallat | A wavelet tour of signal processing[END_REF] for more details on the construction of wavelet bases. This yields an orthonormal basis (ψ λ ) λ∈Λ of L 2 (Ω), where

Λ = (j, l), j ∈ N, 0 ≤ l ≤ 2 j . (2.9)
We let I λ = supp(ψ λ ) and for λ ∈ Λ, we use the notation |λ| = j. Let u and v be two functions in L 2 (Ω), the notation u ⊗ v will be used both to indicate the function w ∈ L 2 (Ω × Ω) defined by

w(x, y) = (u ⊗ v)(x, y) = u(x)v(y), (2.10)
or the Hilbert-Schmidt operator w : L 2 (Ω) → L 2 (Ω) defined for all f ∈ L 2 (Ω) by:

w(f ) = (u ⊗ v)f = u, f v. (2.11)
The meaning can be inferred depending on the context. Let H : L 2 (Ω) → L 2 (Ω) denote a linear integral operators. Its kernel will always be denoted K and its time varying impulse response T . The linear integral operator with kernel T will be denoted J.

The following result is an extension of the singular value decomposition to operators.

H = k≥1 σ k • e k ⊗ f k .
(2.12)

A function u ∈ L 2 (Ω) is denoted in regular font whereas its discretized version u ∈ R n is denoted in bold font. The value of function u at x ∈ Ω is denoted u(x), while the i-th coefficient of vector u ∈ R n is denoted u[i]. Similarly, an operator H : L 2 (Ω) → L 2 (Ω) is denoted in upper-case regular font whereas its discretized version H ∈ R n×n is denoted in upper-case bold font.

Preliminary facts

In this section, we gather a few basic results necessary to derive approximation results.

Assumptions on the operator and examples

All the results stated in this paper rely on the assumption that the TVIR T of H is a sufficiently simple function. By simple, we mean that i) the functions T (x, •) are smooth for all x ∈ Ω and ii) the impulse responses T (•, y) have a bounded support or a fast decay for all y ∈ Ω.

There are numerous ways to capture the regularity of a function. In this paper, we assume that T (x, •) lives in the Hilbert spaces H s (Ω) for all x ∈ Ω. This hypothesis is deliberately simple to clarify the proofs and the main ideas. Definition 3.1 (Class T s ). We let T s denote the class of functions T : Ω × Ω → R satisfying the smoothness condition:

T (x, •) ∈ H s (Ω), ∀x ∈ Ω and T (x, •) H s (Ω) is uniformly bounded in x, i.e: sup x∈Ω T (x, •) H s (Ω) ≤ C < +∞. (3.1)
Note that if T ∈ T s , then H is a Hilbert-Schmidt operator since:

H 2 HS = Ω Ω K(x, y) 2 dx dy (3.2)
= Ω Ω T (x, y) 2 dx dy

(3.3) = Ω T (x, •) 2 L 2 (Ω) dx < +∞. (3.4)
We will often use the following regularity assumption.

Assumption 3.1. The TVIR T of H belongs to T s .

In many applications, the impulse responses have a bounded support, or at least a fast spatial decay allowing to neglect the tails. This property will be exploited to design faster algorithms. This hypothesis can be expressed by the following assumption.

Assumption 3.2. T (x, y) = 0, ∀|x| > κ/2.

Examples

We provide 3 examples of kernels that may appear in applications. Figure 1 shows each kernel as a 2D image, the associated TVIR and the spectrum of the operator J (the linear integral operator with kernel T ) computed with an SVD.

Example 3.1. A typical kernel that motivates our study is defined by:

K(x, y) = 1 √ 2πσ(y) exp - (x -y) 2 2σ 2 (y) . (3.5)
The impulse responses K(•, y) are Gaussian for all y ∈ Ω. Their variance σ(y) > 0 varies depending on the position y. The TVIR of K is defined by:

T (x, y) = 1 √ 2πσ(y) exp - x 2 2σ 2 (y) . (3.6)
The impulse responses T (•, y) are not compactly supported, therefore, κ = 1 in Assumption 3.2. However, it is possible to truncate them by setting κ = 3 sup y∈Ω σ(y) for instance. This kernel satisfies Assumption 3.1 only if σ : Ω → R is sufficiently smooth. In Figure 1, left column, we set σ(y) = 0.08 + 0.02 cos(2πy).

Example 3.2. The second example is given by:

T (x, y) = 2 σ(y) max(1 -2σ(y)|x|, 0). (3.7)
The impulse responses T (•, y) are cardinal B-splines of degree 1 and width σ(y) > 0. They are compactly supported with κ = sup y∈Ω σ(y). This kernel satisfies Assumption 3.2 only if σ : Ω → R is sufficiently smooth. In Figure 1, central column, we set σ(y

) = 0.1 + 0.3(1 -|y|).
This kernel satisfies Assumption 3.1 with s = 1.

Example 3.3. The last example is a discontinuous TVIR. We set:

T (x, y) = g σ 1 (x)1 [-1/4,1/4] (y) + g σ 2 (x)(1 -1 [-1/4,1/4] (y)), (3.8)
where

g σ (x) = 1 √ 2π exp -x 2 σ 2 .
This corresponds to the last column in Figure 1, with σ 1 = 0.05 and σ 2 = 0.1. For this kernel, both Assumptions 3.1 and 3.2 are violated. Notice however that T is the sum of two tensor products and can therefore be represented using only four 1D functions. The spectrum of J should have only 2 non zero elements. This is verified in Figure 1i up to numerical errors.

Product-convolution expansions as low-rank approximations

Though similar in spirit, convolution-product (1.3) and product-convolution (1.4) expansions have a quite different interpretation captured by the following lemma.

Lemma 3.1. The TVIR of the convolution-product expansion T m in (1.3) is given by:

T m (x, y) = m k=1 h k (x)w k (x + y).
(3.9)

The TVIR of the product-convolution expansion T m in (1.4) is given by: Proof. We only prove (3.10) since the proof of (3.9) relies on the same arguments. By definition:

T m (x, y) = m k=1 h k (x)w k (y). ( 3
(H m u)(x) = m k=1 h k (w k u) (x) (3.11) = Ω m k=1 h k (x -y)w k (y)u(y) dy. (3.12)
By identification, this yields:

K m (x, y) = m k=1 h k (x -y)w k (y), (3.13)
so that

T m (x, y) = m k=1 h k (x)w k (y). (3.14)
As can be seen in (3.10), product-convolution expansions consist of finding low-rank approximations of the TVIR. This interpretation was already proposed in [START_REF] Denis | Fast approximations of shift-variant blur[END_REF] for instance and is the key observation to derive the forthcoming results. The expansion (3.9) does not share this simple interpretation and we do not investigate it further in this paper.

Discretization

In order to implement a product-convolution expansion of type 1.4, the problem first needs to be discretized. We address this problem with a Galerkin formalism. Let (ϕ 1 , . . . , ϕ n ) be a basis of a finite dimensional vector space V n of L 2 (Ω). Given an operator H : L 2 (Ω) → L 2 (Ω), we can construct a matrix H n ∈ R n×n defined for all 1 ≤ i, j ≤ n by H n [i, j] = Hϕ j , ϕ i . Let S n : H → H n denote the discretization operator. From a matrix H n , an operator H n can be reconstructed using, for instance, the pseudo-inverse S n,+ of S n . We let H n = S n,+ (H n ). For instance, if (ϕ 1 , . . . , ϕ n ) is an orthonormal basis of V n , the operator H n is given by:

H n = S n,+ (H n ) = 1≤i,j≤n H n [i, j]ϕ i ⊗ ϕ j . (3.15)
This paper is dedicated to analyzing methods denoted A m that provide an approximation H m = A m (H) of type (1.4), given an input operator H. Our analysis provides guarantees on the distance H -H m HS depending on m and the regularity properties of the input operator H, for different methods. Depending on the context, two different approaches can be used to implement A m .

• Compute the matrix H n m = S n (H m ) using numerical integration procedures. Then create an operator H n m = S n,+ (H n m ). This approach suffers from two defects. First, it is only possible by assuming that the kernel of H is given analytically. Moreover it might be computationally intractable. It is illustrated below.

H - A m H m - S n H n m - S n,+ H n m
• In many applications, the operator H is not given explicitly. Instead, we only have access to its discretization H n . Then it is possible to construct a discrete approximation algorithm A m yielding a discrete approximation H n m = A m (H n ). This matrix can then be mapped back to the continuous world using the pseudo-inverse:

H n m = S n,+ (H n m )
. This is illustrated below. In this paper, we will analyze the construction complexity of H n m using this second approach.

H - S n H n - A m H n m - S n,+ H n m
Ideally, we would like to provide guarantees on H -H n m HS depending on m and n. In the first approach, this is possible by using the following inequality:

H -H n m HS ≤ H -H m HS a(m) + H m -H n m HS d (n) , (3.16)
where a (m) is the approximation error studied in this paper and d (n) is the discretization error. Under mild regularity assumptions on K, it is possible to obtain results of type d (n) = O(n -γ ), where γ depends on the smoothness of K. For instance, if K ∈ H r (Ω × Ω), the error satisfies d (n) = O(n -r/2 ) for many bases including Fourier, wavelets and B-splines [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. For K ∈ BV (Ω × Ω), the space of functions with bounded variations, d (n) = O(n -1/4 ), see [START_REF] Mallat | A wavelet tour of signal processing[END_REF]Theorem 9.3]. As will be seen later, the approximation error a (m) behaves like O(m -s ). The proposed approximation technique is valuable only in the case m n, so that the discretization error can be considered negligible compared to the approximation error most practical applications. In all the paper, we assume that d (n) is negligible compared to a (m) without mention.

In the second approach, the error analysis is more complex since there is an additional bias due to the algorithm discretization. This bias is captured by the following inequality:

H -H n m HS ≤ H -H n HS d (n) + H n -A m (H n ) HS a(m) + A m (H n ) -H n m HS b (m,n) . (3.17) The bias b (m, n) = A m (S n,+ (S n (H))) -S n,+ (A m (S n (H)))
HS accounts for the difference between using the discrete or continuous approximation algorithm. In this paper, we do not study this bias error and assume that it is negligible compared to the approximation error a .

Implementation and complexity

Let F n ∈ C n×n denote the discrete inverse Fourier transform and F * n denote the discrete Fourier transform. Matrix-vector products F n u or F * n u can be evaluated in O(n log 2 (n)) operations using the fast Fourier transform (FFT). The discrete convolution product v = h u is defined for all i ∈ Z by v

[i] = n j=1 u[i -j]h[j],
with circular boundary conditions.

Discrete convolution products can be evaluated in O(n log 2 (n)) operations by using the following fundamental identity:

v = F n • ((F * n h) (F * n u)). (3.18)
Hence a convolution can be implemented using three FFTs (O(n log 2 (n)) operations) and a point-wise multiplication (O(n) operations). This being said, it is straightforward to implement formula (1.4) with an O(mn log 2 (n)) algorithm.

Under the additional assumption that w k and h k are supported on bounded intervals, the complexity can be improved. We assume that, after discretization, h k and w k are compactly supported, with support length q k ≤ n and p k ≤ n respectively. Lemma 3.2. A matrix-vector product of type (1.4) can be implemented with a complexity that does not exceed O ( m k=1 (p k + q k ) log 2 (min(p k , q k ))) operations. Proof. A convolution product of type h k (w k u) can be evaluated in O((p k + q k ) log(p k + q k )) operations. Indeed, the support of h k (w k u) has no more than p k + q k contiguous nonzeros elements. Using the Stockham sectioning algorithm [START_REF] Thomas | High-speed convolution and correlation[END_REF], the complexity can be further decreased to O((p k + q k ) log 2 (min(p k , q k ))) operations. This idea was proposed in [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF].

Projections on linear subspaces

We now turn to the problem of choosing the functions h k and w k in equation (1.4). The idea studied in this section is to fix a subspace E m = span(e k , k ∈ {1, . . . , m}) of L 2 (Ω) and to approximate T (x, •) as:

T m (x, y) = m k=1 c k (x)e k (y). (4.1)
For instance, the coefficients c k can be chosen so that

T m (x, •) is a projection of T (x, •) onto E m .
We propose to analyze three different different family of functions e k : Fourier atoms, wavelets atoms and B-splines. We analyze their complexity and approximation properties as well as their respective advantages.

Fourier decompositions

It is well known that functions in H s (Ω) can be well approximated by linear combination of low-frequency Fourier atoms. This loose statement is captured by the following lemma.

Lemma 4.1 ([16, 15]). Let f ∈ H s (Ω) and f m denote its partial Fourier series:

f m = m k=-m f [k]e k , (4.2)
where e k (y) = exp(-2iπky). Then

f m -f L 2 (Ω) ≤ Cm -s |f | H s (Ω) . (4.3)
The so-called Kohn-Nirenberg symbol N of H is defined for all (x, k) ∈ Ω × Z by Therefore:

N (x, k) = Ω T (x,
H m -H 2 HS = Ω Ω (T m (x, y) -T (x, y)) 2 dx dy (4.7) = Ω T m (x, •) -T (x, •) 2 L 2 (Ω) dx (4.8)
≤ κC 2 m -2s dx (4.9) (4.10)

As will be seen later, the convergence rate (4.6) is optimal in the sense that no productconvolution expansion of order m can achieve a better rate under the sole Assumptions 3.1 and 3.2. Corollary 4.2. Let > 0 and set m = C -1/s κ 1/2s . Under Assumptions 3.1 and 3.2, H m satisfies H -H m HS ≤ and products with H m and H * m can be evaluated with no more than O(κ 1/2s n log n -1/s ) operations.

Proof. Since Fourier atoms are not localized in the time domain, the modulation functions w k are supported on intervals of size p = n. The complexity of computing a matrix vector product is therefore O(mn log(n)) operations by Lemma 3.2.

Finally, let us mention that computing the discrete Kohn-Nirenberg N costs O(κn 2 log 2 (n)) operations (κn discrete Fourier transforms of size n). The storage cost of this Fourier representation is O(mκn) since one has to store κn coefficients for each of the m vectors h k .

In the next two sections, we show that replacing Fourier atoms by wavelet atoms or B-splines preserves the optimal rate of convergence in O( √ κm -s ), but has the additional advantage of being localized in space, thereby reducing complexity.

Spline decompositions

B-Splines form a Riesz basis with dual Riesz basis of form [START_REF] Christensen | An introduction to frames and Riesz bases[END_REF]:

( Bα,m (• -k/m)) 0≤k≤m-1 . (4.11)
The projection f m of any f ∈ L 2 (Ω) onto B α,m can be expressed as:

f m = arg min f ∈Bα,m f -f 2 2 (4.12) = m-1 k=0 f, Bα,m (• -k/m) B α,m (• -k/m). ( 4
f -f m 2 ≤ C √ κm -s f W s,2 . (4.14)
The following result directly follows.

Corollary 4.3. Set α ≥ s. For each x ∈ Ω, let (c k (x)) 0≤k≤m-1 be defined as

c k (x) = T (x, •), Bα,m (• -k/m) . ( 4

.15)

Define T m by: Proof. In this approximation, m B-splines are used to cover Ω. B-splines have a compact support of size (α + 1)/m. This property leads to windowing vector w k with support of size p = (α + 1) n m . Furthermore the vectors (h k ) have a support of size q = κn. Combining these two results with Lemma 3.2 and Corollary 4.3 yields the result for the choice α = s.

T m (x, y) = m-1 k=0 c k (x)B α,m (y -k/m). ( 4
The complexity of computing the vectors c k is O(κn 2 log(n)) (κn projections with complexity n log(n), see e.g. [START_REF] Unser | B-spline signal processing. Part I -Theory[END_REF]).

As can be seen in Corollary (4.4), B-splines approximations are preferable over Fourier decompositions whenever the support size κ is small.

Wavelet decompositions

f m = |µ|≤ log 2 (m) c µ ψ µ , (4.20)
where ψ is a Daubechies wavelet with α > s vanishing moments and c µ = ψ µ , f . Then

f m -f L 2 (Ω) ≤ Cm -s |f | H s (Ω) . (4.21)
A direct consequence is the following corollary.

T m (x, y) = |µ|≤ log 2 (m) c µ (x)ψ µ (y), (4.22)
where c µ (x) = ψ µ , T (x, •) . Then, under Assumptions 3.1 and 3.2

H m -H HS ≤ C √ κm -s . (4.23)
Proof. The proof is identical to that of Corollary (4.1).

Proposition 4.1. Let > 0 and set m = C -1/s κ Proof. In (4.22), the windowing vectors w k are wavelets ψ µ of support of size min((2s + 1)n2 -|µ| , n). Therefore each convolution has to be performed on intervals of size |ψ µ | + q + 1. Since there are 2 j wavelets at scale j, the total number of operations is:

µ | |µ|<log 2 (m) (|ψ µ | + q + 1) log 2 (min(|ψ µ |, q + 1)) (4.26) ≤ µ | |µ|<log 2 (m) ((2s + 1)n2 -|µ| + κn) log 2 (κn) (4.27) = log 2 (m)-1 j=0 2 j (2s + 1)n2 -j + κn log 2 (κn) (4.28) = log 2 (m)-1 j=0 (2s + 1)n + 2 j κn log 2 (κn) (4.29)
≤ ((2s + 1)n log 2 (m) + mκn) log 2 (κn) (4.30)

= (2s + 1)n log 2 ( -1/s κ 1/2s ) + -1/s κ 1+1/2s n log 2 (κn). (4.31)

Computing the vectors c µ costs O(κsn 2 ) operations (κn discrete wavelet transforms of size n). The storage cost of this wavelet representation is O(mκn) since one has to store κn coefficients for each of the m functions h k .

As can be seen from this analysis, wavelet and B-spline approximations roughly have the same complexity over the class T s . The first advantage of wavelets compared to B-splines is that the coefficients c µ (x) have a simple analytic expression, while B-splines coefficients c k are found by solving a linear system. This is slightly more complicated to implement.

The second significant advantage of wavelets compared to B-splines with fixed knots is that they are known to characterize much more general function spaces than H s (Ω). For instance, if all functions T (x, •) have a single discontinuity at a given y ∈ Ω, only a few coefficients c µ (x) will remain of large amplitude. Wavelets will be able to efficiently encode the discontinuity, while B-splines with fixed knots -which are not localized in nature -will fail to well approximate the TVIR. It is therefore possible to use wavelets in an adaptive way. This effect is visible on Figure 3c: despite discontinuities, only wavelets localized around the discontinuities yield large coefficients. In the next section, we propose two other adaptive methods, in the sense that they are able to automatically adapt to the TVIR regularity. Notice that the wavelet coefficients in kernel 1 rapidly decay as scales increase. The decay is slower for kernels 2 and 3 which are less regular. The adaptivity of wavelets can be visualized in Kernel 3: some wavelet coefficients are non zero at large scales, but they are all concentrated around discontinuities. Therefore only a few number of couples (c µ , ψ µ ) will be necessary to encode the discontinuities. This was not the case with Fourier or B-spline atoms.

Interpolation VS approximation

In all previous results, we constructed the functions w k and h k in 1.4 by projecting T (x, •) onto linear subspaces. This is only possible if the whole TVIR T is available. In very large scale applications, this assumption is unrealistic, since the TVIR contains n 2 coefficients, which cannot even be stored. Instead of assuming a full knowledge of T , some authors (e.g. [START_REF] Nagy | Restoring images degraded by spatially variant blur[END_REF]) assume that the impulse responses T (•, y) are available only at a discrete set of points y i = i/m for 1 ≤ i ≤ m.

In that case, it is possible to interpolate the impulse responses instead of approximating them. Given a linear subspace E m = span(e k , k ∈ {1, . . . , m}), where the atoms e k are assumed to be linearly independent, the functions c k (x) in (4.1) are chosen by solving the set of linear systems:

m k=1 c k (x)e k (y i ) = T m (x, y i ) for 1 ≤ i ≤ m. (4.32)
In the discrete setting, under Assumption 3.2, this amounts to solving κn linear systems of size m × m. We do not discuss the rates of approximation for this interpolation technique since they are usually expressed in the L ∞ -norm under more stringent smoothness assumptions than T ∈ T s . We refer the interested reader to [START_REF] Jacob | Cardinal spline interpolation[END_REF][START_REF] Devore | Constructive approximation[END_REF][START_REF] David | Interpolating wavelet transforms[END_REF] for results on spline and wavelet interpolants.

On Meyer's operator representation

Up to now, we only assumed a regularity of T in the y direction, meaning that the impulse responses vary smoothly in space. In many applications, the impulse responses themselves are smooth. In this section, we show that this additional regularity assumption can be used to further compress the operator. Finding a compact operator representation is a key to treat identification problems (e.g. blind deblurring in imaging). Since (ψ λ ) λ∈Λ is a Hilbert basis of L2 (Ω), the set of tensor product functions (ψ λ ⊗ψ µ ) λ∈Λ,µ∈Λ is a Hilbert basis of L 2 (Ω × Ω). Therefore, any T ∈ L 2 (Ω × Ω) can be expanded as:

T (x, y) = λ∈Λ µ∈Λ c λ,µ ψ λ (x)ψ µ (y). (4.33)
The main idea of the construction in this section consists of keeping only the coefficients c λ,µ of large amplitude. A similar idea was proposed in the BCR paper [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF] 2 , except that the kernel K was expanded instead of the TVIR T . Decomposing T was suggested by Beylkin at the end of [START_REF] Beylkin | On the representation of operators in bases of compactly supported wavelets[END_REF] without a precise analysis.

In this section, we assume that T ∈ H r,s (Ω × Ω), where

H r,s (Ω×Ω) = {T : Ω×Ω → R, ∂ α 1 x ∂ α 2 y T ∈ L 2 (Ω×Ω), ∀α 1 ∈ {0, . . . , r}, ∀α 2 ∈ {0, . . . , s}}. (4.34)
This space arises naturally in applications, where the impulse response regularity r might differ from the regularity s of their variations. Notice that

H 2s (Ω × Ω) ⊂ H s,s (Ω × Ω) ⊂ H s (Ω).
Theorem 4.2. Assume that T ∈ H r,s (Ω × Ω) and satisfies Assumption 3.2. Assume that ψ has max(r, s)

+ 1 vanishing moments. Let c λ,µ = T, ψ λ ⊗ ψ µ . Define H m 1 ,m 2 = |λ|≤log 2 (m 1 ) |µ|≤log 2 (m 2 ) c λ,µ ψ λ ⊗ ψ µ . (4.35)
Let m ∈ N, set m 1 = m s/(r+s) , m 2 = m r/(r+s) and

H m = H m 1 ,m 2 . Then H -H m HS ≤ C √ κm -rs r+s . (4.36)
Proof. First notice that

T ∞,m 2 = |µ|≤ log 2 (m 2 ) c µ ⊗ ψ µ , (4.37)
where c µ (x) = T (x, •), ψ µ . From Corollary 4.5, we get:

T ∞,m 2 -T L 2 (Ω×Ω) ≤ C √ κm -s 2 . (4.38)
Now, notice that c µ ∈ H r (Ω). Indeed, for all 0 ≤ k ≤ r, we get:

Ω (∂ k x c µ (x)) 2 dx (4.39) = Ω ∂ k x Ω T (x, y)ψ µ (y) dy 2 dx (4.40) = Ω Ω (∂ k x T )(x, y)ψ µ (y) dy 2 dx (4.41) ≤ Ω (∂ k x T )(x, •) 2 L 2 (Ω) ψ µ 2 L 2 (Ω) dx (4.42) = (∂ k x T ) L 2 (Ω×Ω) < +∞. (4.43)
Therefore, we can use Lemma 4.2 again to show:

T ∞,m 2 -T m 1 ,m 2 L 2 (Ω×Ω) ≤ C √ κm -r 1 . (4.44)
Finally, using the triangle inequality, we get:

T -T m 1 ,m 2 HS ≤ C √ κ(m -r 1 + m -s 2 ). (4.45)
By setting m 1 = m s/r 2 , the two approximation errors in the right-hand side of (4.45) are balanced. This motivates the choice of m 1 and m 2 indicated in the theorem.

The approximation result in inequality (4.36) is worst than the previous ones. For instance if r = s, then the bound becomes O( √ κm -s/2 ) instead of O( √ κm -s ) in all previous theorems. The great advantage of this representation is the operator storage: until now, the whole set of vectors (c µ ) had to be stored (O(κnm) values), while now, only m coefficients c λ,µ are required. For instance, in the case r = s, for an equivalent precision, the storage cost of the new representation is O(κm 2 ) instead of O(κnm).

In addition, evaluating matrix-vector products can be achieved rapidly by using the following trick:

H m u = |λ|≤log 2 (m 1 ) |µ|≤log 2 (m 2 ) c λ,µ ψ λ (ψ µ u) (4.46) = |µ|≤log 2 (m 2 )   |λ|≤log 2 (m 1 ) c λ,µ ψ λ   (ψ µ u). (4.47)
By letting cµ = |λ|≤log 2 (m 1 ) c λ,µ ψ λ , we get

H m u = |µ|≤log 2 (m 2 ) cµ (ψ µ u). (4.48)
which can be can be computed in O(m 2 κn log 2 (κn)) operations. This remark leads to the following proposition. Notice that the complexity of matrix-vector products is unchanged compared to the wavelet or spline approaches with a much better compression ability. However, this method requires a preprocessing to compute cµ with complexity -1/s κ 1/2s n.

Adaptive decompositions

In the last section, all methods shared the same principle: project T (x, •) on a fixed basis for each x ∈ Ω. Instead of fixing a basis, one can try to find a basis adapted to the operator at hand. This idea was proposed in [START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF] and [START_REF] Denis | Fast approximations of shift-variant blur[END_REF].

Singular value decompositions

The authors of [START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF] proposed to use a singular value decomposition (SVD) of the TVIR in order to construct the functions h k and w k . In this section we first detail this idea and then analyze it from an approximation theoretic point of view. Let J : L 2 (Ω) → L 2 (Ω) denote the linear integral operator with kernel T ∈ T s . First notice that J is a Hilbert-Schmidt operator since J HS = H HS . By Lemma 2.2 and since Hilbert-Schmidt operators are compact, there exists two Hilbert bases (e k ) and (f k ) of L 2 (Ω) such that J can be decomposed as (5.

J = k≥1 σ k • e k ⊗ f k , ( 5
2)

The following result is a standard.

Theorem 5.1. For a given m, a set of functions (h k ) 1≤k≤m and (w k ) 1≤k≤m that minimizes H m -H HS is given by:

h k = σ k f k and w k = e k . (5.3)
Moreover, if T (x, •) satisfies Assumptions 3.1 and 3.2, we get:

H m -H HS = O √ κm -s . (5.4)
Proof. The proof of optimality (5.4) is standard. Since T m is the best rank m approximation of T , it is necessarily better than bound (4.6), yielding (5.4).

Theorem 5.2. For all > 0 and m < n, there exists an operator H with TVIR satisfying 3.1 and 3.2 such that:

H m -H HS ≥ C √ κm -(s+ ) . (5.5)
Proof. In order to prove, (5.5), we construct a "worst case" TVIR T . We first begin by constructing a kernel T with κ = 1 to show a simple pathological TVIR. Define T by: is real for all (x, y). We now prove that T ∈ T s . The k-th Fourier coefficient of T (x, •) is given by σ k f k (x) which is bounded by σ k for all x. By Lemma 2.1, T (x, •) therefore belongs to H s (Ω) for all x ∈ Ω. By construction, the spectrum of T is (|σ k |) k∈N , therefore for any rank 2m + 1 approximation of T , we get:

T (x, y) = k∈Z σ k f k (x)f k (y), (5.6
T -T 2m+1 2 HS ≥ |k|≥m+1 1 |k| 2s+1+ (5.8) ≥ ∞ m+1 2 t 2s+1+ dt (5.9) = 1 2s + 2 (m + 1) 2s+
(5.10)

= O(m -2s-), (5.11) proving the result for κ = 1. Notice that the kernel K of the operator with TVIR T only depends on x:

K(x, y) = |k|≤N 2σ k cos(2πx).
(5.12)

Therefore the worst case TVIR exhibited here is that of a rank 1 operator H. Obviously, it cannot be well approximated by product-convolution expansions.

Let us now construct a TVIR satisfying Assumption 3.2. For this, we first construct an orthonormal basis ( fk ) k∈Z of L 2 ([-κ/2, κ/2]) defined by: fk

(x) = 1 √ κ f k x κ if |x| ≤ κ 2 , 0 otherwise. (5.13)
The worst case operator considered now is defined by:

T (x, y) = k∈Z σk fk (x)f k (y). (5.14)
Its spectrum is (|σ k |) k∈Z , and we get (5.16)

| T (x, •), f k | = |σ k fk (x)| = 1 κ |σ k |. ( 5
Even if the SVD provides an optimal decomposition, there is no guarantee that functions e k are supported on an interval of small size. As an example, it suffices to consider the "worst case" TVIR given in equation (5.6). Therefore, vectors w k are generically supported on intervals of size p = n. This yields the following proposition.

Corollary 5.1. Let > 0 and set m = C -1/s κ 1/2s . Then H m satisfies H -H m HS ≤ and a product with H m and H * m can be evaluated with no more than O(κ 1/2s n log n -1/s ) operations. Computing the first m singular vectors in (5.2) can be achieved in roughly O(κn 2 log(m)) operations thanks to recent advances in randomized algorithms [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF]. The storage cost for this approach is O(mn) since the vectors e k have no reason to be compactly supported. [START_REF] Denis | Fast approximations of shift-variant blur[END_REF] In [START_REF] Denis | Fast approximations of shift-variant blur[END_REF], the authors propose to construct the windowing functions w k and the filters h k using constrained optimization procedures. For a fixed m, they propose solving:

The optimization approach in

min (h k ,w k ) 1≤k≤m T - m k=1 h k ⊗ w k 2 HS
(5.17) under an additional constraint that supp(w k ) ⊂ ω k with ω k chosen so that ∪ m k=1 ω k = Ω. A decomposition of type 5.17 is known as structured low rank approximation [START_REF] Chu | Structured low rank approximation[END_REF]. This problem is non convex and to the best of our knowledge, there currently exists no algorithm running in a reasonable time to find its global minimizer. It can however be solved approximately using alternating minimization like algorithms.

Depending on the choice of the supports ω k , different convergence rates can be expected. However, by using the results for B-splines in section 4.2, we obtain the following proposition. Finally, let us mention that -owing to Corollary 4.5 -it might be interesting to use the optimization approach (5.17) with windows of varying sizes.

Summary and extensions 6.1 A summary of all results

Table 1 summarizes the results derived so far under Assumptions 3.1 and 3.2. In the particular case of Meyer's methods, we assume that T ∈ H r,s (Ω × Ω) instead of Assumption 3.1. As can be seen in this table, different methods should be used depending on the application. The best methods are:

• Wavelets: they are adaptive, have a relatively low construction complexity, and matrixvector products also have the best complexity.

• Meyer: this method has a big advantage in terms of storage. The operator can be represented very compactly with this approach. It has a good potential for problems where the operator should be inferred (e.g. blind deblurring). It however requires stronger regularity assumptions.

• The SVD and the method proposed in [START_REF] Denis | Fast approximations of shift-variant blur[END_REF] both share an optimal adaptivity. The representation however depends on the operator and it is more costly to evaluate it.

Method Approximation

4.1 O κ 1 2 m -s O κ 1 2s n log(n) -1 s O(κn 2 log(n)) O(mκn) B-Splines 4.2 O κ 1 2 m -s O κ 2s+1 2s n log(n) -1 s O(κn 2 log(n)) O(mκn) Wavelets 4.3 O κ 1 2 m -s O κ

Extensions to higher dimensions

Most of the results provided in this paper are based on standard approximation results in 1D, such as Lemmas 4.1, 4.2 and 4.1. All these lemmas can be extended to higher dimension and we refer the interested reader to [START_REF] Devore | Constructive approximation[END_REF][START_REF] Mallat | A wavelet tour of signal processing[END_REF][START_REF] Devore | Nonlinear approximation[END_REF][START_REF] Pinkus | N-widths in Approximation Theory[END_REF] for more details. We now assume that Ω = [0, 1] d and that the diameter of the impulse responses is bounded by κ ∈ [0, 1]. Using the mentioned results, it is straightforward to show that the approximation rate of all methods now becomes H -H m HS = O(κ d/2 m -s/d ).

(6.1)

The space Ω can be discretized on a finite dimensional space of size n d . Similarly, all complexity results given in Table 1 are still valid by replacing n by n d , -1/s by -d/s and κ by κ d .

Extensions to least regular spaces

Until now, we assumed that the TVIR T belongs to Hilbert spaces (see e.g. Assumption 3.1). This assumption was deliberately chosen easy to clarify the presentation. The results can most likely be extended to much more general spaces using nonlinear approximation theory results [START_REF] Devore | Nonlinear approximation[END_REF].

For instance, assume that T ∈ BV (Ω × Ω), the space of functions with bounded variations. Then, it is well known (see e.g. [START_REF] Cohen | Harmonic analysis of the space BV[END_REF]) that T can be expressed compactly on a Hilbert basis of tensor-product wavelets. Therefore, the product-convolution expansion 1.4 could be used by using the trick proposed in 4.48.

Similarly, most of the kernels found in partial differential equations (e.g. Calderòn-Zygmund operators) are singular at the origin. Once again, it is well known [START_REF] Meyer | Wavelets and operators[END_REF] that wavelets are able to capture the singularities and the proposed methods can most likely be applied to this setting too.

A precise setting useful for applications requires more work and we leave this issue open for future work.

Controls in other norms

In all the paper we only controlled the Hilbert-Schmidt norm • HS . This choice simplifies the analysis and also allows getting bounds for the spectral norm

H 2→2 = sup u L 2 (Ω) ≤1
Hu L 2 (Ω) , (

since H 2→2 ≤ H HS . In applications, it often makes sense to consider other operator norms defined by

H X→Y = sup u X ≤1
Hu Y , (

where • X and • Y are norms characterizing some function spaces. We showed in [START_REF] Escande | Sparse wavelet representations of spatially varying blurring operators[END_REF] that this idea could highly improve practical approximation results. Unfortunately, it is not clear yet how to extend the proposed results and algorithms to such a setting and we leave this question open for the future. Let us mention that our previous experience shows that this idea can highly change the method's efficiency.

Conclusion

In this paper, we analyzed the approximation rates and numerical complexity of productconvolution expansions. This approach was shown to be efficient whenever the time or space varying impulse response of the operator is well approximated by a low rank tensor. We showed that this situation occurs under mild regularity assumptions, making the approach relevant for a large class of applications. We also proposed a few original implementations of this methods based on orthogonal wavelet decompositions and analyzed their respective advantages precisely. Finally, we suggested a few ideas to further improve the practical efficiency of the method.

3 Figure 1 :

 31 Figure 1: Different kernels K, the associated TVIR T and the spectrum of the operator J. Left column corresponds to Example 3.1. Central column corresponds to Example 3.2. Right column corresponds to Example 3.3.

2 . 4 . 1 .

 241 y) exp(-2iπky) dy. (4.4) Illustrations of different Kohn-Nirenberg symbols are provided in Figure Corollary Set e k (y) = exp(-2iπky) and define T m by: T m (x, y) = |k|≤m N (x, k)e k (y). (4.5) Then, under Assumptions 3.1 and 3.2 H m -H HS ≤ C √ κm -s . (4.6) Proof. By Lemma 4.1 and Assumption 3.1, T m (x, •) -T (x, •) L 2 (Ω) ≤ Cm -s for some constant C and for all x ∈ Ω. In addition, by Assumption 3.2, T m (x, •) -T (x, •) L 2 (Ω) = 0 for |x| > κ/2.

Figure 2 :

 2 Figure 2: symbols of the kernels given in Examples 3.1, 3.2 and 3.3 in log 10 scale. Observe how the decay speed from the center (low frequencies) to the outer parts (high frequencies) changes depending on the TVIR smoothness. Note: the lowest values of the Kohn-Nirenberg symbol have been set to 10 -4 for visualization purposes.

Lemma 4 . 2 ([ 28 ,

 4228 Theorem 9.5]). Let f ∈ H s (Ω) and f m denote its partial wavelet series:

Corollary 4 . 5 .

 45 Let ψ be a Daubechies wavelet with α = s + 1 vanishing moments. Define T m by:

Figure 3 :

 3 Figure 3: "Wavelet symbols" of the operators given in Examples 3.1, 3.2 and 3.3 in log 10 scale. The red bars indicate separations between scales.Notice that the wavelet coefficients in kernel 1 rapidly decay as scales increase. The decay is slower for kernels 2 and 3 which are less regular. The adaptivity of wavelets can be visualized in Kernel 3: some wavelet coefficients are non zero at large scales, but they are all concentrated around discontinuities. Therefore only a few number of couples (c µ , ψ µ ) will be necessary to encode the discontinuities. This was not the case with Fourier or B-spline atoms.

3 Figure 4 :

 34 Figure 4: Meyer's representations of the operators in Examples 3.1, 3.2 and 3.3 in log 10 scale.

σ

  k f k (x)e k (y).

  ) where f k (x) = exp(2iπkx) is the k-th Fourier atom, σ 0 = 0 and σ k = σ -k = 1 |k| s+1/2+ /2 for |k| ≥ 1. With this choice, T (x, y) = |k|≤N 2σ k cos(2π(x + y))(5.7)

. 15 )T -T 2m+1 2 HS

 152 By Lemma 4.1, if σk = κ (1+|k| 2 ) s |k| 1+ , then T (x, •) H s (Ω)is uniformly bounded by a constant independent of κ. Moreover, by reproducing the reasoning in (5.8), we get: = O(κm -2s-).

Proposition 5 . 1 .

 51 Set ω k = [(k -1)/m, k/m + s/m] and let (h k , w k ) 1≤k≤m denote the global minimizer of (5.17). Define T m by T m (x, y) = m k=1 h k (x)w k (y). Then:T -T m 2 HS ≤ C √ κm -s . (5.18)Set m = κ 1/2s

- 1 s

 1 High (iterative) O(mκn)

  → L 2 (Ω) denote a compact operator. There exists two finite or countable orthonormal systems {e 1 , . . .}, {f 1 , . . .} of L 2 (Ω) and a finite or infinite sequence σ 1 ≥ σ 2 ≥ . . . of positive numbers (tending to zero if it is infinite), such that H can be decomposed as:

	Lemma 2.2 (Schmidt decomposition [33, Theorem 2.2] or [24, Theorem 1 p. 215]). Let H :
	L 2 (Ω)

  The proof is similar to that of Corollary (4.1). Let > 0 and set m = C -1/s κ 1/2s . Under Assumptions 3.1 and 3.2 H m satisfies H -H m HS ≤ and products with H m and H * m can be evaluated with no more than

	If α ≥ s, then, under Assumptions 3.1 and 3.2, operations. For small and large n, the complexity behaves like H Corollary 4.4. O s + κ 1+1/2s -1/s n log 2 (κn) O κ 1+1/2s n log 2 (κn) -1/s .	.16) (4.18) (4.19)

m -H HS ≤ C √ κm -s .

(4.17)

Proof.

  1/2s . Under Assumptions 3.1 and 3.2 H m satisfies H -H m HS ≤ and products with H m and H * m can be evaluated with no more than

	O sn log 2	-1/s κ 1/2s + κ 1+1/2s n -1/s log 2 (κn)	(4.24)
	operations. For small , the complexity behaves like
	O κ 1+1/2s n log 2 (κn) -1/s .	(4.25)

  C -1/s , then H m -H HS ≤ and the evaluation of a product with H m or H *

		m
	is of order
	O(κ 1+1/2s n log(n) -1/s ).	(5.19)

Proof. First notice that cardinal B-Splines are also supported on [(k -1)/m, k/m + s/m]. Since the method in

[START_REF] Denis | Fast approximations of shift-variant blur[END_REF]

provides the best choices for (h k , w k ), the distance H m -H HS is necessarily lower than that obtained using B-splines in Corollary 4.3.

Table 1 :

 1 Summary of the properties of different constructions. Approximation ≡ approximation rates in terms of m. Product ≡ matrix-vector product complexity to get an approximation. Construction ≡ complexity of the construction of order m representation. Storage ≡ cost of storage of a given representation. Adaptivity ≡ ability to automatically adapt to different input operators.

With the terminology of[START_REF] Busby | Product-convolution operators and mixed-norm spaces[END_REF], the name product-convolution would have been convolution-product and viceversa. The name product-convolution seems more appropriate to describe a product followed by a convolution.

This was also the basic idea in our recent paper[START_REF] Escande | Sparse wavelet representations of spatially varying blurring operators[END_REF].

Acknowledgments

We began investigating a much narrower version of the problem in a preliminary version of [18]. An anonymous reviewer however suggested that it would be more interesting to make a general analysis and we therefore discarded the aspects related to product-convolution expansions from [18]. We thank the reviewer for motivating us to initiate this research. We also thank Sandrine Anthoine, Jérémie Bigot, Caroline Chaux, Jérôme Fehrenbach, Hans Feichtinger, Clothilde Mélot, Anh-Tuan Nguyen and Bruno Torrésani for fruitful discussions on related matters. The

project funded by the PRES of Toulouse University and Midi-Pyrénèes region.