Mario Barbareschi
email: mario.barbareschi@unina.it

Alberto Bosio
email: alberto.bosio@ec-lyon.fr

Ian O'connor
email: ian.oconnor@ec-lyon.fr

Petr Fišer
email: fiserp@fit.cvut.cz

Marcello Traiola
email: marcello.traiola@inria.fr

A Design Space Exploration Framework for Memristor-Based Crossbar Architecture

Keywords: Emerging Technologies, Memristor crossbar, Design Space Exploration, Boolean Functions, Logic Synthesis

In the literature, there are few studies describing how to implement Boolean logic functions as a memristor-based crossbar architecture and some solutions have been actually proposed targeting back-end synthesis. However, there is a lack of methodologies and tools for the synthesis automation. The main goal of this paper is to perform a Design Space Exploration (DSE) in order to analyze and compare the impact of the most used optimization algorithms on a memristor-based crossbar architecture. The results carried out on 102 circuits lead us to identify the best optimization approach, in terms of area/energy/delay. The presented results can also be considered as a reference (benchmarking) for comparing future work.

I. INTRODUCTION

Today's computing devices are based on the CMOS technology, that is subject to the famous Moore's Law [START_REF]ITRS 2013 report[END_REF]. Despite the advantages of the technology shrinking, we are facing the physical limits of CMOS. Among the multiple challenges arising from technology nodes lower than 10 nm, we can highlight the high leakage current (i.e., high static power consumption), reduced performance gain, reduced reliability, complex manufacturing process leading to low yield and complex testing process, and extremely costly masks [START_REF] Hoefflinger | Chips 2020: A guide to the future of nanoelectronics[END_REF].

Several emerging technologies are under investigation, among them the memristor is a promising one [START_REF] Chua | Memristor-the missing circuit element[END_REF]. The memristor is a non-volatile device able to act as both storage and information processing unit that presents many advantages: CMOS process compatibility, lower cost, zero standby power, nanosecond switching speed, great scalability, high density, and non-volatile capability [START_REF] Waser | Redox-based resistive switching memories -nanoionic mechanisms, prospects, and challenges[END_REF], [START_REF]Memristive devices for computing[END_REF]. Thanks to its nature (i.e., computational as well as storage element), the memristor is used in different kinds of applications, such as neuromorphic systems [START_REF]Variation-tolerant computing with memristive reservoirs[END_REF], non-volatile memories [START_REF] Kim | A functional hybrid memristor crossbar-array/cmos system for data storage and neuromorphic applications[END_REF], or computing architecture for data-intensive applications [START_REF] Hamdioui | Memristor based computation-in-memory architecture for data-intensive applications[END_REF].

A fundamental component of any kind of computing architecture is the implementation of Boolean logic functions. In [START_REF] Xie | Fast boolean logic mapped on memristor crossbar[END_REF], the authors proposed a methodology for the synthesis of Boolean logic functions on a memristor crossbar. Their work showed that it is possible to implement any kind of Boolean function on a memristor crossbar. However, the experimental results have been carried out only on a couple of small circuits owing to the lack of a model with high abstraction level and also to the lack of an automated synthesis framework. In fact, they used a SPICE level model of the memristor crossbar that actually limits the possibility of dealing with realistic Boolean functions and of course the logic synthesis.

As defined in [START_REF] Testa | Logic synthesis for established and emerging computing[END_REF], the overall problem of logic synthesis is the one of finding "the best implementation" of a Boolean function. The term "best" corresponds to a trade-off between several metrics such as the area, delay, and power consumption. Today, with the rising of alternative technologies to CMOS, we are facing new challenges for logic synthesis. It is therefore mandatory to well understand the logic gates built on top of emerging technologies and identify the available opportunities. In [START_REF] Traiola | XbarGen: a memristor based boolean logic synthesis tool[END_REF], we presented an automatic tool for mapping a given Boolean function into a memristor crossbar. Thanks to this tool, it is possible to analyze two kinds of logic synthesis (i.e., 2-level and multi-level). In this paper, we extend [START_REF] Traiola | XbarGen: a memristor based boolean logic synthesis tool[END_REF] by considering additional approaches for logic synthesis and present a massive exploration phase. More in detail, the contributions can be summarized as follows:

• We analyze the impact of the different kinds of technology mapping to identify the best one for memristor based crossbars; • We quantify the circuits implementation w.r.t. area/energy/delay and we made a comparison with a CMOS based implementation; • The experimental results can be useful as a reference (benchmarking) for comparing future work. Results are published as open access repository [START_REF]Xbargen website[END_REF]. The remainder of the paper is structured as following. Section II presents the state-of-the-art and provides the required background about the memristor based computation. Section III details the proposed memristor and crossbar model as well as the synthesis framework, while Section IV gives the experimental results. Finally, Section V draws the conclusions.

II. BACKGROUND AND STATE-OF-THE-ART

In this section, we provide an overview of the state-of-theart about memristor-based circuits and the existing automatic synthesis approaches.

A. State-of-the-art overview

Memristors are currently under investigation and used in different types of computing architectures: from analog computation mainly exploited in neuromorphic computing [START_REF] Duan | Memristor-based cellular nonlinear/neural network: Design, analysis, and applications[END_REF], [START_REF] Wen | An EDA framework for large scale hybrid neuromorphic computing systems[END_REF], [START_REF] Shafiee | ISAAC: a convolutional neural network accelerator with in-situ analog arithmetic in crossbars[END_REF], up to digital architectures where the memristors are used to implement Boolean functions.

In [START_REF] Kvatinsky | Memristor-based material implication (IMPLY) logic: Design principles and methodologies[END_REF], the authors described the memristor-based IMPLY logic gate and proposed a methodology for designing it. Then, in [START_REF] Wang | Synthesis of memristive circuits based on stateful IM-PLY gates using an evolutionary algorithm with a correction function[END_REF], the authors presented an approach to synthesize memristor-based combinational logic circuits based on such logic. They proposed an evolutionary algorithm for obtaining a reduced number of working memristors. Also Raghuvanshi et al. considered implication-logicbased memristor circuits in [START_REF] Raghuvanshi | Logic synthesis and a generalized notation for memristor-realized material implication gates[END_REF]. Indeed, they illustrated new logic synthesis methods which aim at minimizing the number of implication gates, and thus the delay. Furthermore, in [START_REF] Kvatinsky | MAGIC -memristor-aided logic[END_REF], the authors addressed the memristor-based stateful logic problems (i.e., sequential voltage activations, result stored by one of the inputs, additional circuit components required) proposing a memristoronly logic family, i.e., memristor-aided logic (MAGIC). The authors of [START_REF] Chakraborti | BDD based synthesis of boolean functions using memristors[END_REF] took into account implication-logic-based memristor circuits. They represented Boolean functions as Reduced Ordered Binary Decision Diagrams (ROBDDs) and mapped them to 2-to-1 multiplexers implemented with memristors.

In [START_REF] Shirinzadeh | Fast logic synthesis for RRAM-based in-memory computing using majority-inverter graphs[END_REF], the authors presented an approach for the synthesis of memristive-based logic circuits using the Majority-Inverter Graphs (MIGs). They proposed optimizations of MIGs to reduce the number of memristive and computational steps in both MAJbased and IMP-based realizations. Xie et al. in [START_REF] Xie | Fast boolean logic mapped on memristor crossbar[END_REF] proposed a methodology that led to the synthesis of any kind of Boolean logic function on a memristor crossbar architecture. Afterwards, in [START_REF] Nguyen | Synthesizing HDL to memristor technology: A generic framework[END_REF], they proposed a generic synthesis framework to map logic circuits on memristor crossbars. In particular, it performs the synthesis of both combinational and sequential circuits. The framework takes HDL descriptions as input and generates circuits that can be simulated using electronic circuits simulators (i.e., SPICE). Furthermore, in [START_REF] Xie | A mapping methodology of boolean logic circuits on memristor crossbar[END_REF] they illustrated a methodology to map large Boolean logic circuits on a memristor crossbar. Specifically, they presented efficient place-and-route solutions and several optimization layouts for area, delay, and power consumption. They are able to apply their methodology to large logic circuits and to evaluate the performance of the design.

The scope of this work is in the latter category: digital crossbararray circuits. Solutions proposed in [START_REF] Xie | Fast boolean logic mapped on memristor crossbar[END_REF], [START_REF] Nguyen | Synthesizing HDL to memristor technology: A generic framework[END_REF], [START_REF] Xie | A mapping methodology of boolean logic circuits on memristor crossbar[END_REF] are very useful for the automation of the Physical Synthesis process that appears in [START_REF] Weste | CMOS VLSI Design: A Circuits and Systems Perspective[END_REF]). However, since there is no well-defined strategy for optimizing at the best a given Boolean function aiming to have the best memristor-based circuit according to design requirements, a massive and quick exploration phase is still needed. Thus, a higher level of abstraction is needed in order to speed up such process. We intend to investigate Logic Synthesis in order to have a comprehensive view of the optimization opportunities.

Next subsections describe the memristor model that we designed in our work. For this purpose, firstly we have to describe how the memristor is exploited in a crossbar architecture to perform a computation.

B. Memristor model

A memristor is a non-linear electrical component whose electrical resistance is not constant but depends on the history of the charge flowed through the device itself. Since we intend to implement a digital circuit, we refer to the memristor Voltage-Current relation depicted in Figure 1, detailed in [START_REF] Snider | Computing with hysteretic resistor crossbars[END_REF], as the best solution for modeling the memristor's behavior (i.e., thanks to the ideal response to a pulse-wave). Thus, as depicted in Figure 1, the voltage applied to the memristor's terminals does not change its resistance until it crosses a threshold. In the adopted ideal model, the upper and the lower thresholds have the same absolute value.

We resort to the Snider Boolean Logic (SBL) [START_REF] Snider | Computing with hysteretic resistor crossbars[END_REF] convention whereby a lower resistance (steeper curve denoted as R ON ) represents a logic '0' while a higher resistance (lower slope curve denoted as R OFF ) represents a logic '1'.

Two basic operations can be performed, defined as SET and RESET. The first one allows to program the memristor to R ON and thus at logic '0', while the second one programs the memristor to R OFF that corresponds to logic '1'. The Figure 1 depicts SET and RESET operations as described by Xie et al. in [START_REF] Xie | Fast boolean logic mapped on memristor crossbar[END_REF]. C. Fast Boolean Logic Circuits Snider proposed in [START_REF] Snider | Computing with hysteretic resistor crossbars[END_REF] a design methodology to implement Boolean functions on memristor-based crossbars. This design was then improved by Xie et al. in [START_REF] Xie | Fast boolean logic mapped on memristor crossbar[END_REF]. Let us briefly recall their proposition, by referring to it as Fast Boolean Logic Circuit (FBLC). Let us resort to an example by considering the following Boolean function:

O = AB + AB + AB = AB • AB • AB (1)
The left member of the Equation 1 can be easily manipulated through transformation rules (i.e., De Morgan's laws). The obtained form (right member of the Equation 1) can be computed by using three Boolean operations: NAND, AND and NOT.

Then, as sketched in Figure 2-a, the FBLC implementation of 1 is divided in blocks:

• Input block: it stores the inputs, in the example A, B and their complements; • Minterm block (M ): it configures and evaluates the minterms, in the example we have three minterms, one per row; • AND block: it stores the results of the minterms and performs the AND between them; • Output block: it stores the results of the AND block and performs the inversion operation to obtain the output value O in our example. To drive and control the different blocks a control circuitry is required. All the details can be found in [START_REF] Xie | Fast boolean logic mapped on memristor crossbar[END_REF].

III. TOWARDS AUTOMATIC CROSSBAR SYNTHESIS

This section presents the synthesis framework and the behavior model of the memristor-based crossbar.

A. XbarGen

XbarGen [START_REF] Traiola | XbarGen: a memristor based boolean logic synthesis tool[END_REF], [START_REF]Xbargen website[END_REF] analyzes the given input function to determine its levels. We refer to a "level" as a subset of the Boolean function having inputs and outputs such that:

• inputs are independent of each other; the inputs of a given level can be either the primary inputs or intermediate values;

• each output depends only on inputs (i.e, there are no intermediate nets); the outputs can be either primary outputs or intermediate values. Secondly, the framework performs a mapping to one or more crossbars, depending on the number of levels the input function is made of. For each level, inputs, outputs, and related minterms are mapped to a FBLC crossbar, as explained in [START_REF] Xie | Fast boolean logic mapped on memristor crossbar[END_REF] and briefly reminded in the previous section. When multiple crossbars are produced, they must be connected together in series, according to what was proposed by Snider in [START_REF] Snider | Computing with hysteretic resistor crossbars[END_REF]. Consequently, the latency (i.e., the computational time) of the circuit grows proportionally to the number of serially connected crossbars.

In this work, we study how different expressions of a Boolean function (i.e., logic synthesis) lead to different FBLC implementations. In turn, this will lead to to different trade-offs between performance (i.e., computational time) versus costs (i.e., area and energy consumption). The goal of the exploration is determining the best logic synthesis approach to obtain a Boolean function expression leading to the best FBLC implementation.

IV. EXPERIMENTS AND RESULTS

To the best of our knowledge, this is the first work comparing the effectiveness of existing logic minimization approaches for the synthesis of Boolean functions as memristor-based crossbars. In this section, we provide experimental results showing the impact of minimization approaches in terms of area/energy/delay.

We performed experiments following the flow described in Figure 4. Specifically, given a Boolean function, we employed different logic minimization (i.e., Logic Synthesis) approaches available in the state-of-the-art synthesis tool ABC [START_REF] Mishchenko | ABC: A system for sequential synthesis and verification[END_REF]. In particular, we have run these processes:

• Collapse: Collapsing into a two-level form by using the ABC collapse command • Collapse-E: Collapsing into a two-level form by using the ABC collapse command, followed by the two-level minimization by Espresso [START_REF] Brayton | Logic Minimization Algorithms for VLSI Synthesis[END_REF] • Collapse-Eso: Collapsing into a two-level form by using the ABC collapse command, followed by the two-level minimization by Espresso [START_REF] Brayton | Logic Minimization Algorithms for VLSI Synthesis[END_REF] with the -Dso option (singleoutput minimization) • LUTa: LUT mapping optimized for area by using the ABC command sequence dch; if -a -K K; mfs -W 10; st; dc2; if -a -K K; mfs -W 10 • LUTd: LUT mapping optimized for delay by using the ABC command sequence dch; if -K K; mfs -W 10; st; dc2; if -K K; mfs -W 10 Both LUT-mapping scripts were iterated 20-times, for different K, 2 ≤ K ≤ 31. Also, note that the synthesis with K = 2 performed mapping using arbitrary 2-input gates. We adopted such solution since a K-LUT can be easily mapped as a crossbar by using XbarGen [START_REF] Traiola | XbarGen: a memristor based boolean logic synthesis tool[END_REF], [START_REF]Xbargen website[END_REF]. The flow is depicted in Figure 3. Then, we employed XbarGen for obtaining statistics: occupied area, delay, and estimated energy consumption (switching activity) [START_REF] Traiola | Estimating dynamic power consumption for memristor-based CiM architecture[END_REF].

K-inputs Xbar

K-inputs Xbar

4-inputs Xbar

inputs outputs

4-inputs Xbar

4-inputs Xbar

3-inputs Xbar

inputs outputs

3-inputs Xbar

3-inputs Xbar

Mapping Mapping Mapping

Synthesis Synthesis

Fig. 3: Synthesis flow using K-input LUTS

We carried out experiments over 102 different combinational circuits from the LGSynth'91 [START_REF] Yang | Logic synthesis and optimization benchmarks user guide: Version 3.0[END_REF], IWLS'93 [START_REF] Mcelvain | IWLS'93 Benchmark Set: Version 4[END_REF], ISCAS'85 [START_REF] Brglez | A neutral netlist of 10 combinational benchmark circuits and a target translator in Fortran[END_REF], ISCAS'89 [START_REF]Combinational profiles of sequential benchmark circuits[END_REF], ITC'99 [START_REF] Corno | RT-level ITC'99 benchmarks and first ATPG results[END_REF], EPFL [START_REF] Amaru | The EPFL combinational benchmark suite[END_REF], and IWLS 2005 [START_REF] Albrecht | IWLS 2005 benchmarks[END_REF] benchmark sets, with up to 452 inputs, 256 outputs, and up to 78 thousands literals. It is worth to mention that some circuits were obtained by extracting the combinational logic of sequential circuits.

We first detail how we compute statistics. Table I reports the three metrics that we used: Area, Delay, and Energy. Area depends on the number of inputs, outputs, and minterms of each level in the Boolean function. These parameters depend on the applied logic synthesis. We do not report absolute values since we do not target any specific memristor technology. Therefore, the Area is measured as "units", where each unit corresponds to the size of a memristor. The Delay depends on the number of steps (N Steps ) required to evaluate the Boolean function implementation. In the case of a single crossbar implementation, it corresponds to 7 (see Section II). The absolute time is related to the memristor technology and it is expressed as t M : the time required to change the status of a single memristor. T C is thus the absolute time required to obtain the output values for a single crossbar. Finally, in Given from XbarGen, details available in [START_REF] Traiola | Estimating dynamic power consumption for memristor-based CiM architecture[END_REF] the case of several crossbars, the overall Delay time is computed as T C multiplied by N C . Again, the number of crossbars is provided by the logic synthesis. The last metric is the Energy. As for previous ones, it is related to the number of memristors used for implementing the given Boolean function. To quantify the energy consumption, we do not perform any simulations, we simply apply the methodology described in [START_REF] Traiola | Estimating dynamic power consumption for memristor-based CiM architecture[END_REF] allowing to estimate the average case (i.e., the average energy consumption).

A v e ra g e e n e rg y (s w it c h in g ) 0.0e+00 In Figure 5, we report the 3D representation of the obtained results w.r.t. the three metrics. To better discuss the results, we present three views obtained from Figure 5 by plotting two at a time. This leads to charts in Figures 6, 7, and8. For all charts, we depicted as blue dots the solutions LUTd (i.e., optimized w.r.t. delay) and as red dots the solutions LUTa (i.e., optimized w.r.t. area). Each dot is further characterized by its K value. In the figures, larger the K value, larger the plotted point dimensions.

The black dot always refers to the "Collapse" solution (i.e., the two-level synthesis). Last information about the charts is related to the fact that each dot (blue, red and black) corresponds to the average of the metrics extracted by synthesizing and mapping to the FBLC all 102 circuits (i.e., the overall number of syntheses is thus 3,366).

A. Area Vs. Delay

Figure 6 reports the trade-off between Area and Delay. First of all, it is clear from the results that the applied logic synthesis has a huge impact on the resulting circuit. Indeed, by using K-LUT mapping with K = 11 we obtained the best area. Increasing the K value makes the area larger while the delay drops. By collapsing the design to 2 levels, the best delay is obtained, but also the larger area. In this sense, we find the well-known areadelay trade-off characterizing the standard CMOS digital circuits. On the contrary, for K-LUT mappings with K < 11, both area and timing increase/decrease at the same time, thus not following the area-delay trade-off.

B. Energy Vs. Area

Figure 7 reports the relation between Energy and Area. First of all, the Collapse logic synthesis appears not to be the ideal solution. On the other hand, the LUT based synthesis with K = 11 provides one order of magnitude less energy consumption than the Collapse solution, with also smaller Area. For K > 11, both area and energy increase. This is expected, since bigger circuits usually consume more energy. On the contrary, for K < 11, energy consumption stays roughly constant around 10 3 while area changes in the range 10 5 -10 6 . depends on the Area (see Figure 6). Previous observations are thus valid in this case too.

C. Energy Vs. Delay

D. Comparison with CMOS

Finallty, Figure 9 wants to give an idea of the obtained results w.r.t. the CMOS counterpart, in terms of area and area-delay product. CMOS implementations were obtained by synthesizing the circuits with a standard 45nm 2-input gate library [START_REF]Freepdk45[END_REF]. To achieve a meaningful comparison, we needed technological parameter for the memristors: we used the following parameters, according to [START_REF] Borghetti | Memristive switches enable stateful logic operations via material implication[END_REF], 3nm × 3nm cell-size and 1 ns switching delay. As shown in the figure, the memristor implementations outperform the CMOS in terms of area. Furthermore, also in terms of areadelay product, we observe a similar result. Finally, since we did not model the system level, we do not compare memristor and CMOS implementations w.r.t. energy consumption and delay. Indeed, both should be measured while taking into account also data movement, which strictly depends on the system level. V. CONCLUSION The memristor is one of the most promising technologies which are able to deal with the CMOS limitations. In particular, since the memristor is inherently able to behave also as a non-volatile device, it allows overcoming the bottleneck of the data transfer to the computational unit and back to storage elements. The research community is facing the synthesis of Boolean functions by utilizing memristor-based crossbars. In order to speed up the exploration of new memristor-based crossbar architectures, our proposition is to move the analysis to a higher level (i.e., Behavioral).

The main goal was the analysis of the well-known logic optimization techniques such as the 2-level and the multi-level optimizations. Moreover, we also identified the best technology mapping process to further improve the optimization of logic synthesis for memristor-based crossbars. To achieve these results, we presented a Synthesis framework flow to explore the optimization process targeting a memristor-based architecture (i.e., FBLC). In the first phase, we applied Logic Synthesis to obtain a two-level (Collapse) or a LUT based implementation. Then, we map the results of Logic synthesis in memristor crossbar(s) and we quantified the final implementation by using Area, Delay and Energy as metrics. The above framework has been applied to a set of 102 benchmarks leading to 3,366 different implementations.

Results prove that K-LUT optimization approach provides the best results when considering Area/Energy/Delay as metrics. Finally, we believe that all the presented results, tool and HDL model can be considered as a valid benchmark set for comparing future work. For this purpose, it is possible to download all the presented material in [START_REF]Xbargen website[END_REF].

Fig. 1 :

 1 Fig.1: Set and Reset operations[START_REF] Xie | Fast boolean logic mapped on memristor crossbar[END_REF]

-

  Fig. 2: FBLC example

Fig. 5 :

 5 Fig. 5: Results: Area vs. Delay vs. Energy

Figure 8 Fig. 7 :

 87 Figure8reports the relation between Energy and Delay. This chart is qualitatively similar to Figure7(i.e., for K < 11, energy consumption stays roughly constant and for K > 11, it increases while the delay decreases). Again, this was expected, since Delay

Fig. 8 :

 8 Fig. 8: Results: Energy vs. Delay

Fig. 9 :

 9 Fig. 9: Comparison w.r.t. CMOS in terms of average area and average area-delay product

  Energy Computed as the worst case of memristors switching activity

		Boolean	Collapse		Circuit VHDL projects	System
		Function				Xbar(s)
			LUTa	Mapping
		Logic Synthesis	LUTd	XbarGen	Statistics	Analysis
			Fig. 4: Synthesis Framework Flow
			TABLE I: Legend
	Symbol		Description			Expressed as
	p ij	Whether the i-th minterm is in the j-th level (Boolean)		Given by Logic Synthesis
	N in (l j )	Num. of input for the j-th level			Given by Logic Synthesis
	N out (l j )	Num. of output for the j-th level			Given by Logic Synthesis
	Area		Total area
	t M	Memristor switching time		Given from technological characterization
	N Steps	Num. of steps of a Crossbar computation		Given by the architecture (equal to 7)
	T C	"Latency" of a single Crossbar			t M * N Steps
	N C	Num. of Crossbar in the circuit			Given by Logic Synthesis
	Delay	Needed Cycles to complete while computation		T C * N C

j [2 * N in (l j ) + 2 * N out (l j )] * 1 + i (p ij ) + N out (l j )


Annachiara Ruospo

Gabriele Gavarini

Ilaria Bragaglia

Marcello Traiola

Alberto Bosio

Ernesto Sanchez

Selective Hardening of Critical Neurons in Deep Neural Networks

Keywords: Deep Neural Network, Reliability, Mitigation, Fault Injection

In the literature, it is argued that Deep Neural Networks (DNNs) possess a certain degree of robustness mainly for two reasons: their distributed and parallel architecture, and their redundancy introduced due to over provisioning. Indeed, they are made, as a matter of fact, of more neurons with respect to the minimal number required to perform the computations. It means that they could withstand errors in a bounded number of neurons and continue to function properly. However, it is also known that different neurons in DNNs have divergent fault tolerance capabilities. Neurons that contribute the least to the final prediction accuracy are less sensitive to errors. Conversely, the neurons that contribute most are considered critical because errors within them could seriously compromise the correct functionality of the DNN. This paper presents a software methodology based on a Triple Modular Redundancy technique which aims at improving the overall reliability of the DNN, by selectively protecting a reduced set of critical neurons. Our findings indicate that the robustness of the DNNs can be enhanced, clearly, at the cost of a larger memory footprint and a small increase in the total execution time. The trade-offs as well as the improvements are discussed in the work by exploiting two DNN architectures: ResNet and DenseNet trained and tested on CIFAR-10.

I. INTRODUCTION

In the last decades, the growing complexity of emerging computing systems has called for enhanced computing paradigms. Among all the existing possibilities, artificial intelligence (AI) based solutions and, specifically, brain-inspired computing models like Artificial Neural Networks (ANNs), and in particular Deep Neural Networks (DNNs), have gained large interest in the industry and academia for their near-human computational capabilities [START_REF] He | Delving deep into rectifiers: Surpassing human-level performance on imagenet classification[END_REF]. DNNs are now considered attractive solutions for different areas including safety-critical applications such as selfdriving cars, radars, flight control, robots, and space applications. Therefore, there is a compelling need for ensuring their reliability and tolerance to faults.

Being brain-inspired models, DNNs are traditionally considered intrinsically fault-tolerant and tightly robust. The first reason is related to their distributed and parallel structure, the second one to the redundancy introduced because of the over-provisioning [START_REF] Lawrence | What size neural network gives optimal generalization? convergence properties of backpropagation[END_REF]. Indeed, neural networks are made of a number of neurons higher than the minimal required to perform a computation. It means that they can tolerate some errors due to the excessive neurons budget [START_REF] Haykin | Neural networks and learning machines[END_REF]. This property implies that DNNs can tolerate a bounded number of potentially failing neurons without impacting the result of the computation; beyond that number, the precision degrades gracefully [START_REF] El Mhamdi | When neurons fail[END_REF]. Nevertheless, when DNNs are deployed on resource-constrained hardware devices, single physical faults on hardware units might jeopardize the activity of multiple neurons, depending on how the DNN workload is scheduled on hardware resources [START_REF] Ruospo | On the reliability assessment of artificial neural networks running on ai-oriented mpsocs[END_REF]. The corruption of multiple neurons may lead to misprediction or, generally, unexpected outcomes. It is worth adding that not all neurons play the same role in a neural network: in the literature, it is claimed that individual network parts differ in their error resilience [START_REF] Venkataramani | Axnn: Energy-efficient neuromorphic systems using approximate computing[END_REF], [START_REF] He | Joint design of training and hardware towards efficient and accuracy-scalable neural network inference[END_REF]. Neurons within a neural network turn up having different fault tolerance and resilience levels. Some of them strongly contribute to the final task and their failures have a large influence on the degradation of the final predictions; others turn up being less important to the output accuracy (the concept is connected to the property of the over provisioning). Thanks to this, the latter are often removed from the network, for example through the pruning technique [START_REF] Yu | Nisp: Pruning networks using neuron importance score propagation[END_REF]. In the literature, the problem of assessing the resilience of artificial neurons and especially identifying the critical ones is addressed in many research projects. The proposals are different and are oriented towards two main directions: studying the contribution of each neuron to the global error during the training phase ( [START_REF] Venkataramani | Axnn: Energy-efficient neuromorphic systems using approximate computing[END_REF], [START_REF] Zhang | Approxann: An approximate computing framework for artificial neural network[END_REF], [START_REF] Wang | Enabling energy-efficient and reliable neural network via neuron-level voltage scaling[END_REF]), or their contribution to the final network prediction during the inference phase ( [START_REF] Ruospo | On the reliability assessment of artificial neural networks running on ai-oriented mpsocs[END_REF], [START_REF] Schorn | Accurate neuron resilience prediction for a flexible reliability management in neural network accelerators[END_REF]).

In this work, we present an approach to selectively protect a subset of neurons in deep neural networks by using a Triple Modular Redundancy technique. Firstly, we use one of the most recent methodologies to classify neurons based on their criticality (i.e., [START_REF] Ruospo | On the reliability assessment of artificial neural networks running on ai-oriented mpsocs[END_REF]). Then, to realize the selective hardening, we triple the computations and the parameters of a set of critical neurons. Hence, under the single fault assumption, faults are corrected thanks to the adoption of a majority voter. Clearly, this implies advantages and drawbacks: the reliability of the DNN is enhanced at the cost of increased memory footprint and Multiply-ACcumulate (MAC) operations. The aim of this work is to investigate (i) the existing trade-offs between memory footprint, execution time, and reliability, and (ii) the best percentage (pNeu%) of protected neurons required to achieve a substantial increase in the overall reliability without heavily impacting the neural network size and execution time.

The rest of the paper is organized as follows. Section II provides the reader with background knowledge on critical neurons and presents the related studies. Section III describes the proposed approach and Section IV outlines the case study. Next, Section V reports on the experimental results. Finally, Section VI draws conclusions and future directions.

II. BACKGROUND

The intent of this section is to provide a basic background knowledge in the field. First, an overview of the existing fault models in artificial neurons is provided. Then, some of the most relevant related works in the literature are discussed.

A. Fault Models in Artificial Neurons

In the neural network field, each artificial neuron is considered as a single entity which can fail independently of the failure of any other [START_REF] El Mhamdi | When neurons fail[END_REF]. Neural networks are viewed as distributed systems consisting of two principal components: neurons and synapses, i.e., the communication channels connecting the neurons. Faults in neurons or synapse are therefore independent of those of others. Specifically, it is worth to point out that, as neuron, we refer to each pixel in the output feature maps of convolutional layers and to each node in the pooling (min, max, average) or fully connected layers. Normally, batch normalization and activation functions (e.g., rectified linear unit, sigmoid, Gaussian) are not considered as independent layers, and thus they do not provide the network with additional neurons. In the literature, two fault models cover the most common faults in neurons:

• Crash: Neurons completely stop their activity. A crashed neuron is modelled by purposely setting its output to zero. This can also be known as dropout fault model [START_REF] Schorn | Accurate neuron resilience prediction for a flexible reliability management in neural network accelerators[END_REF]. • Byzantine: Neurons keep their activity but produce arbitrary values, within their bounded transmission capacity [START_REF] El Mhamdi | When neurons fail[END_REF], [START_REF] Allen | An evaluation of causes for unreliability of synaptic transmission[END_REF]. It is worth underlining that a fault affecting a single neuron may not lead to a failure. This is not only due to the intrinsic definition of a fault [START_REF] Bushnell | Essentials of Electronic Testing for Digital, Memory and Mixed-Signal VLSI Circuits[END_REF], but specifically for the ANN property of being overprovisioned.

B. Related Works

Over the past decades, the understanding of the reliability, fault tolerance, and robustness of neural networks has gained increasing attention around the world [START_REF] Torres-Huitzil | Fault and error tolerance in neural networks: A review[END_REF]. Among all the possible aspects, understanding the importance of individual neurons takes on great relevance when running deep and complex DNN models on systems with limited computing and memory resources. To mitigate the problems related to the size of DNNs, in the literature, many researchers provided pruning techniques to remove either redundant neurons or connections from over-parameterized neural models. However, we advocate that, to properly remove neurons or weights, a proper resiliency analysis is needed. In fact, as mentioned before, neural network parts differ in fault tolerance capability.

Spyrou et al. in [START_REF] Spyrou | Neuron fault tolerance in spiking neural networks[END_REF] propose a neuron fault tolerance strategy for Spiking Neural Networks (SNNs). They propose a fault tolerance mechanisms based on dropout to nullify the effect of certain faults (by temporarily removing neurons during training with some probability) and some active fault tolerance techniques to detect and recover from the remaining faults. In [START_REF] Venkataramani | Axnn: Energy-efficient neuromorphic systems using approximate computing[END_REF], [START_REF] Zhang | Approxann: An approximate computing framework for artificial neural network[END_REF], the problem of identifying critical neurons in neural networks is addressed by using backpropagation of error gradients to find out those neurons that impact the output quality the least. The overall idea is to approximate those neurons that contribute the least to the global error with energy-efficient ones. The algorithms adopt the training set to forward and backward propagate values and then assign resilience scores based on the magnitude of their average error contribution. To the same end, Wang et al. in [START_REF] Wang | Enabling energy-efficient and reliable neural network via neuron-level voltage scaling[END_REF] address the problem of the voltage scaling technique applied to neural networks. This method is widely used for energy saving, but it may cause significant reliability challenges. They propose to build voltage islands with the same voltage to neurons with similar fault tolerance capability. To identify critical neurons, they use the derivative of the cost function to represent the final output quality degradation caused by the error of the neuron's calculation.

Afterwards, authors in [START_REF] Schorn | Accurate neuron resilience prediction for a flexible reliability management in neural network accelerators[END_REF] propose a method that differs from the above-mentioned gradient-based approaches ( [START_REF] Venkataramani | Axnn: Energy-efficient neuromorphic systems using approximate computing[END_REF], [START_REF] Zhang | Approxann: An approximate computing framework for artificial neural network[END_REF], [START_REF] Wang | Enabling energy-efficient and reliable neural network via neuron-level voltage scaling[END_REF]). They characterize the resilience of individual neurons based on their average contribution to the output of the DNN. In other words, they state that neurons which strongly contribute to the output of a neural network have a greater impact on the classification accuracy degradation in case their output is disturbed. Their method bases on a Taylor decomposition and layer-wise relevance propagation (LRP) algorithm to compute the value for each neuron and shows to be more effective than the gradient-based ones.

Based on the same observation (the average contribution to the output of a neural network), authors in [START_REF] Ruospo | On the reliability assessment of artificial neural networks running on ai-oriented mpsocs[END_REF] recently proposed a similar approach to assign resilience scores to neurons while adding a key contribution. Unlike existing works in the literature which can be considered as network-oriented, they claim that in a multioutput neural network, the contribution of a single neuron can have a significant importance for individual classes and not only for the entire neural networks. As a matter of fact, neuroscience theories state that neurons become active in a particular pattern of neuronal activity based on precise input stimuli [START_REF] Squire | Memory systems of the brain: A brief history and current perspective[END_REF]. Therefore, they propose a class-oriented analysis and an algorithm to sort neurons based on both their contribution to individual classes and to the entire neural network.

III. PROPOSED APPROACH

This work presents a methodology to improve the reliability of DNNs by selectively hardening subsets of neurons based on their criticality. It leverages a Triple Modular Redundancy (TMR) mechanism to protect the functionality of the most important computations and prevent possible errors in them from propagating to the DNN output.

The approach builds on a preliminary phase: the ranking of neurons according to their criticality and resilience values. Among all the mentioned techniques (Section II-B), our interest landed on the methodology proposed in [START_REF] Ruospo | On the reliability assessment of artificial neural networks running on ai-oriented mpsocs[END_REF], where a traditional networkoriented analysis (neurons only considered as entities of the entire neural network) is combined with a class-oriented one (every neuron is likewise regarded as an entity of each individual class). In the following, we briefly summarize the approach. Specifically, neurons are sorted based on the magnitude of their average contribution over the training set, as also suggested in [START_REF] Schorn | Accurate neuron resilience prediction for a flexible reliability management in neural network accelerators[END_REF]. The algorithm proposed in [START_REF] Ruospo | On the reliability assessment of artificial neural networks running on ai-oriented mpsocs[END_REF] follows this procedure: (i) class-oriented analysis and (ii) network-oriented analysis. In the former, for each instance in the validation dataset related to the specific output class, a forward propagation cycle is performed, and a score is given to each neuron by averaging the absolute output values In this paper, the proposed mitigation approach starts from the neurons criticality list, ordered based on the neurons criticality, and obtained as specified above. A diagram of the proposed flow is shown in Fig. 1. First, as a measure of the desired level of reliability, a specific percentage of neurons (pNeu) must be selected to perform the following actions (detailed later):

1) Prune non-critical neurons, i.e., eliminate the selected pNeu from the bottom of the criticality list of neurons; 2) Reinforce critical neurons through a TMR-based approach where the pNeu is used to select a subset of critical neurons (TMR-pNeu) from the top of the list.

Although selectively hardening critical neurons allows increasing the reliability level of the DNN model, it comes at the cost of increasing the memory footprint (i.e., increasing weights and biases connected to TMR-pNeu) and also the number of MAC operations, which in turn increases the execution time. At the same time, the introduced computational overhead (in terms of MAC operations) is mitigated by pruning non-critical neurons.

The reader should notice that, as the objective of this work is to improve the fault tolerance of a DNN model without increasing, as far as possible, its computational costs, the proposed solution leverages the removal of the same quantity of protected neurons. This approach is based on the neuron output sparsity: a careful analysis of the list of critical neurons revealed the occurrence of neurons propagating, on average, zero values, that were located at the bottom of the list (non-critical neurons). Thus, the idea is to prune them to reduce the total number of MAC operations without affecting (or while slightly affecting) the DNN model accuracy. A similar approach is also proposed in [START_REF] Reagen | Minerva: Enabling low-power, highly-accurate deep neural network accelerators[END_REF], where the authors use a threshold to remove neuron activities below the threshold and therefore prune non-significant neurons. During this pruning phase, as shown in the diagram (Fig 1 ), not only the number of MAC operations is reduced, but also the memory footprint: indeed, Fig. 2: The proposed methodology applied to a 3-layer feedforward neural network: an illustrative example.

if neurons belonging to fully connected layers are removed, their weights and biases are deleted as well.

Afterwards, given the same initial percentage (pNeu), the selective hardening of critical neurons is applied. In more details, each neuron in the set (TMR-pNeu) is replaced with a TMRbased structure: not only its calculation is tripled, but also all the parameters related to it. Then, before transferring the value to subsequent neurons, a majority voter looks for the consistency of the three results. Basically, it ensures correctness, provided at least two out of three copies of the neuron remain operational and faultfree. Otherwise, if the values of the three copies do not match, the majority voter turns off the whole computation of the critical neuron, to avoid the propagation of wrong or, even worse, out-ofrange values. As a result, the number of MAC operations and the memory footprint of the DNN application unavoidably might rise, but the reliability level of the DNN is improved.

For the sake of clarity, let us consider the simple 3-layer feedforward neural network illustrated in Fig. 2. In the example, pNeu is equal to 20%: a TMR-20 is applied, i.e., 20% of the critical neurons (at the top of the list) are protected and 20% of non-critical neurons (at the bottom of the list) is removed. The reader can note that: first, all the parameters (weights and biases) are tripled (the blue synapse), and second, the connections between the original critical neurons and the subsequent layers are cut and are replaced with those of the majority voter. Similarly, weights associated with non-critical neurons (if they are not shared -such as in fully connected layers), are removed as well.

To conclude, at the end of the proposed flow, if the designers are not satisfied with the resulting parameters (reliability level, MAC operations, and memory footprint), they can adjust the pNeu percentage to meet their goals.

IV. CASE STUDY

To prove the effectiveness of the proposed methodology, we used four different DNNs: ResNet-20, ResNet-32 [START_REF] He | Deep residual learning for image recognition[END_REF], DenseNet-121, and DenseNet-161 [START_REF] Huang | Densely connected convolutional networks[END_REF]. The number of parameters and neurons of such DNN models are shown in Table I. All the DNNs were trained using CIFAR-10 [20], a well-known dataset used for the task of image classification. This dataset is composed of 60k 32x32 coloured instances, 50k training images and 10k test images. The DNN models are trained and tested with a 32-bit single-precision floating-point representation by using the PyTorch framework on a Linux server equipped with a dual Intel Xeon CPU E5-2680 v3 and 256 GB of RAM. To demonstrate the validity of the methodology, fault injection (FI) campaigns are performed. We resorted to the crashed neurons fault model: we model faulty neurons as neurons that completely stop their activity and the propagation of values. In the literature, this is also known as the dropout fault model. In practice, this error is introduced in the DNN by forcing to zero the output of the targeted faulty neuron. In the last experiment, we resorted to the Byzantine neurons fault model to model neurons that propagate random values. To run such FI campaigns, a specific library was built based on PyTorch to allow the injection of single and multiple faults on neurons. In more details, the fault injection is implemented by applying a binary mask to the output feature map of each layer at inference time. Each entry of the mask represents the occurrence of faults on the corresponding element of the target feature map. The computations of the masks are done offline, before running the injection campaigns.

V. EXPERIMENTAL RESULTS

To demonstrate the validity of the adopted ranking of neurons, we performed a first set of cumulative FI experiments where the same growing quantity of neurons was crashed (their outputs are cumulatively set to zero) in two scenarios. The reader should note that in this case the multiple fault assumption better reflects reality: in modern architectures, due to the size and the complexity of state-of-the-art neural networks, the computations of multiple neurons are assigned to a single hardware processing element. It means that a single physical fault affecting a processing element might compromise the activity of multiple neurons.

In the first scenario (Critical), the faulty neurons were selected from the top of the list of the critical ones; in the second (Random) they were randomly chosen from among all neurons. In both cases, the DNNs accuracy is reported in Fig. 3a and 3b for all the examined neural networks. The x-axis represents the increasing percentage of cumulative crashed neurons, whereas the y-axis shows the average accuracy that the DNN under assessment achieves by running all the instances of the CIFAR-10 Testset. As for the Random scenario, the FI experiments are repeated 10 times for each point of the x-axis (i.e., 1, 5, 10, 25 %): every time random neurons are selected (including the possibility of having neurons in the set of critical). As it turns out, the accuracy of the DNNs falls dramatically if faults occur in critical neurons: this is where our mitigation proposal comes from. As stated, it leverages redundancy to (i) prune non-critical neurons and (ii) protect critical ones.

Therefore, to show that pruning does not affect, up to a certain point, the accuracy of the neural networks, we performed a second set of experiments. Fig. 4 shows how DNNs accuracy varies if a certain percentage of non-critical neurons is eliminated from the computations. It emerges that, in line with the redundancy property of DNNs, the deeper the DNN model, the higher is the quantity of non-critical neurons that you can remove without affecting the prediction accuracy of the network. For example, in DenseNet-161 it is possible to remove up to 15% of non-critical neurons while keeping the same accuracy. On the other hand, in reduced-size neural networks, such as ResNet-20, this bound is reduced to 5%. The pruning of non-critical neurons allows to exploit the intrinsic redundancy property of neural networks and, at the same time, avoid an excessive increase in the number of computations: the proposed selective TMR repeats three times all the arithmetic operations of critical neurons. Based on the outcome of this experiment (Fig. 4), in the selective hardening phase pNeu was limited to 5 percent to ensure that, in all our case studies, pruning would not affect the accuracy and, thus, the behaviour of the DNNs.

In this light, in Table II we propose a complete analysis of costs: for every DNN and TMR-pNeu, both the variation of the memory footprint of the DNN application and the execution time (in terms of number of MAC operations) is reported. Specifically, the selective hardening methodology was applied with three levels of protection: TMR-1, TMR-2.5, and TMR-5. Data in Table II demonstrate that, thanks to the pruning of redundant neurons, the proposed methodology leads, on average, to a minimal increase in execution time (the maximum is reached with ResNet-20 and corresponds to +0.12%). It is interesting to point out that, in some scenarios, despite the unavoidable redundancy introduced by the TMR, the total number of MAC operations reduces. To explain this phenomenon, we resort to the formula used to compute the number of MACs per neuron (Eq. 1).

M AC neu = C in * K h * K w (1)
where C in is the number of channels in the input feature map (its depth) while K h and K w the height and the width of the kernel of weights, respectively. This number is intuitively proportional to the size of the input channel and the kernel. In our DNN models, all the kernels have the same size (3x3). It means that Eq. 1 becomes only proportional to the first parameter (C in). We have a reduction of MAC operations (Table II) if the TMR-pNeu is applied to neurons belonging to layers with a small C in. II), in a further set of experiments we show the improvements in terms of reliability and robustness. Fig. 5 illustrates the outcomes of FI experiments, reproducing the worst-case scenarios: when critical neurons are crashed and completely stop their transmission activity. In the graphs, the selective hardening is used with three values of pNeu: 1, 2.5 and 5%. Aiming at reproducing the worst-case scenarios, the injections are not performed randomly, but following the order of critical neurons. Clearly, the accuracy of the DNNs remains stable when the same percentage of protected neurons crash. After that point, the accuracy of the DNNs slightly and smoothly decrease, specifically when a TMR-5 is applied. It means that, on the critical list, subsequent unprotected neurons are not as critical as previous ones. DenseNet-161, for instance, thanks to the TMR-5 maintains a prediction accuracy greater than 80% up to 8% of crashed neurons (the top 5% is protected by the selective hardening), in the worst case. Clearly, both in TMR-1, TMR-2.5, and TMR-5, the pruning of non-critical neurons is applied together with the hardening.

Cumulative Byzantine Neurons: As discussed in Section II, faults in neurons can be also modelled by forcing random values at their output (Byzantine fault model). We performed a further experiment to illustrate the behaviour of the DNNs when (i) random faults affect the unprotected DNN and (ii) random faults affect the protected DNN with different levels of TMR. The entire range of values that can be represented on 32 bits was used to force Byzantine faults. Experimental results are shown in Fig. 6 for ResNet-32. The x-axis reports the number of cumulative Byzantine faults that are introduced in random neurons (the experiments are repeated 10 times); the y-axis depicts the average accuracy that ResNet-32 achieves in four cases: with its original shape (i.e., unprotected), and with three different TMRs. For each point of the x-axis, the entire CIFAR-10 testset is run. Even with a reduced set of cumulative faults, it is clear that the DNNs with protected neurons have an enhanced ability to tolerate faults.

VI. CONCLUSIONS

This paper describes a technique to improve the reliability of DNNs by protecting the correct functionality of critical neurons. It is done by exploiting a Triple Modular Redundancy (TMR) mechanism combined with the pruning of non-critical neurons. Indeed, this work builds on two important and well-known factors: the over-provisioning of DNNs (i.e., their redundancy), and the different fault tolerance capability of neurons inside DNNs. The proposed approach comes with advantages and disadvantages, but, overall, it shows that the reliability and the robustness of the neural network can be improved if the computation of a reduced percentage of critical neurons is protected by a software TMR.

Fig. 1 :

 1 Fig. 1: Proposed Flow

Fig. 3 :

 3 Fig. 3: Critical Neurons Vs. Random Neurons.

Fig. 4 :

 4 Fig. 4: Variation in DNNs accuracy due to the pruning of noncritical neurons.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Proposed Methodology with Crashed Neurons

TABLE I :

 I DNNs Description

		Total Parameters [k] Total Neurons [k]
	ResNet-20	270	188
	ResNet-32	464	303
	DenseNet-121	6,956	578
	DenseNet-161	26,483	926

  As an example, most of the MAC operations added by the TMR in DenseNet-161 are due to critical neurons in the first layer: up to TMR-5, only in the first one where C in is equal to 3. In deeper layers of DenseNet-161, C in reaches a depth of 2,064 channels. In contrast, neurons removed by pruning are withdrawn from various layers (of different size, typically from deeper layers with deeper channels). This justifies why, despite the redundancy introduced by the proposed technique, the execution time might decrease. The same applies to the other DNNs under study. However, critical neurons in ResNet-20 and ResNet-32 are principally concentrated in layers with a small C in, but also spread in deeper ones. Similarly to DenseNet-121 and DenseNet-161, the pruning acts in different layers. It is worth to point out that the maximum C in in ResNet-20 and ResNet-32 is at most 64 (against a 2,064 of DenseNet-161). As regards the memory footprint, Table II evidences a general increase for every DNN and every TMR-pNeu. As explained in Section III, to protect the computation of a specific pNeu of critical neurons, all their weights and, overall, parameters, are replicated. Typically, kernels' size in initial layers is lower compared to that of the last ones and the most critical neurons are mostly located in them.Cumulative Crashed Neurons: After discussing costs and gains of the proposed technique (Table

TABLE II :

 II Overhead of the proposed technique.

	TMR-%	ResNet-20 Memory Footprint [kB] MACs [M]	ResNet-32 Memory Footprint [kB] MACs [M]	DenseNet-121 Memory Footprint [kB] MACs [M]	DenseNet-161 Memory MACs Footprint [kB] [M]
	TOTAL*	1,070	40.55	1,844	68.86	27,450	223.51	104,962	619.11
	TMR-1	+7.67%	+0.04%	+1.66%	-0.01%	+0.03%	-0.16%	+0.01%	-0.20%
	TMR-2.5	+10.06%	+0.05%	+7.13%	+0.03%	+0.03%	-0.37%	+0.01%	-0.42%
	TMR-5	+22.56%	+0.12%	+15.45%	-0.08%	+0.04%	-0.87%	+0.01%	-0.76%

* Without the proposed techniques (pruning+hardening).


Gabriel Bathie

Gaétan Berthe

Yoann Coudert-Osmont

David Desobry

Amadeus Reinald

Mathis Rocton

PACE solver description: DreyFVS

Keywords: 2012 ACM Subject Classification Theory of computation → Graph algorithms analysis; Theory Directed Feedback Vertex Set, Heuristic, Sinkhorn algorithm, Local search

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In the following, a directed graph, or digraph G = (V, E) is a set of vertices V , and a set of arcs E, consisting of ordered pairs of vertices. For a vertex v ∈ V (G), we denote its out-neighbourhood by N + (v) = {u ∈ V (G) : (v, u) ∈ E}, and its in-neighbourhood by N -(v) = {u ∈ V (G) : (u, v) ∈ E}. A directed cycle is a sequence (u 1 , ..., u k , u 1 ) such that (u i , u i+1 [k] ) ∈ E(G). A directed acyclic graph, DAG for short, is a digraph containing no cycles. Given X a set of vertices, we let G[X] be the graph induced by X. A strongly connected component of G is a vertex set X such that there is a path between any pair of vertices in G[X]. We say that (u, v) is a digon if it is a (directed) cycle of length two. We say that a set C ⊆ V is a d-clique whenever |C| = d and each pair (c, c ′ ) ∈ C forms a digon. The Directed Feedback Vertex Set problem takes as input G = (V, E) and asks for a set X ⊆ V such that G[V \X] contains no directed cycle. We consider the minimization version of the problem, asking for X to be of minimal size.

Our algorithm first applies reduction rules to obtain a smaller instance, guesses a first solution using a Sinkhorn-Knopp algorithm, then uses two heuristics to improve the solution in a local manner. A C++ implementation is available on a public git repository [1].

In Section 2, we present the steps needed to guess our first solution. We first describe our reduction rules 2.1, based on those described by Lin and Jou [START_REF] Lin | On computing the minimum feedback vertex set of a directed graph by contraction operations[END_REF], and provide an extension to their CORE rule. Then, we give an overview of the Sinkhorn-Knopp method in 2.2. In Section 3, we describe our two local search heuristics.

2

Guessing a solution

Reductions

At any step of the reduction, one can deal with each strongly connected component independently, as solutions to the instance are exactly unions of solutions to each strongly connected component. The main reductions rules applied in our algorithm are those developed in [START_REF] Lin | On computing the minimum feedback vertex set of a directed graph by contraction operations[END_REF], namely LOOP, IN0, IN1, OUT0, OUT1, PIE, CORE and DOME. Given a vertex v ∈ V (G), a bypass of v consists in the deletion of v, followed by the addition of arcs (u, w) for any

u ∈ N -(v), w ∈ N + (v).
The goal of the bypass operation is to reduce instances where v can always be replaced by another vertex in the solution. Towards motivating our extension of the CORE rule, we first describe rule IN1, note that the OUT1 rule is symmetric when v admits a single out-neighbour.

▶ Reduction Rule 1 (IN1). If a vertex v ∈ V (G) admits a single in-neighbour, bypass v.
Indeed, the in-neighbour of v is part of any cycle containing v. A solution using v in the original graph yields a solution of the same size using its in-neighbour in the reduced graph. Thus, a solution in the reduced graph is also a solution in the original one.

In [START_REF] Lin | On computing the minimum feedback vertex set of a directed graph by contraction operations[END_REF], the CORE rule looks to bypass a vertex v that is a CORE of some d-clique, meaning that it forms a d-clique with all its (in or out) neighbours. We extend this rule by applying it whenever its in-neighbours or its out-neighbours form a d-clique.

▶ Reduction Rule 2 (CORE'). If a vertex v ∈ V (G) is such that N -(v) or N + (v) forms a d-clique, bypass v.
Indeed, considering the case where N -(v) forms a d-clique, any solution to DFVS must contain at least one vertex per digon in N -(v), thus at least d -1 vertices of N -(v) in total. Then, by the same arguments as rule IN1, taking the remaining in-neighbour of v in a solution is at least as good as taking v. The case is symmetrical when N + (v) forms a d-clique, using the OUT1 rule.

The Sinkhorn-Knopp Algorithm

Given a graph G, we construct our first solution X by the heuristic of Shook and Beichl [START_REF] Shook | Matrix scaling: A new heuristic for the feedback vertex set problem[END_REF]. We call a set of vertex disjoint cycles a disjoint cycle union, or DCU for short. Given an n × n matrix A, the permanent of A is defined as perm(A) = σ∈Sn n i=1 a iσ(i) . The idea of the heuristic is to iteratively add vertices v that are contained in many DCUs to the solution X. While computing all DCUs is not efficient, a first observation is that the permanent of A allows us to count the number of spanning DCUs in G, where a spanning DCU is a DCU covering all vertices of the graph. To lift the spanning constraint and count all DCUs, it is possible to take the permanent of the matrix M = A + I. If we denote by M ij the matrix M with row i and column j removed, we define the m-balance matrix m-bal(M ) as in [START_REF] Beichl | Approximating the permanent via importance sampling with application to the dimer covering problem[END_REF]:

m-bal(M ) ij = m ij • perm(M ij ) perm(M ) ,
Then, m-bal(M ) vv is the fraction of DCUs that do not contain v. Thus, small values lead us to consider the addition of v to the solution. Computing m-bal(M ) is still hard, nevertheless [START_REF] Beichl | Approximating the permanent via importance sampling with application to the dimer covering problem[END_REF] provides a way to approximate m-bal(M ). This is done by computing the Sinkhorn balance matrix s-bal(M ) through the Sinkhorn-Knopp algorithm [START_REF] Sinkhorn | Concerning nonnegative matrices and doubly stochastic matrices[END_REF] applied on M . Sinkhorn-Knopp is an iterative algorithm where each iteration has a O(|E|) complexity. We decided to apply log |V | iterations leading to a O(|E| log |V |) complexity for the computation of s-bal(M ). We can now describe our iterative construction of solution X. As long as graph G is not empty, we compute s-bal(M ) and add the vertex v minimising s-bal(M ) vv to X, then, we apply our reduction rules to the remaining graph G \ {v} and re-iterate.

Local search

The solution obtained in the previous section is further improved using two greedy local search heuristics sequentially. In the following, G refers to a strongly connected component of the reduced instance, and X to our current solution for DFVS.

Vertex Swapping

The first local operation consists of a simple vertex swap. Namely, we remove a vertex from the current solution X, and if a cycle is created we attempt to replace it with another vertex of the cycle at random. We iterate through every vertex v of X. If G \ X ∪ {v} is acyclic, we remove v from X such that the solution size decreases by 1. Otherwise, we consider a cycle C in G \ X ∪ {v}, and choose some other vertex w ∈ C \ {v}, replacing v with w in X if G \ X ∪ {v} \ {w} is acyclic.

Vertex Ordering Perturbation

Recall that a digraph is acyclic if and only if its set of vertices admits a topological ordering, that is, an order π over V such that for every arc (u, v), π(u) < π(v). Our idea, similar to that of [START_REF] Galinier | Applying local search to the feedback vertex set problem[END_REF], is to define a unique DFVS X π from any ordering π of V , and then swap the positions of random pairs of vertices in π, aiming to reduce the size of the solution X π . We first show how to construct X π from π. Start with X 0 π = ∅, and then define X n π inductively as follows:

X n+1 π = X n π if N + (v n+1 π ) ∩ V n π ⊆ X n π X n π ∪ v n+1 π otherwise Where V n π = {v ∈ V | π(v) ⩽ n} and v n π is the vertex such that π(v n π ) = n. We then set X π = X |V | π .
Notice that X π can be computed in linear time using a traversal of G in the order given by π.

We can use this structure to derive a simple local search algorithm, which repeatedly chooses a uniformly random pair of vertices (u, v) and swaps their positions to obtain a new ordering π ′ such that π ′ (u) = π(v), π ′ (v) = π(u) and π ′ (w) = π(w) for w ̸ ∈ {u, v}. We then compute |X π ′ |, and if |X π ′ | ⩽ |X π |, we validate the swap and set π to π ′ , otherwise we cancel the swap and leave π unchanged. In practice, iterating this process quickly reduces the size of the solution X π on most instances.

Given a solution X obtained after applying the local search of the previous section 3.1, we can compute the initial value of π using a topological sort of the acyclic digraph G[V \X], to which we append a random permutation of X.