Mališa Vučinić
email: malisa.vucinic@imag.fr

Bernard Tourancheau
email: bernard.tourancheau@imag.fr

Andrzej Duda
email: andrzej.duda@imag.fr

Performance Comparison of the RPL and LOADng Routing Protocols in a Home Automation Scenario

RPL, the routing protocol proposed by IETF for IPv6/6LoWPAN Low Power and Lossy Networks has significant complexity. Another protocol called LOADng, a lightweight variant of AODV, emerges as an alternative solution. In this paper, we compare the performance of the two protocols in a Home Automation scenario with heterogenous traffic patterns including a mix of multipoint-to-point and point-to-multipoint routes in realistic dense non-uniform network topologies. We use Contiki OS and Cooja simulator to evaluate the behavior of the ContikiRPL implementation and a basic non-optimized implementation of LOADng. Unlike previous studies, our results show that RPL provides shorter delays, less control overhead, and requires less memory than LOADng. Nevertheless, enhancing LOADng with more efficient flooding and a better route storage algorithm may improve its performance.

I. INTRODUCTION

IETF standardization efforts for sensor networks have enabled their full interoperability and have made the envisioned "Internet of Things" (IoT) a reality. 6LoWPAN, the IPv6 over Low power Wireless Personal Area Networks protocol [START_REF] Montenegro | Transmission of IPv6 packets over IEEE 802.15.4 networks[END_REF] is the adaptation layer that allowed tiny devices to become reachable on the global IP network. In parallel, the IETF ROLL working group specified RPL, a Routing Protocol for Low Power and Lossy Networks (LLNs) [START_REF] Winter | RPL: IPv6 routing protocol for low power and lossy networks[END_REF]. Many argued that even if RPL met the initial goal of providing the industry actors with a fully-fledged routing protocol, its significant complexity poses a threat for implementation interoperability of constrained devices. Another routing proposal is LOADng (The Lightweight On-demand Ad hoc Distance-vector routing protocol -Next Generation) [START_REF] Clausen | The Lightweight On-demand Ad hoc Distance-vector Routing Protocol -Next Generation (LOADng)[END_REF], a lightweight variant of AODV [START_REF] Perkins | Ad hoc On-Demand Distance Vector (AODV) Routing[END_REF]. The aim of this paper is to analyze the performance characteristics and issues of the two routing protocols for a specific scenario of Home Automation networks.

Routing in LLNs is one of the key challenges for the IoT emergence. The constraints of LLNs have a significant impact on the protocol design. Small memory limits the number of stored route entries. Limited energy supply dictates minimal radio usage and optimized control overhead. The increasing scale of IoT networks calls for scalable solutions. Finally, Home Automation interactive applications may require low latency communications.

Much research work focused on the performance of RPL [START_REF] Tripathi | Performance evaluation of routing protocol for low power and lossy networks[END_REF], [START_REF] Ko | Evaluating the performance of RPL and 6LoWPAN in TinyOS[END_REF], [START_REF] Saad | Simulation of the RPL Routing Protocol for IPv6 Sensor Networks: two cases studies[END_REF], [START_REF] Accettura | Performance analysis of the RPL routing protocol[END_REF] leading to the common observation that RPL performs well in case of multipoint-to-point traffic, but induces a large overhead in scenarios where point-tomultipoint traffic is non-negligible. We show how the Contiki RPL implementation behaves in a realistic Home Automation scenario. LOADng is a more recent protocol compared to RPL and less studied. Our paper compares both protocols and gives new insights into their performance, existing issues, and draws possible research directions.

We have used the framework for low power IPv6 routing simulation, experimentation, and evaluation based on Contiki OS and its Cooja simulator/emulator [START_REF] Tsiftes | A framework for low-power IPv6 routing simulation, experimentation, and evaluation[END_REF]. Cooja uses a hardware emulator to run the code, so the simulation results perfectly reflect the behavior of the underlying protocols. The only part of the simulation that approximates real-world conditions is the radio propagation model with the assumption of the Unit Disk Graph. We strictly follow the Home Automation scenario [START_REF] Brandt | Home Automation Routing Requirements in Low-Power and Lossy Networks[END_REF] by mimicking the logical roles of nodes and appropriate traffic patterns, as well as by generating realistic topologies on a virtual home plan.

The paper is organized as follows. Sections II and III give an overview of RPL and LOADng, respectively. Sections IV and V discuss the simulation scenario and the performance results. Section VI presents the related work. Finally, we provide some concluding remarks, observations, and perspectives in Section VII.

II. RPL-ROUTING PROTOCOL FOR LOW POWER AND LOSSY NETWORKS

RPL is a Distance Vector protocol that specifies how to construct a Destination Oriented Directed Acyclic Graph (DODAG) with a defined objective function and a set of metrics and constraints. RPL uses a proactive approach: it finds and maintains routes without any traffic considerationsroutes are created even if not used.

RPL specifies a set of new ICMPv6 control messages to exchange information related to a DODAG:

• DODAG Information Solicitation (DIS) messages proactively solicit the DODAG related information from neighboring nodes.

• DODAG Information Object (DIO) defines and maintains upward routes.

• DODAG Destination Advertisement Object (DAO) advertizes prefix reachability towards the leaf nodes of a DODAG enabling downward traffic. A root starts the DODAG building process by transmitting a DIO. Neighboring nodes process DIOs and make a decision on joining the DODAG based on the objective function and/or local policy. A node computes its Rank with respect to the root and starts advertising DIO messages to its neighbors with the updated information. As the process converges, each node in the network receives one or more DIO messages and has a preferred parent towards the sink. Hence, RPL optimizes the upward routes for multipoint-to-point traffic that accounts for most of the traffic in LLNs.

To support downward routes, RPL uses DAO control messages that give the prefix information, the route lifetime, and other information about the distance of the prefix. RPL RFC [START_REF] Winter | RPL: IPv6 routing protocol for low power and lossy networks[END_REF] defines the storing and non-storing modes. In the nonstoring mode, packets use source-routing for downward traffic. In our study, we focus on the storing mode in which each node keeps track of all accessible downlink prefixes.

The Trickle algorithm [START_REF] Levis | The Trickle Algorithm[END_REF] governs the emission interval of DIOs. The idea is to reduce the control overhead of the protocol by sending DIOs less frequently when there is no change in the topology. In case of a change in the network, trickle forces more frequent emissions of DIOs. The RPL RFC [START_REF] Winter | RPL: IPv6 routing protocol for low power and lossy networks[END_REF] does not specify the mechanisms for the DAO emission (it is left to implementation). ContikiRPL emits DAOs with a similar approach to the trickle algorithm based on the DIO transmission timers.

III. LOADNG

The LOADng protocol [START_REF] Clausen | The Lightweight On-demand Ad hoc Distance-vector Routing Protocol -Next Generation (LOADng)[END_REF] uses a reactive approach based on the idea that LLNs are idle most of the time so a proactive approach would generate unnecessary overhead. Thus, LOADng establishes a route towards a given destination only on demand when there is some data to send. As IETF did a lot of work on designing AODV, a reactive protocol for MANETs [START_REF] Perkins | Ad hoc On-Demand Distance Vector (AODV) Routing[END_REF], a logical consequence was to adapt it for LLNs to make it implementable on memory constrained devices.

When a device has a packet to send towards a given destination, it consults a routing table and invokes LOADng in case of an invalid entry. The protocol floods a Route-Request (RREQ) message through the network to reach all nodes. A node receiving a RREQ checks if it is the message destination. If not, it forwards RREQ to its neighbors. The node also learns the reverse path towards the originator of the RREQ message and adds it to the routing table. Eventually, the destination node receives RREQ and responds by unicasting a Route-Reply (RREP) message towards the request originator. RREP follows the stored reverse route. At the same time, intermediate nodes learn the forward route towards the destination. When RREP reaches the request originator, the bidirectional route is installed at intermediate nodes.

One of the main drawbacks of LOADng is the route discovery delay. During the discovery process, outgoing packets are buffered, which may cause losses in memory constrained devices. Moreover, flooding is highly energy inefficient so nodes may suffer from energy depletion. Another issue is related to the collisions of control messages due to flooding, which may lead to unnecessary retransmissions.

We have used a LOADng implementation developed in our lab based on the AODV implementation in Contiki.

IV. SIMULATION SCENARIO

Home automation is one of the key applications of the envisioned "Internet of Things". Nodes in such networks are typical highly constrained devices. Good examples are light dimmers, window shades, motion sensors, typical monitoring sensors, remote control units. When it comes to the network architecture, two approaches are generally used [START_REF] Brandt | Home Automation Routing Requirements in Low-Power and Lossy Networks[END_REF]:

• Centralized architecture: each networked device communicates with a central node that controls the network.

• Distributed architecture: different devices within the network may cooperate to locally control the network. The majority of existing deployments uses the centralized approach, but the distributed approach is also gaining popularity, mainly because of better latency. However, in a futuristic scenario, "smart homes" are expected to be a part of a "smart grid" and a "smart city" with dynamical adjustments of their energy consumption and production. Thus, the central unit playing the role of a coordinator among home devices and the external ones, has to exist. For this reason, we have decided to use the centralized architecture for our study.

A. Topology

A very common simplification in the literature is the assumption of a grid-like network topology. Its main consequence is a near-uniform distribution of the number of node neighbors (node degree) that does not reflect common "smart home" deployment scenarios. Furthermore, the simplification strongly impacts the performance of routing protocols. To obtain realistic results, we have designed a virtual 130 m 2 home map and developed a realistic topology generator. The home map is the background for generating topologies of different sizes. Nowadays, networks of around 20 nodes are quite common. However, this number is expected to grow and may exceed 250 in the near future [START_REF] Brandt | Home Automation Routing Requirements in Low-Power and Lossy Networks[END_REF]. Due to the fixedsize area, topologies with a larger number of devices in a network become more dense and consequently, the number of neighbors increases.

We have generated topologies under the assumption of having 80 % of devices uniformly placed on the walls of a room, while the remaining 20 % are uniformly distributed around. The number of devices in each room is proportional to its area. The controlling unit (sink) is placed in the living room. 70 % of the devices have the logical role of sensors, while the remaining 30 % are assigned the role of actuators. 70 % of all sensors perform monitoring while the other 30 % generate events (switches, motion detectors).

Fig. 1 shows an example network topology with 25 nodes in the virtual house. The average number of neighbors for the generated topologies ranged from 4.13 for a 15 node topology, 6.24 for a 25 node network, up to 12.2 for the largest simulated network composed of 40 nodes.

B. Traffic Pattern

The traffic pattern of each node in the network depends on its assigned logical role. We have defined four logical types summarized in Table I. It is important to note that we have chosen the reporting periods and the general traffic patterns as suggested by IETF [START_REF] Brandt | Home Automation Routing Requirements in Low-Power and Lossy Networks[END_REF]. Nodes use UDP as the transport layer protocol and the application process acknowledges each received message. Therefore, traffic in the network is evenly distributed between point-to-multipoint and multipointto-point.

C. Simulation Parameters

To closely reflect real deployments, nodes run Contiki OS with the whole protocol stack. We have compiled Contiki for the Tmote Sky platform based on MSP430 microcontroller with 10 KB of RAM and an IEEE 802.15.4 compliant radio [START_REF] Moteiv | Tmote Sky -ultra low power IEEE[END_REF]. By using the MSP430 hardware emulator, Cooja thus takes into account all the hardware constraints of the devices under study. Table II Notice a fairly low communication range at the PHY layer. Indeed, we have chosen this setup to reflect the assumption that nodes will operate with a very low power, i.e. with a low power amplifier gain mainly to reduce interference. Furthermore, radio-propagation obstacles present in a typical home environment, coupled with 2.4 GHz frequency, limit the communication range. With our parameters, the whole house is covered within four hops as suggested before [START_REF] Brandt | Home Automation Routing Requirements in Low-Power and Lossy Networks[END_REF]. We have modified Contiki CSMA to buffer multicast packets and retransmit them when the channel was found busy at the initial attempt. Without this modification, the performance of LOADng was highly degraded due to a large number of dropped RREQ messages.

Another realistic assumption taken into account in our simulation was the use of ContikiMAC, a Radio Duty Cycling (RDC) protocol [START_REF] Dunkels | The ContikiMAC radio duty cycling protocol[END_REF] based on preamble sampling and possibly sleeping receivers. Note that ContikiMAC operates on top of the IEEE 802.15.4 non-slotted CSMA and it may lead to some lost frames, which alleviates the assumption of the Unit Disk Graph that does not introduce any losses for nodes in the radio range. Table III summarizes RPL and LOADng protocol parameters. Notice that we vary the Route Hold Time (RHT) of LOADng to study the protocol performance as a function of the route lifetime. RPL uses the Trickle algorithm [START_REF] Levis | The Trickle Algorithm[END_REF] to emit DIO messages with default parameters and we have set the DIO Max Interval used in the steady state to approximately 17 minutes (1048 s). Tuning of RHT and the DIO Max Interval is an engineering challenge that should take into account the dynamic behavior of a given network, loss rate, and the probability of a node failure. Our choices for the LOADng parameters cover a broad range of cases. We can still decrease the overhead of RPL by increasing the DIO Max Interval.

We average the simulation results over 5 simulation runs and show 95 % confidence intervals. We plot CDF graphs using the cumulative results of 5 simulation runs.

V. PERFORMANCE EVALUATION

We first focus on a 25 node network (cf. Fig. 1) and evaluate LOADng and RPL in terms of packet delays, hop counts, and the routing table size. Furthermore, we study the control plane overhead and the average routing table size as a function of the number of nodes to show how the protocols scale for larger networks.

A. Packet Delay

The latency studied in this section represents the total delay that a packet experiences from the instant it is passed to the UDP layer until it is successfully received at the destination. Fig. 2 presents the Cumulative Distribution Function (CDF) of packet delays for different hop counts. LOADng performs worse and packets experience significantly higher delays. Even though some LLNs can be delay insensitive, the interactive nature of home automation applications requires a routing protocol with delays of less than 0.5 seconds [START_REF] Brandt | Home Automation Routing Requirements in Low-Power and Lossy Networks[END_REF]. Due to the route discovery procedure of LOADng, the delay CDF converges slowly and 56.8 % of packets are within 0.5 seconds in case of a 10 minute RHT. Flooding is less frequent when RHT increases to one hour and the delay becomes smaller. However, the delay CDF still converges slowly and approximately 26 % of packets experience a delay greater than 0.5 seconds (cf. Fig. 2a). The pro-active route discovery of RPL results in shorter delays and 90.9 % of packets experience a delay less than 0.5 seconds in case of the one hop distance. Moving farther away from the sink does not significantly increase packet delays of RPL. However, we have observed a significant delay increase for LOADng. More precisely, 48.2 % of packets experience a delay greater than 0.5 seconds (cf. Fig. 2c). In case of a lower, 10 minute RHT, approximately 62 % of packets undergo a significant delay. The reason for the delay increasing with the hop distance is the flooding operation of LOADng as well as the instillation of forward routes. Basically, a node in the network sends a RREQ with the target address of the sink. Due to flooding, the sink receives more than one RREQ, each containing different path metrics. It responds with a RREP to the first RREQ and to all other RREQs with lower metrics. Nodes that form the reverse route and forward the RREPs towards its final destination instill the route towards the sink. The sink, however, does not have a route towards the intermediate nodes and upon the reception of a packet from them, has to flood the network with another RREQ before sending the acknowledgment. This results in an increased delay with respect to packets destined to nodes physically closer to the sink.

B. Hop Distance

We have studied the optimality of the constructed routes through path hop and packet hop distances. The path hop distance represents the average number of hops between a source and a destination. As the traffic in the network is either multipoint-to-point or point-to-multipoint, the sink node is either a source or a destination of each packet. It is interesting to note that while RPL has built only bidirectional paths, some paths in LOADng networks are unidirectional, because of the stochastic nature of the flooding process. More precisely, most nodes add a route to the sink during the first sink-originated RREQ flooding. However, the reverse path is found at a later time, when the sink needs to address the node. Fig. 3a presents the CDF of path hop distances in the network. Notice that the RPL routes are shorter than those of LOADng. As mentioned in Section IV-A, the whole house is covered within four hops. The figure shows that the average path hop distances of RPL are bounded within four hops, which is not the case for LOADng. Moreover, 29 % of RPL routes are within one hop, while it is around 21 % for the LOADng, in the case of 10 minute RHT. This gap is even larger for two hop paths: 67 % of paths in case of RPL, while only 43 % in case of LOADng. Note that LOADng provides its most optimized routes in the case of a shorter RHT. There are two main reasons for this: i) due to the flooding operation the channel is saturated and some RREQ messages are dropped, because they reach the maximal number of attempts of the CSMA protocol, ii) as the nodes are memory constrained, they are only able to store two packets at a time, so RREQ messages that should be forwarded are often dropped. As a consequence, more frequent flooding (shorter RHT) supersedes the lost RREQs and gives more optimized routes.

Another issue is the following. If a node has a packet to send and the route is not ready, the most desired behavior, the one implemented in our case, is to buffer the packet. Once the RREP arrives, the node transmits the packet. However, the first arriving RREP is not necessarily the optimal one in terms of the route length. Therefore, the packet will follow a non-optimized route. We can notice in Fig. 3b that a nonnegligible number of packets go over more than four hops. Upon reception of a subsequent RREP, the node will update its routing table and the following packets will follow the updated route. On the other hand, routes constructed with RPL are optimal: not a single packet goes over more than 4 hops.

C. Routing Table Size

Due to the memory constraints of LLN devices, the routing table size is an important aspect. Fig. 4 shows the results for a 25 node network. LOADng requires much higher number of entries with respect to RPL. Again, as a consequence of flooding, each node receiving a RREQ or on a route of a RREP, instills a route towards the sender, which results in a large number of unnecessary routes. Protocol implementors should thus care about the priority of RREP routes. This fact could endanger the operation of the protocol if a node runs out of the available memory. On the other hand, most nodes in the RPL network just have a default entry towards the preferred parent. However, depending on their position in the network, some nodes also have a significant number of entries. Namely, as the RPL operates in the storing mode, intermediate nodes selected as preferred parents by others, have to store downward routes. This is a critical issue, because if such a node runs out of memory, a loop may be formed. To see how the protocols scale for a larger network, we have studied the average number of route entries as a function of an increasing number of nodes (cf. Fig. 5a). We can notice that the main consequence of having a long RHT is a significantly increasing number of route entries. With a shorter lifetime, the protocol scales better as the routes expire faster and the average number of required entries slowly increases. For RPL, the average number of routes slowly increases. However, it has a significantly higher number of entries for 40 nodes, a consequence of the limited number of neighbor entries that nodes can store. During the simulation time, we have observed oscillations in the RPL graph as parent nodes were constantly overwritten by other neighbors so that a number of nodes played the role of the preferred parent and stored downward routes. This behavior may impact the RPL operation for more dense networks.

D. Overhead

We have evaluated the routing overhead of the protocols as a function of an increasing number of nodes (cf. Fig. 5b presenting the overhead in bytes during the simulation time). Although the previous work [START_REF] Herberg | A comparative performance study of the routing protocols LOAD and RPL with bi-directional traffic in lowpower and lossy networks[END_REF] revealed a significantly higher overhead with RPL, our results show that RPL benefits from fairly low overhead compared to LOADng. We evaluate the total overhead as a sum of the lengths of all control messages passed from the routing protocols to the layers bellow. Note that by doing so, we avoid taking into account the overhead introduced by ContikiMAC.

We can notice that the implementation of ContikiRPL results in a lower overhead than observed previously [START_REF] Herberg | A comparative performance study of the routing protocols LOAD and RPL with bi-directional traffic in lowpower and lossy networks[END_REF]. As the RPL RFC under-specifies the DAO control message emission, the overhead strongly depends on a given implementation. Furthermore, it is important to note that our simulation scenario does not trigger any global repairs of the RPL DODAG during the eight hours of simulation. In this way, our results give an idea about the performance of the two protocols in the steady state. Thus, we observe that most of the RPL overhead appears at the beginning of the simulation, reduced later by Trickle. The overhead of LOADng strongly depends on RHT: for a shorter RHT, nodes flood more frequently. We can notice from the figure that the total amount of the LOADng overhead sharply increases with the number of nodes due to the high density of nodes in the network, which may however, be reduced with an optimized flooding algorithm at the cost of a higher complexity.

The lifetime of battery-powered nodes directly relates to the use of the radio transceiver. More precisely, in case of the Tmote Sky platform, power consumption due to the radio use is three orders of magnitude larger than that due to CPU processing [START_REF] Moteiv | Tmote Sky -ultra low power IEEE[END_REF]. Thus, the control plane overhead and average number of hops of a protocol are the determining factors of the expected node lifetime. Our study shows that the control overhead of LOADng for short hold times does not scale well with dense deployments of smart home applications. As we expect an increasing number of devices in home networks, this effect may significantly impact the overall network lifetime.

VI. RELATED WORK

Several authors analyzed the performance of RPL. An IETF draft [START_REF] Tripathi | Performance evaluation of routing protocol for low power and lossy networks[END_REF] evaluates the protocol by considering several routing metrics in real-life deployment scenarios. Nuvolone studied stability delays of RPL using OmNet++ [START_REF] Nuvolone | Stability analysis of the delays of the routing protocol over low power and lossy networks[END_REF]. Clausen et al. studied the multipoint-to-point performance of RPL using NS2 [START_REF] Clausen | Multipoint-to-point and broadcast in RPL[END_REF]. The authors [START_REF]Comparative study of RPL-enabled optimized broadcast in wireless sensor networks[END_REF] evaluated different optimized broadcast techniques for the use in RPL [START_REF]Comparative study of RPL-enabled optimized broadcast in wireless sensor networks[END_REF]. Ko et al. presented the implementation of RPL for TinyOS and discussed its performance [START_REF] Ko | Evaluating the performance of RPL and 6LoWPAN in TinyOS[END_REF]. Ben Saad et al. used the Contiki framework to evaluate the performance of a PLC network serving as a LLN backbone [START_REF] Saad | Simulation of the RPL Routing Protocol for IPv6 Sensor Networks: two cases studies[END_REF]. Other authors studied the routing overhead and delay induced by ContikiRPL implementation [START_REF] Accettura | Performance analysis of the RPL routing protocol[END_REF]. Finally, an IETF draft [START_REF] Clausen | Experiences with RPL: IPv6 routing protocol for low power and lossy networks[END_REF] presented a thorough study of several critical issues in the RPL protocol.

Gomez et al. [START_REF] Gomez | Adapting AODV for IEEE 802.15.4 mesh sensor networks: theoretical discussion and performance evaluation in a real environment[END_REF] compared various AODV-based routing protocols, including LOAD and discussed design tradeoffs. Iliev et al. [START_REF] Iliev | Mesh routing for low-power mobile adhoc wireless sensor networks using LOAD[END_REF] evaluated the route discovery delay as well as the round trip time induced by the TinyOS implementation of LOAD in real-life deployment scenarios.

Herberg et al. presented a comparative performance study of RPL and LOADng in case of bidirectional traffic using simulations in NS2 [START_REF] Herberg | A comparative performance study of the routing protocols LOAD and RPL with bi-directional traffic in lowpower and lossy networks[END_REF]. The paper shows a significantly larger control overhead of RPL caused by the maintenance of downward routes. It also compares the two protocols with an ideal routing protocol to show that RPL provides near-optimal routes while LOADng results in a certain gap. As the RPL RFC does not specify the period or the mechanism to use for maintaining downward routes, the study assumed an interval of 15 seconds. This choice is questionable, because it is the main cause of the high control overhead of RPL. Furthermore, the application-layer scenario used for the comparison is not the same for the two protocols. Thus, the question remains how LOADng and RPL perform under the same scenario.

VII. CONCLUSION

In the light of our study, RPL results in good overall performance for our Home Automation scenario while LOADng could serve better sparse LLN deployments with low-priority traffic in which the Route Hold Time can attain a large value. For Home Automation applications in which the response time is important, our results suggest that LOADng is not the best candidate. However, it is important to recall that we have used a simple flooding scheme and a better mechanism may reduce the control overhead. Furthermore, an intelligent route storing algorithm may reduce memory requirements of LOADng. These aspects are open research issues and their solutions will improve the efficiency and scalability, but at the cost of a higher implementation complexity.

We have confirmed the previous findings that RPL provides shorter routes and a smaller spanning tree depth compared to LOADng in our rather dense Home Automation topologies. We have also showed that RPL delays stay small and its overhead strongly depends on the implementation and a careful choice of parameters. Moreover, RPL has shown better results and less memory requirements than LOADng, but at the price of a higher implementation complexity.

To sum up, several aspects remain interesting research challenges: reducing the specification complexity of RPL, optimizing flooding and route storage of LOADng as well as mixing the proactive and reactive approaches to maximize some performance criteria for a given application domain.

Fig. 1 .

 1 Fig. 1. An example 25 node topology generated in a virtual house.

Fig. 2 .

 2 Fig. 2. Packet delay CDFs in a 25 node network. Packets going through one intermediate node, i.e. two hops away from the sink, (cf. Fig 2b) experience higher delay, as expected. In case of RPL, 15.4 % of packets experience a delay greater than 0.5 seconds. The maximal delay in this case is 2.4 seconds. With respect to the 1 hop distance delays, the delay CDF of RPL converges slightly slower due to the processing time and the channel contention at the intermediate node, because nodes have already a route. In case of LOADng and 1 hour RHT, 31 % of packets experience a delay greater than 0.5 sec.

  LOADng 30 min hold LOADng 10 min hold (a) CDF of path hop distances.

  LOADng 30 min hold LOADng 10 min hold (b) CDF of packet hop distances.

Fig. 3 .

 3 Fig. 3. Distance CDFs in the 25 node network.

Fig. 4 .

 4 Fig. 4. Average number of routing table entries during the simulation.

  Average number of route entries.

  Control plane overhead (bytes).

Fig. 5 .

 5 Fig. 5. Memory and overhead for an increased number of nodes.

TABLE I TRAFFIC

 I PATTERN OF DIFFERENT NODE TYPES IN THE NETWORK.

	Node Type	Traffic Pattern
	Monitoring Sensor	Periodic reporting in [8, 12] minute interval
	Event Sensor	Poisson process with a mean of 10 packets per
		hour for the whole house
	Actuator	Periodic reporting in [8, 12] minute interval and
		sending acknowledgment frames
	Main Controller Unit	Acknowledgment of received frames; sending a
		5 packet burst to different actuators as a Poisson
		process with a mean of 10 bursts per hour and
		upon reception of a frame from an Event Sensor

  summarizes the Contiki OS and Cooja setups.

TABLE II CONTIKI

 II OS AND COOJA PARAMETER SETUP.

	Settings	Value
	Wireless channel model	UDG Model with Distance Loss
	Communication range	5 m
	Mote type	Tmote Sky
	Transport and network layers	UDP + µIPv6 + 6LoWPAN
	Max number of queued packets	2
	MAC layer	non-slotted CSMA + ContikiMAC
	Radio interface	CC2420 2.4 GHz (IEEE 802.15.4)
	Simulation time	8h

ACKNOWLEDGEMENTS

This work was partially supported by the European Commission project CALIPSO under contract 288879, the French National Research Agency (ANR) projects ARESA2 under contract ANR-09-VERS-017, and IRIS ANR-11-INFR-016. Many thanks to Chi-Anh La for making available the LOADng implementation.


Steve Y Oudot
email: steve.oudot@inria.fr

Donald R Sheehy
email: don.r.sheehy@gmail.com

Zigzag Zoology: Rips Zigzags for Homology Inference

Keywords: F.2.2 [Analysis of Algorithms and Problem Complexity]: Nonnumerical Algorithms and Problems-Computations on discrete structures, Geometrical problems and computations Algorithms, Theory Zigzag Persistence, Persistent Homology, Persistence Barcodes, Rips Filtration, Weak Feature Size

For points sampled near a compact set X, the persistence barcode of the Rips filtration built from the sample contains information about the homology of X as long as X satisfies some geometric assumptions. The Rips filtration is prohibitively large, however zigzag persistence can be used to keep the size linear. We present several species of Ripslike zigzags and compare them with respect to the signalto-noise ratio, a measure of how well the underlying homology is represented in the persistence barcode relative to the noise in the barcode at the relevant scales. Some of these Rips-like zigzags have been available as part of the Dionysus library for several years while others are new. Interestingly, we show that some species of Rips zigzags will exhibit less noise than the (non-zigzag) Rips filtration itself. Thus, Rips zigzags can offer improvements in both size complexity and signal-to-noise ratio.

Along the way, we develop new techniques for manipulating and comparing persistence barcodes from zigzag modules. We give methods for reversing arrows and removing spaces from a zigzag while controlling the changes occurring in its barcode. We also discuss factoring zigzags and a kind of interleaving of two zigzags that allows their barcodes to be compared. These techniques were developed to provide our theoretical analysis of the signal-to-noise ratio of Rips-like zigzags, but they are of independent interest as they apply to zigzag modules generally.

INTRODUCTION

The goal of homology inference is to extract the homology of a space from a finite sample. The problem is ill-posed in general, but under the right geometric assumptions about the input and the underlying space, one can compute an object called a persistence barcode which provably contains information about the underlying homology. Indeed, homology inference was and continues to be one of the main motivations for topological persistence theory.

The barcode is computed from a sequence of simplicial complexes, for which two main challenges arise. The first challenge is to guarantee that the simplicial complexes remain small. Commonly used methods produce complexes that quickly become too large to fit in memory. The second challenge is to decrease noise in the barcode while preserving the signal, i.e. the information about the underlying space. We confront both challenges, analyze several approaches that give linear size data structures, and provide guarantees on the signal-to-noise ratio in the barcodes.

Context. Persistent homology applies to nested, parameterized families of simplicial complexes called filtrations. The persistence algorithm takes a filtration and produces a barcode describing all the changes in homology as one goes from one complex to the next in the filtration [START_REF] Edelsbrunner | Topological persistence and simplification[END_REF][START_REF] Zomorodian | Computing persistent homology[END_REF].

Persistent homology has an important connection with geometric inference results that describe conditions when homology inference is possible using a union of balls centered at the sample points-see the survey by Chazal and Cohen-Steiner [START_REF] Chazal | Geometric inference[END_REF]. The (Vietoris-)Rips filtration {Rα} α≥0 is useful when these conditions are met. It is defined to have a simplex in Rα for every subset of points with diameter at most α. So, the filtration parameter is the geometric scale and the theory guarantees the existence of some range of scales for which the barcode encodes the homology of the underlying space. The barcode of this filtration thus has an elegant multi-scale interpretation as "the homology of the input point cloud across scales."

The immediate drawback to using the Rips filtration is its size. The scale at which it exceeds the available memory varies with the input data, the filtration, and the computer used. However, it is observed to happen early enough so that not all the interesting homological information hidden in the data can be discovered-see Section 6 for a compelling ex-ample. Recent research looks at how to reduce the size of the complexes in the filtration to postpone the breaking point. The most notable example is the witness complex [START_REF] De | A weak characterisation of the Delaunay triangulation[END_REF], which was introduced for this purpose. However, such attempts are limited as the complexes have controlled size only at small scales and still incur a blowup at large scales.

A different approach was proposed by Chazal and Oudot [START_REF] Chazal | Towards persistence-based reconstruction in Euclidean spaces[END_REF], who used truncated Rips filtrations on a nested sequence of subsets of the input points corresponding to samplings at different scales. Their method computes the barcodes of the Rips filtration of each subset restricted to a range of scales near the sampling scale of the subset. This can prevent the size blowup in the Rips filtrations because every subset looks like a uniform sample at the relevant scale. The lingering challenge from this work is to relate the bars in the resulting barcodes for different scales.

Taking advantage of the recent introduction of zigzag persistence by Carlsson and de Silva [START_REF] Carlsson | Zigzag persistence[END_REF], Morozov suggested a simple way to connect the truncated Rips filtrations of consecutive subsamples together, to obtain a single long sequence of simplicial complexes connected by inclusionscalled the Morozov zigzag (M-ZZ) hereafter. Zigzag persistence relaxes the condition that the family of complexes be a filtration and instead allows consecutive spaces to be included in either direction, forwards or backwards, so the sequence is a zigzag diagram rather than a filtration. The M-ZZ has been integrated into the Dionysus library [START_REF]By Dmitriy Morozov[END_REF] since early 2009, and as reported by its author from preliminary experiments [START_REF] Morozov | [END_REF], it has given surprisingly good results in practice. However, to date it comes with no theoretical guarantees, so a primary motivation of our paper is to assess the theoretical quality of the results provided by this zigzag.

Existing methods for building sparse approximations to the Rips filtration have all focused on the size question, but have ignored the question of noise. For example, even the Rips filtration can have noise in the barcode at the scales where it represents the underlying topology. We show that some variants of the M-ZZ not only recover the topological signal but also provably eliminate noise in the relevant range.

Contributions. We provide the following theoretical guarantees for the Morozov zigzag:

• When the input point cloud P is sufficiently close (in the Hausdorff distance) to a compact set X with positive weak feature size in R d , there is a sweet range of geometric scales over which the persistence barcode of the Morozov zigzag exhibits the homology of X (technically, the offsets X λ for an arbitrarily small λ > 0). That is, the barcode has long intervals spanning the entire sweet range, and their number is at least the dimension of the homology group HX λ (Theorem 4.2). • If X has positive µ-reach, then there is a smaller (sweeter) range over which the number of spanning intervals is exactly dim HX λ and no other intervals are present. This motivates the study of more elaborate variants of the Morozov zigzag that are less likely to carry topological noise in the sweet range, even when the underlying space X has zero µ-reach and positive weak feature size. We analyze three variants in the paper:

• The first one, called the discretized Mozorov zigzag, consists in considering only subsamples whose corresponding geometric scales are of the form ζ i for a fixed constant ζ and an integer i. This discretization makes sure that the geometric scale drops significantly (by a factor of ζ) from one subsample to the next, so there is enough room in each connection between truncated filtrations to kill the noise. • The second one, called the oscillating Rips zigzag, consists in somewhat relaxing the truncation parameter in the Rips filtrations before connecting them together. The effect is to leave enough room in every truncated filtration for the noise to be killed. • The third one, called the image Rips zigzag, consists in taking a nested pair of Morozov zigzags with different filtration parameters, and in connecting them by canonical inclusions to obtain an image zigzag module at the homology level. Taking a pair of zigzags instead of a single zigzag kills the noise in the same way as taking a pair of Rips complexes instead of a single Rips complex did in [START_REF] Chazal | Towards persistence-based reconstruction in Euclidean spaces[END_REF]. Each of these variants comes with the desired guarantee that the sweet and sweeter ranges are equal, meaning that there is guaranteed to be only ephemeral noise in the sweet range even when the underlying space X merely has positive weak feature size. Thus, Rips zigzags offer improvements in both size complexity and signal-to-noise ratio compared to the Rips filtration. The price to pay compared to the basic Morozov zigzag is a somewhat increased time or space complexity (Theorems 5.1 and 5.2). The overhead depends on the variant considered but it always remains bounded, so the variants are tractable alternatives to the Morozov zigzag.

To prove the aforementioned results, we develop new techniques for manipulating zigzag modules and comparing their persistence barcodes:

• We show how arrows in a zigzag module can be reversed while preserving the persistence barcode (Theorem 3.1). • We give a method for removing spaces from a zigzag module while tracking the intervals in its barcode (Theorem 3.2). These low-level manipulations make it possible to transform one module into another while controlling the changes in its barcode, a strategy at the core of the proofs of our main theorems.

Related work.

A different approach to the problem of building sparse filtrations for offsets of point clouds in Euclidean space was presented by Hudson et al. [START_REF] Hudson | Topological inference via meshing[END_REF]. They used ideas from Delaunay refinement mesh generation to build linear size filtrations that provide provably good approximations to the persistence diagram of the offsets. However, that approach requires building a complex that covers the ambient space and includes simplices up to its dimension. Moreover, the construction requires the use of high degree predicates. In contrast, the new methods described here only depend on an intrinsic dimension of data and can be built using only distances comparisons.

Recently, Sheehy [START_REF] Sheehy | Linear-size approximations to the vietoris-rips filtration[END_REF] proposed a method for building a sparse zigzag filtration whose barcode is provably close to that of the Rips filtration as well as a non-zigzagging variant achieving similar guarantees. Also, Dey et al. gave an alternative persistence algorithm for simplicial maps rather than inclusions, which is closely related to zigzag persistence [START_REF] Tamal | Computing topological persistence for simplicial maps[END_REF]. Their approach, when applied to Rips filtrations, similarly gives barcodes that are provably close to that of the Rips filtration. We obtain comparable space/time bounds to these results but get stronger guarantees regarding noise. Methods that approximate the Rips filtration directly can, in principle, have noise that is as large as the noise in the Rips filtration itself (or worse).

BACKGROUND

We use singular homology with coefficients in a fieldomitted in our notations. The homology functor is denoted by H. See [START_REF] Hatcher | Algebraic Topology[END_REF] for an introduction to homology theory.

Distance Functions

The geometric part of our analysis takes place in Euclidean space R d , where • denotes the Euclidean norm. The distance from a point y to a set X ⊂ R d is d(y, X) = infx∈X x-y . When X is compact, the infimum becomes a minimum and we let dX denote the function distance to X. The λ-offset of X, noted X λ , is the sublevel set d -1

X ([0, λ]).

Given another compact set Y ⊂ R d , the Hausdorff distance dH (X, Y ) is the smallest λ such that Y ⊆ X λ and X ⊆ Y λ .
Although dX may not be differentiable everywhere in R d , its gradient ∇X can be extended to be well-defined over all R d [START_REF] Chazal | Geometric inference[END_REF]. A critical point of dX is a point p ∈ R d \ X such that ∇X (p) = 0. A critical value of dX is a value r ≥ 0 such that r = dX (p) for some critical point p. The weak feature size of X, noted wfs(X), is the smallest positive critical value of dX .

Given a finite set P ⊂ R d and a parameter α ≥ 0, the Čech complex Cα(P ) is the nerve of the collection of Euclidean balls of same radius α centered at the points of P . Given another parameter β ≥ α, the canonical inclusion Cα(P ) ⊆ C β (P ) induces a homomorphism HCα(P ) → HC β (P ) at the homology level. Its image is denoted by HC β α (P ) hereafter. The following theorem follows directly from the sampling theory for compact sets in R d [START_REF] Chazal | Geometric inference[END_REF], so its proof is omitted.

Theorem 2.1. (see [START_REF] Oudot | Zigzag Zoology: Rips Zigzags for Homology Inference[END_REF] for the proof ) Let P ⊆ Q be finite sets in R d , such that dH (P, X) < ε and dH (Q, X) < ε for some compact set X ⊂ R d . • If ε < 1 4 wfs(X), then for any α, β ∈ [ε, wfs(X)ε] such that βα ≥ 2ε, for any λ ∈ (0, wfs(X)), the spaces HC β α (P ) and HX λ are isomorphic.

• If ε < 1 6 wfs(X), then for all α, α ′ , β, β ′ ∈ [3ε, wfs(X) -ε] such that β -α ≥ 2ε, β ′ -α ′ ≥ 2ε, α ′ ≥ α and β ′ ≥ β, the homomorphism HC β α (P ) → HC β ′ α ′ (Q)
induced by the following commutative diagram (where the maps are induced by inclusions) is an isomorphism.

HC β (P ) → HC β ′ (Q) ↑ ↑ HCα(P ) → HC α ′ (Q)

Zigzag persistence

The algebraic part of our analysis relies on zigzag persistence theory. We use the terminology introduced by Carlsson and de Silva [START_REF] Carlsson | Zigzag persistence[END_REF]. A zigzag module V is a finite diagram of finite-dimensional vector spaces over a fixed field k:

V = V1 v 1 ←→ V2 v 2 ←→ • • • v n-1 ←→ Vn,
where the notation Vi v i ←→ Vi+1 indicates that the linear map vi can be oriented either forwards (vi : Vi → Vi+1) or backwards (vi : Vi ← Vi+1). An equivalent notation is vi : Vi ↔ Vi+1. The sequence of map orientations defines the type of the module V. A persistence module, as defined in the standard (non-zigzag) persistence literature [START_REF] Zomorodian | Computing persistent homology[END_REF], is a zigzag module in which all the maps are oriented forwards. Thus, all persistence modules of the same length have the same type.

A homomorphism Φ between two zigzag modules V and W of the same type, noted Φ : V → W, is a collection of linear maps φi : Vi → Wi such that the following diagram commutes for all i = 1, • • • , n -1:

Vi φ i o o v i / / Vi+1 φ i+1 Wi o o w i / / Wi+1 (1)
Φ is called an isomorphism if every map φi : Vi → Wi is an isomorphism.

A submodule W of a zigzag module V is defined by subspaces Wi ⊆ Vi such that for all i we have vi(Wi) ⊆ Wi+1 if vi : Vi ↔ Vi+1 is a forward map and vi(Wi+1) ⊆ Wi if vi is a backward map. The maps in W are the restrictions of the maps in V to the Wis, so V and W have the same type. W is called a summand of V if there exists another submodule X of V such that Vi = Wi ⊕ Xi for all i. In that case, we say that V is the direct sum of W and X, written V = W ⊕ X. As pointed out in [START_REF] Carlsson | Zigzag persistence[END_REF], all summands are submodules but not all submodules are summands.

A zigzag module V is called indecomposable if it admits no other summand than the zero module and itself. It has been known since Gabriel [START_REF] Gabriel | Unzerlegbare Darstellungen I[END_REF] that the indecomposable zigzag modules are the so-called interval modules. Given a module type τ and an integer interval [b, d], the interval τ -module with birth time b and death time d is written Iτ [b, d] and defined with spaces Ii such that Ii = k if i ∈ [b, d] and Ii = 0 otherwise, and with identity maps between adjacent copies of the base field k and zero maps elsewhere (the maps are oriented according to τ ). A consequence of Gabriel's result is that every τ -module is isomorphic to a direct sum of finitely many τ -intervals. Moreover, the Krull-Schmidt principle guarantees that this decomposition is unique up to a reordering of the terms-see Proposition 2.2 in [START_REF] Carlsson | Zigzag persistence[END_REF].

Theorem 2.2 (Interval Decomposition).

For every τ -module V there exists a unique finite multiset of interval modules {Iτ [bi, di]} and an isomorphism

Φ : V → i Iτ [bi, di].
Thus, as in standard (non-zigzag) persistence theory, the structure of V is fully and uniquely described by the multiset of integer intervals 1 {[bi, di]}, called the persistence barcode of V and denoted Pers(V).

Carlsson and de Silva gave a constructive proof of Theorem 2.2-see [START_REF] Carlsson | Zigzag persistence[END_REF]Thm. 4.1], which led to an algorithm for computing the decomposition of a zigzag module. Among the concepts and results presented in their paper, the following one plays an important part here.

Given a zigzag module

V = V1 v 1 ←→ • • • v n-1 ←→ Vn and two integers p ≤ q ∈ [1, n], let V[p, q] denote the restriction of V to the index set [p, q]. That is, V[p, q] = Vp vp ←→ • • • v q-1 ←→ Vq. Let Pers(V)| [p,q] denote the restriction of Pers(V) to [p, q]: Pers(V)| [p,q] = {[b, d] ∩ [p, q] | [b, d] ∈ Pers(V)}.
The Restriction Principle [2, Prop. 2.12] connects the two types of restrictions together:

Theorem 2.3 (Restriction). Pers(V[p, q]) = Pers(V)| [p,q] .

MANIPULATING ZIGZAG MODULES

Arrow reversal

Suppose we have a zigzag module V = V1 ↔ • • • ↔ Vn and we want to reverse the map V k ↔ V k+1 for some arbitrary index k in the range [1, n -1], while preserving the persistence barcode of V. The following theorem states that this is always possible, moreover with a reverse map that is closely related to the original map.

Theorem 3.1 (Arrow Reversal). Let V = V1 ↔ • • • ↔ V k f ←→ V k+1 ↔ • • • ↔ Vn be a zigzag module. Then, there is a map g : V k ↔ V k+1 oriented oppo- site to f , such that f • g| im f = 1 im f and g • f |im g = 1im g , and the zigzag module V * obtained from V by replacing the submodule V k f ←→ V k+1 by V k g ←→ V k+1 has the same per- sistence barcode as V.
When f is injective, g is surjective and g • f is the identity over the domain of f . Conversely, when f is surjective, g is injective and f • g is the identity over the codomain of f . These properties are useful when V is part of a commutative diagram, as they sometimes imply commutativity is preserved after the arrow reversal-as in the proof of Theorem 4.1.

Proof. Let Φ : V → i Iτ [bi, di] be the decomposition of V given by the Interval Decomposition Theorem 2.2. Denote the spaces of Iτ [bi, di] by I i 1 . . . , I i n , and let Φ = (φ1, . . . , φn) where each φj : Vj → i I i j is an isomorphism. We assume without loss of generality that f is oriented forwards; the case when f is oriented backwards is symmetric. The map

f ′ = i (I i k → I i k+1 ) makes the following diagram commute. V k φ k f / / V k+1 φ k+1 i I i k f ′ / / i I i k+1 (2)
To reverse f , we first reverse each map I i k → I i k+1 separately. Recall that this map is either the identity or zero. The reversal of an identity map is an identity map and the reversal of a zero map is a zero map, so we take

I i k 1 ← I i k+1 if I i k 1 → I i k+1 and I i k 0 ← I i k+1 if I i k 0 → I i k+1 . We thus get a new interval module Iτ * [bi, di], where τ * is the zigzag type of V * . Let now g ′ = i (I i k ← I i k+1
) be the reverse of f ′ , and let g = φ -1 k • g ′ • φ k+1 be the reverse of f , which gives the following commutative diagram.

V k φ k V k+1 g o o φ k+1 i I i k i I i k+1 g ′ o o (3)
The commutativity of this diagram and the definition of Φ imply that the the same isomorphisms φj induce an isomorphism Φ * : V * → i Iτ * [bi, di]. So, Φ * yields an explicit interval decomposition of V * with the same intervals as that of V, and thus Pers(V * ) = Pers(V).

It only remains to prove that f • g| im f = 1 im f and g • f |im g = 1im g . By symmetry, it suffices to prove just one of these, say f • g| im f = 1 im f . Using the isomorphism Φ, it suffices to prove the equivalent statement for f ′ and g ′ . The statement now follows from the observation that

f ′ • g ′ = i ((I i k → I i k+1 ) • (I i k ← I i k+1 )) = 1 im f ′ ⊕ 0 ker g ′ ,
where the final step separates the identity maps and the zero maps into two groups.

Space Removal

Suppose we want to remove the space V k from a zigzag module V = V1 ↔ • • • ↔ Vn while preserving most of the persistence barcode. The following theorem shows that this is always possible, and that every interval

[b, d] in Pers(V) becomes [b, d] \ {k} in the new index set {1, • • • , k -1, k + 1, • • • , n}. Note that this is still a single interval, denoted by [b, d] k for clarity. Theorem 3.2 (Space Removal). Let V be a zigzag module containing V k-1 f ←→ V k g ←→ V k+1 .
There exists a map h : V k-1 ↔ V k+1 such that the zigzag module V * formed by removing V k , f , and g from V and replacing them with h has barcode

Pers(V * ) = {[b, d] k | [b, d] ∈ Pers(V)} .
Furthermore, if any map V k-1 ↔ V k+1 commutes with f and g, then so does h.

The first step towards the proof of this theorem is the following special case for removing a space between maps oriented in the same direction.

Lemma 3.3 (Composition)

.

Given V = V1 ↔ • • • ↔ V k-1 f -→ V k g -→ V k+1 ↔ • • • ↔ Vn, let V * = V1 ↔ • • • ↔ V k-1 g•f -→ V k+1 ↔ • • • ↔ Vn. Then, Pers(V * ) = {[b, d] k | [b, d] ∈ Pers(V)}. Proof. Let Φ : V → i Iτ [bi, di]
be the decomposition of V given by the Interval Decomposition Theorem 2.2. Denoting the spaces of Iτ [bi, di] by I i 1 . . . , I i n , and letting

f ′ = i (I i k-1 → I i k ), g ′ = i (I i k → I i k+1
), and Φ = (φ1, . . . , φn) where each φj : Vj → i I i j is an isomorphism, we have the following commutative diagram.

V k-1 f / / φ k-1 V k g / / φ k V k+1 φ k+1 i I i k-1 f ′ / / i I i k g ′ / / i I i k+1 (4)
By removing the kth space from each interval module Iτ [bi, di], we get a new interval decomposition i Iτ * [bi, di] k , where τ * is the type of V * . Observe that the commutativity of (4) implies that the following diagram commutes as well.

V k-1 g•f / / φ k-1 V k+1 φ k+1 i I i k-1 g ′ •f ′ / / i I i k+1 (5)
Thus, Φ * = (φ1, . . . , φ k-1 , φ k+1 , . . . , φn) is an isomorphism from

V * to i Iτ * [bi, di] k . It follows that Pers(V * ) = Pers( i Iτ * [bi, di] k ) = {[b, d] k | [b, d] ∈ Pers(V)}
as claimed in the lemma.

We now combine this special case with the Arrow Reversal Theorem 3.1 to prove the general Space Removal Theorem.

Proof of Theorem 3.2. There are four cases to consider depending on the orientations of the maps f and g. When there exists a map h ′ : V k-1 ↔ V k+1 that commutes with f and g, then either

f = g • h ′ or f = h ′ • g or g = f • h ′ or g = h ′ • f.
The proofs of these four cases are similar, so we only give the proof of the first case here, leaving the others as an exercise.

Assuming that f = g • h ′ , we have im f ⊆ im g. So, applying the Arrow Reversal Theorem 3.1 on the map g results in a map g ′ such that g • g ′ | im f = 1 im f . Then, the Composition Lemma 3.3 lets h = g ′ • f , so we obtain

g • h = g • g ′ • f = g • g ′ | im f • f = 1 im f • f = f as claimed.

RIPS ZIGZAGS

Let P be a finite point cloud in some metric space, and suppose that the matrix of pairwise distances between the points of P is known. Given any ordering (p1, • • • , pn) on the points of P , let Pi := {p1, . . . , pi} denote the ith prefix, and define the ith geometric scale as εi = dH (Pi, P ). Since Pi grows as i increases, we have ε1

≥ ε2 ≥ • • • ≥ εn = 0.
Given a choice of multipliers η ≤ ρ, Chazal and Oudot [START_REF] Chazal | Towards persistence-based reconstruction in Euclidean spaces[END_REF] proposed to do homological inference from P using the sequence of short filtrations Rηε i (Pi) ֒→ Rρε i (Pi). Zigzag persistence makes it possible to replace this sequence of short filtrations by a single zigzag filtration, a representative por-tion of which is depicted below.

Rρε i-1 (Pi) Rρε i (Pi+1) Rηε i-1 (Pi-1) A A Ò Ò Ò Ò Ò Ò Ò Rηε i (Pi) ] ] ```````A A Ò Ò Ò Ò Ò Ò Ò Rηε i+1 (Pi+1)
] ] ```````(

The zigzag module induced at the homology level by this diagram is referred to as the oscillating Rips zigzag (oR-ZZ for short) hereafter.

Analysis of the oscillating Rips zigzag

The following result gives conditions on η and ρ for the persistence barcode of the oR-ZZ to exhibit the homology of the shape underlying an input point cloud P ⊂ R d . The proof relies on the fact that Čech and Rips complexes are interleaved as follows 2 in R d , where

ϑ d = d 2(d+1) ∈ [ 1 2 , 1 √ 2 ]. ∀α ≥ 0, C α 2 (P ) ⊆ Rα(P ) ⊆ C ϑ d α (P ). ( 7
)
Theorem 4.1. Let ρ and η be multipliers such that ρ > 10 and 3 ϑ d < η < ρ-4 2ϑ d . Let X ⊂ R d be a compact set and let P ⊂ R d be such that dH (P, X) < ε with

ε < min ϑ d η-3 6ϑ d η , η-3/ϑ d 3ρ+η , ρ-2ϑ d η-4 6(ρ-2ϑ d η) , ρ-2ϑ d η-4 (4ϑ d +1)ρ-2ϑ d η wfs(X).
Then, for any l > k such that

ε l ≥ max 3ε ϑ d η -3 , 4ε ρ -2ϑ d η -4
, and

ε k < min 1 6 wfs(X) -ε, 1 ϑ d ρ + 1 (wfs(X) -ε) ,
the oR-ZZ restricted to HRρε k (P k+1 ) ← • • • ← HRηε l (P l ) has a persistence barcode made only of full-length intervals and ephemeral (length zero) intervals, the number of fulllength intervals being equal to the dimension of HX λ for any λ ∈ (0, wfs(X)).

Proof. Let ρ = ρ 2 and η = ϑ d η. Our hypotheses imply ρ 2 ≥ ϑ d η, so using [START_REF] Tamal | Computing topological persistence for simplicial maps[END_REF] we can factor the inclusion maps in [START_REF] De | Coverage in sensor networks via persistent homology[END_REF] through Čech complexes with multipliers η and ρ as follows.

Rρε i (Pi+1)

Cρε i-1 (Pi) Cρε i (Pi) o o = = | | | | | | | | / / Cρε i (Pi+1) a a f f f f f f f f Cηε i-1 (Pi) O O Cηε i (Pi) O O o o / / Cηε i (Pi+1) O O Rηε i (Pi) a a f f f f f f f f = = |
HC ρεr ηεr (Ps)-of the homomorphisms induced at the homology level by the vertical arrows in the above diagram.

HRρε i (Pi+1)

HC ρε i-1 ηε i-1 (Pi) HC ρε i ηε i (Pi) o o > > } } } } } } } } / / HC ρε i ηε i (Pi+1)
èe e e e e e e

ρε k ηε k (P k ) → • • • ← HC ρε l
ηε l (P l ), and let V be the mixed zigzag

HC ρε k ηε k (P k ) → HRρε k (P k+1 ) ← HC ρε k ηε k (P k+1 ) ← HRηε k+1 (P k+1 ) → HC ρε k+1 ηε k+1 (P k+1 ) → • • • → HC ρε l-1 ηε l-1 (P l-1 ) → HRρε l-1 (P l ) ← HC ρε l-1
ηε l-1 (P l ) ← HRηε l (P l ) → HC ρε l ηε l (P l ). Our goal is to relate Pers(V) to both Pers(U) and Pers(W), which we will do by turning V successively into U and W via arrow reversals and space removals while tracking the changes in its persistence barcode. We will admit the following easy consequences of Theorem 2.1 under our geometric hypothesessee [START_REF] Oudot | Zigzag Zoology: Rips Zigzags for Homology Inference[END_REF] for a proof:

(i) all spaces in W are isomorphic to HX λ , and all maps in W are isomorphisms, (ii) the map HRρε

i (Pi+1) ← HC ρε i ηε i (Pi+1) is injective for any i ∈ [k, l -1], (iii) the map HRηε i (Pi) → HC ρε i ηε i (Pi) is surjective for any i ∈ [k, l].
To turn V into W, we first apply Theorem 3.1 on every injective map HRρε i (Pi+1) ← HC ρε i ηε i (Pi+1) and on every surjective map HRηε i (Pi) → HC

ηε i-1 (Pi) HC ρε i ηε i (Pi) }} } } } } } } o o > > } } } } } } } } / / HC ρε i ηε i (Pi+1) HRηε i (Pi)
èe e e e e e e Now, we remove the Rips complexes from V * by composing all the adjacent maps with same orientation. Since composition preserves commutativity of the subdiagrams, the following diagram involving W (straight path) and the newly obtained zigzag W * (curved path) commutes:

HC ρε i-1 ηε i-1 (Pi) HC ρε i ηε i (Pi) o o _ _ / / HC ρε i ηε i (Pi+1)
Hence, the zigzags W and W * are identical. It suffices now to compare Pers(V * ) to Pers(W * ). Recall that W * is obtained from V * by removing the Rips complexes. In the process, every interval of Pers(V * ) either vanishes or turns into some interval of Pers(W * ). By the Space Removal Theorem 3.2, only the ephemeral intervals of Pers(V * ) may vanish because there are no consecutive Rips complexes in V * . Moreover, the full-length intervals of Pers(V * ) are mapped bijectively to those of Pers(W * ) because V * and W * have the same endpoints. It follows then from (i) that Pers(V * )-and thus Pers(V)-contains only full-length intervals of multiplicity dim HX λ and possibly some ephemeral intervals.

We turn V into U by removing the Čech complexes. First, we restrict V to HRρε k (P k+1 ) ← HC

ρε k ηε k (P k+1 ) ← • • • ← HC ρε l-1
ηε l-1 (P l ) ← HRηε l (P l ), thus removing the Čech complexes at either ends of the zigzag. Since k < l, the Restriction Theorem 2.3 tells us that the full-length intervals in the barcode of the shortened zigzag V * are in bijection with the ones in the barcode of V, while the other intervals in Pers(V * ) come from length-zero intervals in Pers(V) and cannot be longer. We then compose the incoming and outgoing maps at Čech complexes in the sequence to obtain U. By the Space Removal Theorem 3.2, only the intervals starting or ending at a Čech complex can be affected by this operation, and these can only be shortened. Therefore, the number of full-length intervals remains the same as in the barcode of V * , and the other intervals remain ephemeral.

Morozov zigzag

The following limit case of the oscillating Rips zigzag where parameters η, ρ are equal has been integrated into the Dionysus library [START_REF]By Dmitriy Morozov[END_REF] since early 2009.

• • • ← Rρε i (Pi) → Rρε i (Pi+1) ← Rρε i+1 (Pi+1) → • • • (8)
We call Morozov zigzag (M-ZZ for short) the zigzag module induced at the homology level by this diagram. The motivation for letting η = ρ is obvious from a computational point of view: the closer η to ρ, the fewer simplex additions and deletions during the zigzag calculation. Moreover, as reported by its author [START_REF] Morozov | [END_REF], the M-ZZ has given surprisingly good results in preliminary experiments, despite the fact that η = ρ clearly violates the hypotheses of Theorem 4.1. We are able to provide the following weaker guarantee. Theorem 4.2. Let ρ > 10 be a multiplier. Let X ⊂ R d be a compact set and let P ⊂ R d be such that dH (P, X) < ε with ε < ρ-10 (3+10ϑ d )ρ wfs(X). Then, for any l > k such that

ε l ≥ 10ε ρ -10
, and

ε k < min 1 6 wfs(X) -ε, 5 (1 + 5ϑ d )ρ + 5 (wfs(X) -ε) , the M-ZZ restricted to HRρε k (P k+1 ) ← • • • ← HRρε l (P l )
has a number of full-length intervals that is at least the dimension of HX λ for any λ ∈ (0, wfs(X)).

Proof. Let V denote the restriction of the M-ZZ to

HRρε k (P k+1 ) ← • • • ← HRρε l (P l ).
Consider the image Čech zigzags (as introduced in the proof of Theorem 4.1) of parameters (η1, ρ1) and (η2, ρ2) respectively, where η1

= ρ 2 -2(1 + ε ε l ), ρ1 = ρ 2 , η2 = ϑ d ρ, ρ2 = ϑ d ρ + 2(1 + ε ε l )
, and call U and W respectively their restrictions to the same index set as V.

Since η1 ≤ ρ1 = ρ 2 ≤ ϑ d ρ = η2 ≤ ρ2
, the canonical inclusions between Čech complexes induce homomorphisms between the spaces of U and W of same index. This family of homomorphisms forms a homomorphism U → W. By [START_REF] Tamal | Computing topological persistence for simplicial maps[END_REF], every inclusion Cρ

1 ε i (Q) ֒→ Cη 2 ε i (Q) factors through Rρε i (Q), so U → W itself factors through V. Let U Φ -→ V Ψ
-→ W be the factorization. Under our geometric hypotheses, it follows from Theorem 2.1 that Ψ • Φ is an isomorphism. Hence, we have V = im Φ ⊕ ker Ψ, where im Φ is the submodule of V formed by the images of the homomorphisms comprising Φ, and ker Ψ is the submodule of V formed by the kernels of the homomorphisms comprising Ψ. Hence, by the Interval Decomposition Theorem 2.2 (more precisely the uniqueness part of it), we have Pers(im Φ) ⊆ Pers(V). Moreover, since the homomorphisms comprising Φ are injective, Φ is an isomorphism onto its image, so Pers(im Φ) = Pers(U). The conclusion of the theorem follows because the spaces in U are all isomorphic to HX λ and the maps in U are all isomorphisms, as guaranteed by Theorem 2.1.

Thus, the topological signal of X persists in the M-ZZ throughout a sweet range of the form [O(ε), Ω(wfs(X))]. The question of the resilience of the topological noise within this range remains open, and there currently is no theoretical evidence that the noise should not persist, even under the assumption that X has positive weak feature size.

We now further assume that X has positive µ-reach, denoted rchµ(X) > 0, for some sufficiently large µ. Recall that rchµ(X) is the infimum of distances from X to points outside of X where the gradient of the distance to X is less than µ [START_REF] Chazal | Geometric inference[END_REF]. Attali et al. showed that if dH (P, X) < ε with ε sufficiently small compared to rchµ(X) and µ sufficiently large, then for some values of α, the Rips complex Rα(P ) is homotopy equivalent to X λ for λ ∈ (0, rchµ(X))-see [START_REF] Attali | Vietoris-Rips complexes also provide topologically correct reconstructions of sampled shapes[END_REF]Theorem 14]. The immediate consequence of their result is that for a multiplier ρ and an index i, Rρε i (Pi) and Rρε i (Pi+1) are both homotopy equivalent to X λ for λ ∈ (0, rchµ(X)) whenever

µ(2 -µ)(2ρεi -2ϑ d εi -2(εi + ε)) 1 + µ(1 -µ) -1 -µ(2 -µ) (2ϑ d ρ+1)ε i +ε rchµ(X) 2 > rchµ(X).
This condition depends on rchµ(X), its parameter µ, the multiplier ρ, the Hausdorff distance of the sample ε, and εi. Attali et al. showed that there do exist values for which the condition is satisfied. We do not derive the space of valid assignment of constants here, but merely note that this result implies that there is a multiplier ρ and a range of scales of the form [O(ε), Ω(rchµ(X))] for which the M-ZZ exhibits no noise at all. This holds because Theorem 4.2 implies that the signal is present in the sweet range and the Attali et al. result shows that every space in the strictly smaller range has the same homology as X λ . We call this the sweeter range. Note that the quantities hidden in the big-O and big-Ω notations depend on µ, so the sweeter range requires a sufficiently large µ to be non-empty. Moreover, since the upper bound depends on rchµ(X) rather than wfs(X), the sweeter range can be arbitrarily smaller than the sweet range.

Discretized Morozov Zigzag

We now describe a discretization scheme for the Morozov zigzag that ensures that the noise gets killed within the sweet range. Given a map ζ : R>0 → R>0, referred to as the scale drop function hereafter, we select a subset of the indices 1, • • • , n in the ordered point cloud P = {p1, • • • , pn} by the following iterative procedure: let n1 = 1, and

∀i ≥ 1, let ni+1 = min {j > ni | εj ≤ ζ(εn i ) • εn i } .
Note that nr = n for some index r since εn = 0 < εn-1. We then build the following discretized version of the zigzag diagram of (8):

• • • ← Rρε n i (Pn i ) → Rρε n i (Pn i+1 ) ← Rρε n i+1 (Pn i+1 ) → • • •
where Pn 1 = P1 = {p1} and Pn r = Pn = P . The zigzag induced at the homology level is called the discretized Morozov zigzag (dM-ZZ) hereafter.

Theorem 4.3. Given a choice of multiplier ρ > 10, suppose P ⊂ R d and there is some compact set X ⊂ R d such that dH (P, X) < ε with ε < ρ-10 (3+10ϑ d )ρ wfs(X). Then, for any choice of scale drop function ζ that satisfies

∀i ∈ [1, r -1], ζ(εn i ) ≤ ρ 2ϑ d ρ + 4 - 2 ϑ d ρ + 2 ε εn i , (9)
for any n l > n k such that

εn l ≥ 10ε ρ -10
, and

εn k < min 1 6 wfs(X) -ε, 5 (1 + 5ϑ d )ρ + 5 (wfs(X) -ε) , the dM-ZZ restricted to HRρε n k (Pn k+1 ) ← • • • ← HRρε n l (Pn l )
has a barcode with only three classes of intervals:

• full-length intervals, whose number is equal to the dimension of HX λ for any λ ∈ (0, wfs(X)), • ephemeral (length zero) intervals, • intervals of the form [HRρε n i (Pn i ), HRρε n i (Pn i+1 )],

which are ephemeral on the scale of the geometric scales.

Proof. The proof of Theorem 4.2 applies verbatim to the dM-ZZ, with indices 1, 2, • • • , n replaced by n1, n2, • • • , nr. Thus, the restriction V of the dM-ZZ to HRρε n k (Pn k+1 ) ← • • • ← HRρε n l (Pn l ) has a persistence barcode with at least dim(HX λ ) full-length intervals. Now, hypothesis [START_REF] Edelsbrunner | Topological persistence and simplification[END_REF] implies that for all i we have

εn i ≥ (2ϑ d + 4 ρ )εn i+1 + 4 ρ ε, therefore the inclusion Rρε n i+1 (Pn i+1 ) → Rρε n i (Pn i+1 ) fac- tors through C ϑ d ρεn i+1 (Pn i+1 ) → C ρ 2 εn i (Pn i+1
). Furthermore, under our geometric hypotheses, it follows from Theorem 2.1 that

rank HC ϑ d ρεn i+1 (Pn i+1 ) → HC ρ 2 εn i (Pn i+1 ) = dim(HX λ ),
therefore by composition we have rank HRρε n i (Pn i+1 ) ← HRρε n i+1 (Pn i+1 ) ≤ dim(HX λ ).

Intuitively, this means that only the signal can go through the link HRρε n i (Pn i+1 ) ← HRρε n i+1 (Pn i+1 ), and that the noise gets killed. More formally, by the Restriction Theorem 2.3 the total multiplicity of the intervals spanning [HRρε n i (Pn i+1 ), HRρε n i+1 (Pn i+1 )] in Pers(V) is at most dim(HX λ ). It follows that among these intervals only the full-length one has non-zero multiplicity. Thus, Pers(V) contains three types of intervals: full-length ones, ephemeral ones, and ones of the form [HRρε n i (Pn i ), HRρε n i (Pn i+1 )]. These are not ephemeral, however they become so once represented on the scale of the geometric scales.

Note that ε usually remains unknown in practice, so the user cannot merely set ζ(εn i ) to be the quantity

ρ 2ϑ d ρ+4 - 2 ϑ d ρ+2 ε εn i
. The bounds given in Theorem 4.3 suggest to let ζ be the constant map

ζ = 3ρ + 20 10(ϑ d ρ + 2) , (10)
so that ζ(εn i ) satisfies hypothesis (9) as long as εn i ≥ 10ε ρ-10 . Thus, the conclusion of the theorem continues to hold within the same sweet range. In fact, any smaller value could be chosen for ζ(εn i ) without affecting the sweet range. Nevertheless, the larger ζ(εn i ) the better in general, since the smaller the gaps on the geometric scale the more chances there are that the set of discretization values εn i intersects the sweet range, and furthermore the smaller the gaps the smaller the complexes involved in the construction of the discretized Morozov zigzag.

Image Rips zigzag

We end this section with another Rips zigzag construction, consisting of a nested pair of Morozov zigzags with multipliers ρ ≥ η ≥ 0. Canonical inclusions between Rips complexes give the following commutative diagram at the homology level.

• • • ← HRρε i (Pi) → HRρε i (Pi+1) ← HRρε i+1 (Pi+1) → • • • ↑ ↑ ↑ • • • ← HRηε i (Pi) → HRηε i (Pi+1)← HRηε i+1 (Pi+1) → • • •
The vertical arrows in this diagram form a homomorphism from the M-ZZ with multiplier η to the M-ZZ with multiplier ρ. The image Rips zigzag (iR-ZZ) is defined as the image of this homomorphism. More precisely, it is written as follows, where each HR ρεr ηεr (Ps) denotes the image of the map HRηε r (Ps) → HRρε r (Ps).

• • • ← R ρε i ηε i (Pi) → R ρε i ηε i (Pi+1) ← R ρε i+1
ηε i+1 (Pi+1) → • • • Image Rips zigzags have been available in the Dionysus library [START_REF]By Dmitriy Morozov[END_REF] since early 2009, with no theoretical guarantee on their behavior. Here we provide a guarantee on the output that is similar to (and even slightly better than) the one obtained in Theorem 4.1 for the oscillating Rips zigzag. Theorem 4.4. Given a choice of multipliers ρ, η such that ρ > 10 and 3 ϑ d < η < ρ-4 2ϑ d , suppose P ⊂ R d and there is some compact set X ⊂ R d such that dH (P, X) < ε with

ε < min ϑ d η-3 6ϑ d η , η-3/ϑ d 3ρ+η , ρ-2ϑ d η-4 6(ρ-2ϑ d η) , ρ-2ϑ d η-4 (4ϑ d +1)ρ-2ϑ d η wfs(X).
Then, for any l > k such that

ε l ≥ max 3ε ϑ d η -3 , 4ε ρ -2ϑ d η -4
, and

ε k < min 1 6 wfs(X) -ε, 1 ϑ d ρ + 1 (wfs(X) -ε) ,
the iR-ZZ restricted to HR

ρε k ηε k (P k+1 ) ← • • • ← HR ρε l
ηε l (P l ) contains only isomorphisms, and its spaces are isomorphic to HX λ for any λ ∈ (0, wfs(X)). Therefore, its persistence barcode is made only of full-length intervals, whose number equals the dimension of HX λ .

Proof. Our hypotheses imply ρ 2 ≥ ϑ d η, so we can use (7) to obtain the following diagram where all arrows are canonical inclusions.

• • • ← C ϑ d ρε i (Pi) → C ϑ d ρε i (Pi+1) ← C ϑ d ρε i+1 (Pi+1) → • • • ↑ ↑ ↑ • • • ← Rρε i (Pi) → Rρε i (Pi+1) ← Rρε i+1 (Pi+1) → • • • ↑ ↑ ↑ • • • ← C ρ 2 ε i (Pi) → C ρ 2 ε i (Pi+1) ← C ρ 2 ε i+1 (Pi+1) → • • • ↑ ↑ ↑ • • • ← C ϑ d ηε i (Pi) → C ϑ d ηε i (Pi+1) ← C ϑ d ηε i+1 (Pi+1) → • • • ↑ ↑ ↑ • • • ← Rηε i (Pi) → Rηε i (Pi+1) ← Rηε i+1 (Pi+1) → • • • ↑ ↑ ↑ • • • ← C η 2 ε i (Pi) → C η 2 ε i (Pi+1) ← C η 2 ε i+1 (Pi+1) → • • • This diagram
induces a homomorphism Φ from the image Rips zigzag U of parameters η, ρ to the image Čech zigzag V of parameters ϑ d η, ϑ d ρ. Call U * (resp. V * ) the restriction of U (resp. V) to the scale range [(ε k , P k+1 ), (ε l , P l )]. It turns out that the restriction of Φ to U * is an isomorphism onto V * . To see this, pick an arbitrary index i in the range [k + 1, l] and consider the following sequence of inclusions between Čech complexes:

C η 2 ε i (Pi) a -→ C ϑ d ηε i (Pi) b -→ C ρ 2 ε i (Pi) c -→ C ϑ d ρε i (Pi+1).
Under our geometric hypotheses, it follows from Theorem 2.1 that the maps c • b • a and b induce homomorphisms of same rank at the homology level. Since a and c factor through Rips complexes as shown above, rank inequalities induced by composition imply that the map HR ρε i ηε i (Pi) → HC

ϑ d ρε i ϑ d ηε i (Pi
) is an isomorphism. The same argument holds for the map HR ρε i ηε i (Pi+1) → HC

ϑ d ρε i ϑ d ηε i (Pi+1), with i ∈ [k, l -1]
. Thus, the restriction of Φ to U * is indeed an isomorphism onto V * . The conclusion of the theorem follows then from the fact that the spaces in V * are isomorphic to HX λ and the maps in V * are isomorphisms, as guaranteed by Theorem 2.1.

COMPLEXITY BOUNDS

In this section we assume that the ordering of the points of P is by furthest point sampling 3 . Then, every prefix Pi is an εi-sparse εi-sample of P , and so standard ball packing arguments like the one used in [START_REF] Guibas | Reconstruction using witness complexes[END_REF]Lemma 4.1] can be applied to bound the memory usage and running times of the various Rips zigzags.

Memory usage. The most relevant parameter is the multiplier ρ, which determines the size of the biggest complex in a Rips zigzag.

Theorem 5.1. Suppose P is a finite metric space of doubling dimension d. Then, for any k ≥ 0, the number of k-simplices in the current complex at any time of the construction of the M-ZZ of parameter ρ is at most 2 O(kd log ρ) n, where n is the cardinality of P . The same bound applies to the oR-ZZ and iR-ZZ of parameter ρ, regardless of the value of η ≤ ρ. Finally, given a scale drop function ζ bounded from below by some quantity ζ0 > 0, the size bound for the dM-ZZ is 2 O(kd log ρ ζ 0

) n. 3 Arbitrary orderings may lead to local oversampling and thus to an uncontrolled local growth of the complexes.

We use the doubling dimension here because it may be smaller than the Euclidean dimension and gives a sharper bound.

Since the theoretical lower bounds on ρ derived in Section 4 are constant (ρ > 10), one is allowed to set ρ to some constant value in practice and benefit from our guarantees on the quality of the output. Meanwhile, Theorem 5.1 ensures that the number of k-simplices in the current complex remains at most 2 O(kd) n throughout the construction of the M-ZZ, oR-ZZ or iR-ZZ. When using the dM-ZZ, one can also set ζ to be a constant map equal to some constant value as in [START_REF] Gabriel | Unzerlegbare Darstellungen I[END_REF], thus benefiting from the theoretical guarantees on the quality of the output while maintaining the number of ksimplices in the current complex below 2 O(kd) n throughout the construction of the zigzag. Note however that the exact complex size is bigger than the one achieved with the other types of Rips zigzags when ζ < 1. These asymptotic bounds are as good in order of magnitude as the ones achieved with previous lightweight structures [START_REF] Tamal | Computing topological persistence for simplicial maps[END_REF][START_REF] Sheehy | Linear-size approximations to the vietoris-rips filtration[END_REF].

Running Time. In experiments, we observed a significant slowdown in the running time of the oR-ZZ compared to the M-ZZ, iR-ZZ and dM-ZZ. This is a result of the oR-ZZ inserting and removing the same simplices many times. Thus, the total number of insertions is a relevant indicator of running time.

Theorem 5.2. Suppose P is sitting in some metric space of doubling dimension d. Then, for any k ≥ 0, the total number of k-simplices inserted in the construction of the M-ZZ of parameter ρ is at most 2 O(kd log ρ) n, where n is the cardinality of P . The same bound applies to the iR-ZZ of parameter ρ, regardless of the value of η ≤ ρ. For the oR-ZZ of parameters η ≤ ρ, the bound becomes 2 O(kd log ρ) n 2 . Finally, given a scale drop function ζ bounded from below by some quantity ζ0 > 0, the bound for the dM-ZZ is 2 O(kd log ρ ζ 0

) n.

EXPERIMENTS

The Rips package of the C++ library Dionysus [START_REF]By Dmitriy Morozov[END_REF] provides efficient implementations of the M-ZZ and iR-ZZ. Our implementations of the dM-ZZ and oR-ZZ are built around this package and will be integrated into it in the near future.

As a proof of concept, we ran our code on the so-called Clifford data set from [START_REF] Hudson | Topological inference via meshing[END_REF], which was obtained by evenly spacing 2, 000 points along the line l : y = 20x mod 2π in the 2-d flat torus (R mod 2π) 2 , then mapping the points onto the Clifford torus in R 4 via the embedding f : (u, v) → (cos u, sin u, cos v, sin v). This data set admits three nontrivial candidate underlying spaces: at small scales, the image of l through f , which is a closed helicoidal curve on the torus; at larger scales, the torus itself; at even larger scales, the 3-sphere of radius √ 2 on which the torus is sitting. Table 1: Maximum complex size (in number of simplices) and total number of simplex insertions.

We ran the M-ZZ, dM-ZZ, iR-ZZ and oR-ZZ using the parameter values given in Table 1. Although these values lie outside the intervals prescribed by the theory, they were sufficient to obtain good results in practice. The corresponding results for the M-ZZ and dM-ZZ are reported in Figures 1 and2 respectively. The barcodes are represented on a logarithmic geometric scale (i.e. the horizontal axis shows the value of log 2 εi), with ephemeral (length zero) intervals removed for clarity. The results obtained with the iR-ZZ and oR-ZZ are similar to Figure 2 and therefore omitted.

The three spaces underlying the input data (curve, torus, 3-sphere) appear in all these barcodes, meanwhile the topological noise remains small (M-ZZ) or even ephemeral (dM-ZZ, iR-ZZ, oR-ZZ). Of particular interest is the 3-sphere, whose corresponding 3-homology cycle appears only at large scales and for a short while due to the fact that the 3sphere is not densely sampled by the point cloud. This delicate 3-cycle certainly exists in the barcode of the standard Rips filtration, however it cannot be observed in practice. As mentioned in [START_REF] Hudson | Topological inference via meshing[END_REF] the union of balls of same radius α around the data points covers the entire sphere only for α ≥ 4 -2 √ 2, so in view of (7) the corresponding 3homology cycle can appear in the Rips filtration only at a parameter value α ≥ √

4-2 √ 2 ϑ d > √ 2
. Now, a simulation reveals that the 4-skeleton of the Rips complex at such parameter values contains more than 31 billion simplices, a size that lies at least 2 orders of magnitude beyond the sizes currently handled by existing implementations. On a 24-GB machine we were able to store the 4-skeleton of the Rips filtration and compute its persistent homology within the main memory up to α ≈ 0.625 using Dionysus. The corresponding truncated barcode is represented on a log 2 scale in Figure 3. As expected, it shows only the curve and the torus, not the 3-sphere. by Hudson et al. using their mesh-based filtration [START_REF] Hudson | Topological inference via meshing[END_REF]. Although the scale has been adapted to be the same as in the previous figures, any direct comparison of the barcodes should be made with the caveat that the intervals are halfopen as in [START_REF] Hudson | Topological inference via meshing[END_REF], rather than closed as in the present paper. Nevertheless, the general trend here is that although the three spaces underlying the input data do appear in the barcode, the amount of noise is significant and its structure is quite irregular, despite the regularity of the data. Most notably, the amplitude of the noise in 1-or 2-homology is larger than the amplitude of the signal in 3-homology, which tends to obscure the information carried within the barcode from the user's point of view. The superiority of Rips zigzags over mesh-based filtrations is also in terms of efficiency: as reported in [START_REF] Hudson | Topological inference via meshing[END_REF] the mesh-based filtration contains 12 million simplices, whereas our Rips zigzags contain less than 200,000 simplices at any given time-see Table 1. The oR-ZZ did perform over 7 million simplex insertions and deletions in total, however the other zigzags performed much fewer such operations as predicted by Theorem 5.2. In practice, computing any of our Rips zigzags together with its barcode took only a few minutes on a single Intel Xeon CPU core running at 2.40 GHz, whereas computing the mesh-based filtration and its barcode took hours on a similar machine.

CONCLUSION

In this paper, we explored several Rips-like zigzags that achieve both small size and bounded topological noise for homology inference. While our size bounds are comparable to the ones obtained in previous work, our bounds on the noise are clearly better since they claim the noise is either totally absent or just ephemeral.

Our proofs rely on general new techniques for manipulating and comparing zigzag modules. We hope that these techniques will find further use and stimulate new research and applications of zigzag persistence. We refer the interested reader to [START_REF] Oudot | Zigzag Zoology: Rips Zigzags for Homology Inference[END_REF]Section 6] for an in-depth discussion of the potential and limitations of our results. For now we wish to emphasize that the Arrow Reversal Theorem 3.1 may open the door to a new approach to the open question of defining interleavings between zigzag modules and deriving stability results for their persistence barcodes. This approach is based on a reduction to non-zigzag persistence via repeated arrow reversals, and it provides both a notion of interleaving and a corresponding stability result. It remains unclear how canonical this approach is, and we intend to investigate this question further in the near future.

  handled by reversing one of the maps using the Arrow Reversal Theorem 3.1 and then applying the Composition Lemma 3.3. The arrow reversal preserves the persistence diagram and the composition removes the kth index from each interval as desired. It only remains to show that the commutativity property holds.

  Let U be the restriction of the oR-ZZ to HRρε k (P k+1 ) ← • • • ← HRηε l (P l ), let W be the restriction of the image Čech zigzag to HC

  ρε i ηε i (Pi), to get a new zigzagV * = HC ρε k ηε k (P k ) → HRρε k (P k+1 ) → HC ρε k ηε k (P k+1 ) ← HRηε k+1 (P k+1 ) ← HC ρε k+1 ηε k+1 (P k+1 ) → • • • ← HC ρε l-1ηε l-1 (P l-1 ) → HRρε l-1 (P l ) → HC ρε l-1 ηε l-1 (P l ) ← HRηε l (P l ) ← HC ρε l ηε l (P l ) that has the same persistence barcode as V. Moreover, the reverse maps provided by Theorem 3.1 make the triangles commute in the resulting diagram interleaving V * and W.HRρε i (Pi+1)

Figure 1 :

 1 Figure 1: Barcode of the M-ZZ.

Figure 2 :

 2 Figure 2: Barcode of the dM-ZZ.

Figure 3 :

 3 Figure 3: Barcode of the standard Rips filtration.For comparison we recall in Figure4the barcode obtained

Figure 4 :

 4 Figure 4: Barcode of the mesh-based filtration.

Note that we follow[START_REF] Carlsson | Zigzag persistence[END_REF] and depart from the traditional persistence barcode representation by using closed intervals instead of half-open intervals.

This commutative diagram induces the following interleaving between the oscillating Rips zigzag and the so-called image Čech zigzag, whose spaces are the images-denoted 2 See e.g.[START_REF] De | Coverage in sensor networks via persistent homology[END_REF] for a proof. Our definition of the Rips complex differs from the one in[START_REF] De | Coverage in sensor networks via persistent homology[END_REF] by a factor of 2 in the parameter value. This explains the slight discrepancy between our chain of inclusions and the one in[START_REF] De | Coverage in sensor networks via persistent homology[END_REF].

Acknowledgements

The authors wish to thank Marc Glisse for suggesting the approach used in the proof of Theorem 4.2. They also thank Frédéric Chazal for helpful discussions, and the anonymous referees for their insightful comments. This work was sup-ported by the European project No. 255827 (CG-Learning) and by the ANR project No. 09-BLAN-0331-01 (GIGA).