Vehicle platooning schemes considering V2V communications: A joint communication/control approach Tiago R. Gonc ¸alves1 , Vineeth S. Varma2 , and Salah E. Elayoubi 1

Abstract-This article addresses communication and control aspects of platooning systems with the related challenges introduced by the overlap of both areas. The main objective is to provide a dynamic control mechanism where the parameters of the well-known Predicted Cooperative Adaptive Cruise Control (PCACC) are adapted based on the observed quality of the V2V (Vehicle-to-Vehicle) communication links. Different from the state of the art, our main design goal is the minimization of intervehicular distances while being robust in terms of an extremely low probability of emergency braking. A new adaptive control scheme based on the offline optimization of the control gains is proposed. We evaluate the new approach in a highway scenario and show the improvements obtained by the dynamic adaptation of the control parameters over static control strategies.

Index Terms-Vehicle platoons, cooperative adaptive cruise control (CACC), wireless communication.

I. INTRODUCTION

During the last decades, several approaches have been developed to deal with autonomous vehicles. Among some of the control technologies that were deployed to help the drivers' safety and increase their driving experience, one can cite the Cruise Control (CC) as a precursor of autonomous cars. However, the first meaningful step to allow the implementation of platoons in the vehicular environment was the establishment of the Adaptive Cruise Control (ACC). It was first introduced by Ioannou and Chien [START_REF] Ioannou | Autonomous intelligent cruise control[END_REF] and it consists of an autonomous control scheme based on constant time headway safety distance. Essentially, it allows autonomously to keep a certain desired distance from the preceding vehicle due to onboard sensors such as camera, radar, and lidar. A more sophisticated approach based on a coordinated exchange of information supported by wireless communication was early provided by Rajamani et al. [START_REF] Rajamani | Design and experimental implementation of longitudinal control for a platoon of automated vehicles[END_REF] and it is known as Cooperative Adaptive Cruise Control (CACC). However, as CACC technology relies on other's vehicle information it is vulnerable to inherent communication aspects such as packet loss and latency. Our particular interest in the platooning system is because it is designed to take advantage of the particular distribution of a convoy in order to increase road capacity and to decreased fuel consumption that is achievable by gathering vehicles close together in order to reduce the air resistance of the platoon's members. The exchange of information is crucial to the deployment of platoons as it allows taking control actions based on the most up-to-date information about the road and traffic status.

We consider in this paper the V2V communication method in platoons as illustrated in Fig. 1. V2V enables communication between vehicles as long as they are in a certain range. This is the case of the IEEE 802.11p and 3GPP C-V2X Mode 4 technologies. The former is an adapted version of the IEEE 802.11a to suit vehicular applications [START_REF] Jiang | IEEE 802.11p: Towards an international standard for wireless access in vehicular environments[END_REF], while the latter is the direct communication interface standardized by the 3GPP. Notice that the network-based communication approach of the 3GPP standard, known as 3GPP C-V2X Mode 3, is not considered here. [START_REF] Molina-Masegosa | LTE-V for sidelink 5G V2X vehicular communications[END_REF] conducted simulations to show a comprehensive analysis of the advantages of the performance of C-V2X Mode 4 over the 802.11p. Similarly, [START_REF] Wang | Comparison of LTE and DSRC-based connectivity for intelligent transportation systems[END_REF] used a simulation environment to compare the communications performance of both modes 3 and 4 of C-V2X with the 802.11p standards. More recently, [START_REF] Gonzalez-Martín | Analytical Models of the Performance of C-V2X Mode 4 Vehicular Communications[END_REF] adopted an analytic approach to describe the C-V2X Mode 4 performance. However, these works were limited to the telecommunication aspects and did not consider the control aspects of the platooning problem.

Another set of works considered the platooning scenario under different communication approaches. For instance, [START_REF] Vinel | Vehicle-to-vehicle communication in C-ACC/platooning scenarios[END_REF] used the 802.11p technology to evaluate the communication performance under a CACC controller in platoons. Likewise, [START_REF] Vukadinovic | 3GPP C-V2X and IEEE 802.11 p for Vehicleto-Vehicle communications in highway platooning scenarios[END_REF] has adopted both wireless technologies (802.11p and C-V2X) and compared their performances in terms of the intervehicular distance of the platoon. [START_REF] Ploeg | Graceful degradation of cooperative adaptive cruise control[END_REF] proposed a control strategy for graceful degradation based on estimating the preceding vehicle's acceleration in case of packet losses, but it mainly deals with extreme cases like complete link failure or lack of a wireless device on one of the vehicles. Different from the aforementioned article, we propose an online adaptation of the control parameter based on the observed quality of the communication link determined by the distance to the transmitter and the level of interference caused by other vehicles. More recently, [START_REF] Sybis | Communication Aspects of a Modified Cooperative Adaptive Cruise Control Algorithm[END_REF] introduced a string stability analysis for a CACC alternative control design, in this paper called Predicted Cooperative Adaptive Cruise Control (PCACC).

The novelty in this paper is the introduction of a dynamic control mechanism where some of the parameters of the PCACC controller are adapted based on information about V2V communication. In particular, we adapt the parameter that is responsible to weigh the influence of the leader's broadcasted messages in the control algorithm, as well as the target distance between vehicles, based on the communication links qualities. Keeping in mind the above discussion and the results available in the literature on cooperative platooning systems, the following are the main contributions of this paper:

1) Evaluation of the robustness of the platooning mechanism under severe conditions for V2V communications, expressed in long bursts of losses and in difficult traffic jamming conditions on the road. 2) Offline optimization of the platooning control parameters based on extensive simulations of a highway scenario. 3) Online adaptation of the control parameters based on the observed communication link quality and on the results of the offline optimization. 4) Adoption of safety as a primary performance metric, quantified in terms of avoiding emergency braking. This translates to robustness constraints, where the intervehicle distance in the platoon is set so that emergency braking is avoided in 99,999% of the cases. The paper is organized as follows. An overview of the control and communication aspects is presented in Section II, and the proposed scheme is described. Section III presents the system model, introduces the robustness case scenario and discusses the performance evaluation. Section IV includes considerations for future work and concludes the paper.

II. CONTROL AND COMMUNICATION PLATOONING SCENARIOS

The objective of this section is to present the control schemes and their interaction with algorithms performance, and to introduce the proposed scheme based on PCACC.

A. Adaptive Cruise Control

The ACC scheme autonomously allows the equipped vehicle to keep a certain desired distance apart from the preceding vehicle. This is possible because of the adoption of cameras, radars, and lidars that measure in real-time the preceding vehicle's position and velocity. Different spacing policies can be used such as Constant Spacing (CS) and Constant Time-Gap (CTG) policy, for instance, see [START_REF] Rajamani | Vehicle dynamics and control[END_REF]. Briefly, the controller on the former scheme aims to keep a constant distance of the preceding vehicle while on the latter one it aims to control the clearance, or time gap, between the host vehicle and the preceding vehicle. In this paper, the CTG policy is applied with the following control law introduced by Ioannou and Chien [START_REF] Ioannou | Autonomous intelligent cruise control[END_REF]

a i des = ẍi des = - 1 h ( εi + λδ i ) (1)
where

ε i = x i -x i-1 + L i-1 (2)
δ i = ε i + h ẋi (3)
are the inter-vehicle spacing and the spacing error, respectively. The index i symbolizes the vehicle index, the first vehicle being numbered 0. x i denotes the position of vehicle i, L i its length and a i its acceleration. h is the time-gap parameter and λ is the design gain parameter. The control input is calculated based on the difference of its own velocity and position with the preceding vehicle, ( ẋi , x i ) and ( ẋi-1 , x i-1 ) respectively. Note that in platooning systems, the ACC control law is always adopted by the leader since it is preceded by a vehicle that is not subject to the platooning controller.

B. Predictive Cooperative Adaptive Cruise Control

PCACC, introduced in [START_REF] Sybis | Communication Aspects of a Modified Cooperative Adaptive Cruise Control Algorithm[END_REF], implies that the control effort, the desired acceleration, of the leader and of the preceding vehicle are available to the following vehicle and its control law is given by

Vehiclei Controlleri Radar V2V xi, ẋi xi-1, ẋi-1 ẍl des, ẋl des, ẍ(i-1) des ẍi des Vehiclei-1 Controlleri-1 Radar V2V xi-1, ẋi-1 xi-2, ẋi-2 ẍ(i-1) des Wireless ẍ(i-
a i des = ẍi des = (1 -C 1 )ẍ (i-1) des + C 1 ẍl des -(2ξ -C 1 (ξ + ξ 2 -1))ω n ǫi -(ξ + ξ 2 -1)ω n C 1 ( ẋi -ẋl ) -ω 2 n ǫ i (4)
where

ǫ i =x i -x i-1 + L i + D des (5) ǫi = ẋi -ẋi-1 . ( 6
)
L i is the length of the vehicle and D des is the desired intervehicular distance that we want to minimize. The control parameters to be tuned are C 1 , ξ and ω n . The parameter C 1 takes on values 0 ≤ C 1 < 1 and is responsible to weigh the contribution of the leader's speed and acceleration. ξ is the controller damping ratio and ω n is the controller bandwidth. The adopted control parameters are shown in Table I. Such improvement of this fully predictive cooperative control is due to the fact that it allows the communication between all the vehicles in the platoon including the leader as shown in Fig. 2. It also illustrates the input signals required by the PCACC controller and the respective segments (V2V radio, vehicle dynamics, and radar equipment) that provide them. In other words, the CACC sends the actual acceleration, which is measured after the actuation lag, while the PCACC is able to propagate the actual values that will become effective after the actuation lag. The PCACC is expected to be superior to CACC because the actuation lag of the system does not affect directly the control effort, which is a big limiting factor for achieving short inter-vehicle distances [START_REF] Vukadinovic | 3GPP C-V2X and IEEE 802.11 p for Vehicleto-Vehicle communications in highway platooning scenarios[END_REF].

C. Semi-Autonomous Control

This is a particular class of control that allows the system to achieve string stability with a constant space control strategy without direct communication from the leading vehicle. The control strategy can be obtained with a simple change in one of the control parameters of the CACC control scheme. Therefore, C 1 = 0 yields

a i des = ẍi des = ẍi-1 -2ξω n ǫi -ω 2 n ǫ i . (7)
The control law is based on one vehicle look-ahead communication topology, which means that only preceding vehicle's information is required. In this case, its acceleration, velocity, and position (ẍ i-1 , ẋi-1 , x i-1 ) are necessary, but without any information from the leader. As proven in [START_REF] Swaroop | Constant spacing strategies for platooning in automated highway systems[END_REF], this particular control law is only string stable when the control effort is always available which means when there is no delay in the process. In the control field, the concept of string stability is recurrently mentioned when it comes to platoon systems. A string stable platoon means that any acceleration or braking in the first vehicle is not going to cause an amplification of the error along the tail of the platoon. In other words, as long as the first vehicle is able to avoid a collision all others will be able too. In [START_REF] Rajamani | Vehicle dynamics and control[END_REF], there is more information about mathematical definitions and conditions to ensure the string stability.

D. Proposed dynamic scheme based on PCACC

In contrast to existing works that assume a fixed control strategy, our main contribution is to propose an approach that will adapt the control parameters based on the communication link quality characterized by the packet error rate from the leader to the last vehicle (defined as PER LLV ). We aim to keep a homogeneous control strategy in the whole platoon. In particular, among the control inputs, the pair (C 1 , D des ) has the most substantial impact on the performance of the system. As mentioned before, C 1 is responsible to weigh the influence of the leader's message in the control algorithm while D des is the desired inter-vehicular distance that we want to set, but due to actuator lag and delay in the process it does not correspond to the actual average inter-vehicular distance (D avg ). The following algorithm based on PCACC controller is proposed.

• Step 1: Update the traffic density range limits (PER LLV ). of Step 3. Therefore, we conduct an offline heuristic optimization to determine the best control parameters (C 1 , D des ), in terms on minimum inter-vehicular distance without collisions, for any given value of PER LLV (which is the result of the traffic density and the resulting interference). We build a (C * 1 (P ER LLV ), D * des (P ER LLV )) lookup table that will serve as an optimum reference for each PER LLV value.

Many existing works on CACC and PCACC mention the minimum distance possible assuming a certain level of interference or traffic on parallel lanes, however, it is not clear what should be done in practice where these parameters will evolve over time in an unknown manner. We propose to apply an online adaptation of the parameters (C 1 , D des ) based on the observed PER LLV and on the results of the offline optimization. Note that due to actuator lag and delay the string stability is not guaranteed for all platoon sizes other than those evaluated.

III. SIMULATION ASSUMPTIONS AND PLATOON MODEL

In this section, we present the system model, the control and communication parameters, and the simulation tool adopted.

A. Vehicle dynamics

The vehicle dynamics is modeled as first-order low pass filter due to the actuator lag. So the transfer function is modeled in the frequency-domain as

G(s) = A i (s) A i des (s) = 1 τ s + 1 (8)
where τ is the time constant of the first-order low pass filter.

A i is the output, which can be interpreted as the actual vehicle acceleration, whereas A i des is the vehicle input, which can be seen as the desired acceleration. Note that, •(s) denotes the Laplace transform of the corresponding time-domain variable.

The idea is to approximate the dynamics of the throttle body and vehicle inertia in order to avoid instantaneous response.

In this paper, we assumed a lag of τ = 0.5 s as in [START_REF] Rajamani | Vehicle dynamics and control[END_REF].

B. Platoon Scenario and Robustness Criterion

The system consists of a platoon of 10 automated vehicles following the leader. The simulated scenario takes place in a four-lane highway, as illustrated in Fig. 1, with a maximum traffic density of 20 interfering vehicles/km. The vehicles of all other lanes are not in platoons (blue vehicles) and a jammer (in red) precedes the platoon leader. In the offline optimization section, the speed of the jammer (vehicle outside of the platoon) follows a preset sequence, adapted to only two cycles from [START_REF] Vukadinovic | 3GPP C-V2X and IEEE 802.11 p for Vehicleto-Vehicle communications in highway platooning scenarios[END_REF], as shown in Fig. 3. The main reason for this jammer velocity profile is to study the capacity of the platoon system to avoid a collision in risky scenarios such as when the vehicle outside of the convoy suddenly applies the maximum brake capacity. Furthermore, different from the previous articles so far that just do Monte-Carlo simulations (100 or 1000 iterations of normal conditions), we have considered a burst of packet losses. While bursts are rare events, they may occur and impact the safety of the platoon, but are not well reflected in the numerical analysis of most previous works like [START_REF] Vinel | Vehicle-to-vehicle communication in C-ACC/platooning scenarios[END_REF]- [START_REF] Cecchini | Performance comparison between IEEE 802.11 p and LTE-V2V in-coverage and out-ofcoverage for cooperative awareness[END_REF]. In this sense, the robustness treated here is related to the following worst-case event: the jammer brakes at some time (t = 60 s in our simulation) and this braking coincides with a burst of packet losses (complete interruption of the transmitted signal) during the following interval t + ∆ also illustrated in Fig. 3.

In order to be conservative, we consider long bursts of packet losses that occur with a probability of 10 -5 . Denoting by PER the probability of packet loss taking into account the channel model and packet collisions and T the time sampling interval for vehicle information, the burst size ∆ (in seconds) can be calculated as

∆ = -5/(log 10 PER)T (9)
In our simulations, we have considered the sampling rate of T = 100 ms as advocated by the ETSI EN 302 637-2 standard.

C. Simulation tool

We used the MATLAB/Simulink environment to model the vehicle dynamics and to implement the control law. Furthermore, we adopted the WLAN Toolbox of MATLAB to implement the channel configuration for a 802.11p transmission in order to obtain the Packet Error Rate (PER) taking into account the V2V fading channel aspects, the additive white Gaussian noise, the packet collisions and others communication parameters as in Table I. The mobility behavior of vehicles is also observed in the MATLAB/Simulink environment as we consider a traffic scenario as in Fig. 1. So at the beginning of each simulation step, from the "Radar" module the "Controller i " is able to update the velocity and the position of the preceding vehicle ( ẋi-1 , x i-1 ) as in Fig. 2. Also, the block "Vehicle i " provides the vehicle's dynamics, thus its own velocity and position ( ẋi , x i ), to the "Controller i " module. Furthermore, the controller reads the leader's acceleration and velocity and the preceding vehicle acceleration (ẍ l , ẋl , ẍi-1 ) from the "V2V" module. Recall that in the leader vehicle we implement an ACC controller, so the "V2V" module is responsible only to broadcast its acceleration and velocity since the ACC controller does not require any other inputs as those provided by the modules "Radar" and "Vehicle". Based on the inputs mentioned, the "Controller i " module is able to calculate the desired acceleration (ẍ i des ) in order to keep a certain desired distance (D des ) from the preceding vehicle. Thus, the "Vehicle i " module, responsible to approximate the dynamics of the vehicle, applies the desired acceleration and provides as output the vehicle's position and velocity at the next simulation step.

The WLAN Toolbox is used to simulate the wireless condition for a 802.11p transmission. We adopted two different packet error rate parameters. The first one is related to the packet error rate between two successive vehicles and defined as PER i,i+1 along the article. We adopted an Highway lineof-sight (LOS) profile as vehicles in the platoon are close enough to justify the usage. The second one is the packet error rate between the leader to the last vehicle that was defined as PER LLV . In this case, we used the Highway non line-of-sight (NLOS) profile in the toolbox since in this case the leader is less likely to be in the LOS with the last vehicle in the platoon. Furthermore, we stochastically calculated the SINR (Signal to Interference and Noise Ratio) for each of the vehicles in the platoon considering Winner-II Path Loss Model (B1 scenario) [START_REF]1.1.2 V1.2, WINNER II channel models[END_REF] and the interference from other vehicles.

IV. PERFORMANCE EVALUATION

The leader is equipped with an ACC control to be in accordance with the recommended safety time interval gap of the respective local law while the platoon is equipped with fully predictive cooperative control. The focus is to apply a longitudinal control in the platoon through V2V communication technology and analyze the system stability by means of vehicle collisions in some robust and worst-case scenarios. We considered the zero-order hold mechanism as the holding strategy for the control signal during the periods of packet losses. Furthermore, in all simulations, we focus on obtaining the minimum feasible inter-vehicular distance in the platoon with a emergency breaking probability no more than 10 -5 . Note that we implemented a safety gap distance of 0.5 m for the emergency braking actuation to avoid collisions in practical settings.

The control strategy demands relative position and longitudinal velocity of the preceding vehicle so we assumed that the measurements are sampled each 60 ms with 1 ms process delay and done by a long-range radar as in [START_REF] Vukadinovic | 3GPP C-V2X and IEEE 802.11 p for Vehicleto-Vehicle communications in highway platooning scenarios[END_REF]. All the vehicles in the platoon broadcast a 500 bytes message on a 10 MHz channel bandwidth. Neighbouring vehicles are subject to a Highway LOS channel model [START_REF] Alexander | Cooperative intelligent transport systems: 5.9-GHz field trials[END_REF] and our simulations provide a PER i,i+1 ∈ {0.006, 0.0245} for the cases without interference and with interference, respectively. We have adopted PER i,i+1 = 0.0245 as the default value. As of the leader communication, the leader broadcasts a message to all other vehicles that is subject to a Highway NLOS channel [START_REF] Alexander | Cooperative intelligent transport systems: 5.9-GHz field trials[END_REF] with PER LLV ∈ {0.2, . . . , 0.7}, depending on the interference conditions.

A. Offline optimization

We start by performing an offline optimization of the control parameters (C 1 , D des ). Fig. 4 illustrates the substantial impact of C 1 parameter on the average inter-vehicular distance (D avg ) for different PER LLV values. Thus, from Fig. 4 we can retrieve the C * 1 optimum that minimizes the inter-vehicular distance for each PER LLV evaluated, which is illustrated in Fig. 5 on the left. Similarly, the optimum D des values were established as shown in Fig. 5 on the right. In this figure, C 1 = 0.3 represents the best alternative in the robust scenario for low PER LLV values. While, C 1 = 0.2 indicates to be the best value for mid-range values as 0.2 ≤ PER LLV ≤ 0.4. For higher PER LLV ≥ 0.5 the best parameter values is C 1 = 0 which is the case of semi-autonomous control.

B. Online adaptation of the control parameters

We now move to the online adaptation of the control parameters, where the whole platoon adapts the control inputs (C 1 , D des ) corresponding to the optimal values computed in the offline optimization based on the packet loss observed on the communication link by the leader and the last vehicle (PER LLV ). For comparison purposes, we also simulated cases with fixed control parameters. The idea is to inspect the occurrence of collisions and to compare the inter-vehicular distance of the platoon in a long simulation of 25 minutes in three different cases defined as follows:

• Case 1 -PCACC with fixed C 1 = 0.2. In this case, we apply the D des correspondent to the average PER LLV . We aim to treat the case where the system operates in a moderate scenario where it's not too conservative.

• Case 2 -Semi-autonomous (C 1 = 0). In order to have a robust controller, we apply the D des corresponding to the worst-case PER LLV .

• Case 3 -The proposed dynamic scheme where the control parameters are based on the observed PER LLV . In all cases, we considered the jammer profile as the pattern from Fig. 3 repeated 50 times. Another important factor is the traffic density that generates interference and changes the PER. In order to have a fair comparison between all cases, we apply a predefined traffic density in the simulation. The range of the traffic densities varies from 0 vehicle/km to 20 vehicle/km Notice that the traffic density changes linearly and less frequently than the jammer incidents. We applied the burst as in ( 9) in each 6 min of simulation in the 9 th vehicle in the platoon to simulate the burst of packet losses in the most critical moment.

From Table II, we notice that Case 1 exhibited 8 collisions while cases 2 and 3 had none. The former case comprises the event where the system is assumed to operate in an average traffic condition such as PER LLV = 0.4. However, it can be seen that this is not a safe approach since it does not guarantee a secure outcome. Case 2 corresponds to the case where the optimum (C * 1 , D * des ) control parameters are set for higher PER LLV values (worst case). Note that this corresponds to the adoption of the semi-autonomous control as C 1 = 0. The idea behind this scenario is to have a robust and safe outcome. Despite no collisions, it exhibited an increase of 21% in the inter-vehicular distance when compared to the suggested Case 3 approach. Fig. 6 shows the average inter-vehicular distance and the correspondent PER LLV value of the suggested approach for the considered simulation. Therefore, the proposed method is demonstrated to be the best option so that platoon remains safe and robust while reducing the inter-vehicular distance.

V. CONCLUDING REMARKS This paper studies the design of the platoon control algorithm based on the V2V communication quality. We started by devising the optimal parameters of the controller for different communication qualities, namely the weight given to the information broadcast by the platoon leader (C 1 parameter), and the desired distance (D des ) between vehicles. We then proposed a new dynamic approach based on the offline optimization of the control parameters (C 1 , D des ). In this dynamic scheme, the quality of the communication link is continuously monitored and the control parameters updated accordingly based on the results of the offline optimization. Our simulation results show that, if the control parameters are not adapted to the channel quality, the semi-autonomous control performs best. However, with the proposed adaptive control, using leader information results in a better performance. In future work, combining different wireless technologies, such as V2N (Vehicleto-Network) and VLC (Visible Light Communication), is a promising means for enhancing the robustness of the platoon while reducing the inter-vehicle distance.

V2V

Fig. 1 .

 1 Fig. 1. Traffic scenario including a platoon with V2V communication approach.

Fig. 3 .

 3 Fig.3. Illustration of the jammer velocity profile adopted with 2 cycles.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Illustration of the average inter-vehicular distance with PER i,i+1 = 0.0245 for different C 1 values.

  PER LLV = 0.2 and PER LLV = 0.7 for a 11 platoon vehicle in the considered scenario, respectively. The predefined traffic densities starts from PER LLV = 0.2 and during each 2 minutes (4 cycles) it adds 0.1 in the PER LLV up to the maximum PER LLV = 0.7. Then, with the same logic, it decreases to the minimum PER LLV = 0.2 and then increases again until PER LLV = 0.4 where its elapsed time reaches 25 min.

Fig. 6 .

 6 Fig. 6. Illustration of the average inter-vehicular distance of the Case 3 simulation and the correspondent boundaries of PER LLV .

TABLE I COMMUNICATION

 I AND CONTROLLER PARAMETERS

	Communication	Controller
	Parameter	Value	Parameter	Value
	Channel	Highway NLOS	Leader factor (C1)	Adaptive
	Path loss	Winner+B1 LOS	Desired distance (D des )	Adaptive
	Noise power	-174 dBm/Hz	Damping ratio (ξ)	2
	Tx power	22.5 dBm	Bandwidth (ωn)	0.5 Hz
	MCS	QPSK, R=1/2	Time gap (h)	1.4 s
	CAM size	500 bytes	Vehicle length (L)	16.5 m
	CAM interval	100 ms	Actuator lag (τ )	500 ms
	Radar interval	60 ms	Lane width	5 m
	Process delay	1 ms

  1) des ẍl des, ẋl des, ẍ(i-2) des

		ẍl des, ẋl des
	ẍl des,	ẍl des,
	ẋl des	ẋl des

ẍ(i-1) des . . . Wireless . . .

ẍi des Wireless ẍ(i-2) des . . . ẍ(i-2) des Fig. 2. Block diagram of the platoon system with a PCACC control between vehicle i -1 and vehicle i.

  Vary the C 1 parameter while minimizing D des and register the average inter-vehicular distance D avg for each PER LLV inside the range of Step 1.

• Step 2: • Step 3: Consider the minimum D avg result with no collision and identify its correspondent pair of (C 1 (P ER LLV ), D des (P ER LLV )) to build the optimum lookup table (C * 1 (P ER LLV ), D * des (P ER LLV )). • Step 4: Observe the current communication link and adapt the control inputs (C 1 , D des ) accordingly based on the optimum lookup table

T. R. Gonc ¸alves, and S. E. Elayoubi are with Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des signaux et systèmes, 91190, Gif-sur-Yvette, France. E-mail: {tiago.rochagoncalves; salaheddine.elayoubi}@centralesupelec.fr

. 2 V. S. Varma is with Université de Lorraine, CRAN, CNRS, UMR 7039, France. E-mail: vineeth.satheeskumar-varma@univ-lorraine.fr.


Tuan-Binh Nguyen
email: tuan-binh.nguyen@inria.fr

Jérôme-Alexis Chevalier

Bertrand Thirion

Sylvain Arlot

Aggregation of Multiple Knockoffs

. This new method, called aggregation of multiple knockoffs (AKO), addresses the instability inherent to the random nature of knockoff-based inference. Specifically, AKO improves both the stability and power compared with the original knockoff algorithm while still maintaining guarantees for false discovery rate control. We provide a new inference procedure, prove its core properties, and demonstrate its benefits in a set of experiments on synthetic and real datasets.

INTRODUCTION

In many fields, multivariate statistical models are used to fit some outcome of interest through a combination of measurements or features. For instance, one might predict the likelihood for individuals to declare a certain type of disease based on genotyping information. Besides prediction accuracy, the inference problem consists in defining which measurements carry useful features for prediction. More precisely, we aim at conditional inference (as opposed to marginal inference), that is, analyzing which features carry information given the other features. This inference is however very challenging in high-dimensional settings.

Among the few available solutions, knockoff-based (KO) inference [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF][START_REF] Candès | Panning for gold: 'model-X' knockoffs for high dimensional controlled variable selection[END_REF] consists in introducing noisy copies of the original variables that are independent from the outcome conditional on the original variables, and comparing the coefficients of the original variables to those of the knockoff variables. This approach is particularly attractive for several reasons: i) it is not tied to a given statistical model, but can work instead for many different multivariate functions, whether linear or not; ii) it requires a good generative model for features, but poses few conditions for the validity of inference; and iii) it controls the false discovery rate (FDR, [START_REF] Benjamini | Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing[END_REF]), a more useful quantity than multiplicity-corrected error rates.

Unfortunately, KO has a major drawback, related to the random nature of the knockoff variables: two different draws yield two different solutions, leading to large, uncontrolled fluctuations in power and false discovery proportion across experiments (see Figure 1 below). This makes the ensuing inference irreproducible. An obvious way to fix the problem is to rely on some type of statistical aggregation, in order to consolidate the inference results. Such procedures have been introduced by [START_REF] Gimenez | Improving the Stability of the Knockoff Procedure: Multiple Simultaneous Knockoffs and Entropy Maximization[END_REF] and by [START_REF] Emery | Controlling the FDR in variable selection via multiple knockoffs[END_REF], but they have several limitations: the computational complexity scales poorly with the number B of bootstraps, while the power of the method decreases with B. In highdimensional settings that we target, these methods are thus only usable with a limited number of bootstraps.

In this work, we explore a different approach, that we call aggregation of multiple knockoffs (AKO): it rests on a reformulation of the original knockoff procedure that introduces intermediate p-values. As it is possible to aggregate such quantities even without assuming independence [START_REF] Meinshausen | p-Values for High-Dimensional Regression[END_REF], we propose to perform aggregation at this intermediate step. We first establish the equivalence of AKO with the original knockoff aggregation procedure in case of one bootstrap (Proposition 1). Then we show that the FDR is also controlled with AKO (Theorem 1). By construction, AKO is more stable than (vanilla) knockoff; we also demonstrate empirical benefits in several examples, using simulated data, but also genetic and brain imaging data. Note that the added knockoff generation and inference steps are embarrassingly parallel, making this procedure no more costly than the original KO inference.

Notation. Let [p] denote the set {1, 2, . . . , p}; for a given set given set A, |A| ∆ = card(A); matrices are denoted in bold uppercase letter, while vectors in bold lowercase letter and scalars normal character. An exception for this is the vector of knockoff statistic W, in which we follow the notation from the original paper of [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF].

BACKGROUND

Problem Setting. Let X ∈ R n×p be a design matrix corresponding to n observations of p potential explanatory variables x 1 , x 2 , . . . , x n ∈ R p , with its target vector y ∈ R n . To simplify the exposition, we focus on sparse linear models, as [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF] and [START_REF] Candès | Panning for gold: 'model-X' knockoffs for high dimensional controlled variable selection[END_REF]:

(1) y = Xβ * + σ where β * ∈ R p is the true parameter vector, σ ∈ R + the unknown noise magnitude, ∈ R n some Gaussian noise vector. Yet, it should be noted that the algorithm does not require linearity or sparsity. Our main interest is in finding an estimate S of the true support set S = {j ∈ [p] : β * j = 0}, or the set of important features that have an effect on the response. As a consequence, the complementary of the support S, which is denoted S c = {j ∈ [p] : β * j = 0}, corresponds to null hypotheses. Identifying the relevant features amounts to simultaneously testing

H j 0 : β * j = 0 versus H j a : β * j = 0, ∀j = 1, . . . , p.
Specifically, we want to bound the proportion of false positives among selected variables, that is, control the false discovery rate (FDR, [START_REF] Benjamini | Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing[END_REF]) under certain predefined level α:

FDR = E | S ∩ S c | | S| ∨ 1 ≤ α ∈ (0, 1) .
Knockoff Inference. Introduced originally by [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF], the knockoff filter is a variable selection method for multivariate models with theoretical control of FDR. [START_REF] Candès | Panning for gold: 'model-X' knockoffs for high dimensional controlled variable selection[END_REF] expanded the method to work in the case of (mildly) high-dimensional data, with the assumption that x = (x 1 , . . . , x p ) ∼ P X such that P X is known. The first step of this procedure involves sampling extra null variables that have a correlation structure similar to that of the original variables, with the following formal definition.

Definition 1 (Model-X knockoffs, [START_REF] Candès | Panning for gold: 'model-X' knockoffs for high dimensional controlled variable selection[END_REF]). The model-X knockoffs for the family of random variables x = (x 1 , . . . , x p ) are a new family of random variables x = (x 1 , . . . , xp ) constructed to satisfy the two properties:

1. For any subset K ⊂ {1, . . . , p}, (x, x) swap(K) d = (x, x), where the vector (x, x) swap(K)

denotes the swap of entries x j and xj for all j ∈ K, and

d = denotes equality in distribution. 2. x ⊥ ⊥ y | x .
A test statistic is then calculated to measure the strength of the original variables versus their knockoff counterpart. We call this the knockoff statistic W = {W j } p j=1 , that must fulfill two important properties.

Definition 2 (Knockoff statistic, [START_REF] Candès | Panning for gold: 'model-X' knockoffs for high dimensional controlled variable selection[END_REF]). A knockoff statistic W = {W j } j∈[p] is a measure of feature importance that satisfies the two following properties:

1. It depends only on X, X and y W = f (X, X, y) .

2. Swapping the original variable column x j and its knockoff column xj switches the sign of W j :

W j ([X, X] swap(S) , y) = W j ([X, X], y) if j ∈ S c -W j ([X, X], y) if j ∈ S .
Following previous works on the analysis of the knockoff properties [ACC17, RRJW20], we make the following assumption about the knockoff statistic. This is necessary for our analysis of knockoff aggregation scheme later on.

Assumption 1 (Null distribution of knockoff statistic). The knockoff statistic defined in Definition 2 are such that {W j } j∈S c , are independent and follow the same distribution

P 0 .
Remark 1. As a consequence of [CFJL18, Lemma 2] regarding the signs of the null W j as i.i.d. coin flips, if Assumption 1 holds true, then P 0 is symmetric around zero.

One such example of knockoff statistic is the Lasso-coefficient difference (LCD). The LCD statistic is computed by first making the concatenation of original variable and knockoff variables [X, X] ∈ R n×2p , then solving the Lasso problem [Tib96]:

(2)

β = argmin β∈R 2p 1 2 y -[X, X]β 2 2 + λ β 1
with λ ∈ R the regularization parameter, and finally to take:

(3) ∀j ∈ [p] , W j = | β j | -| β j+p | .
This quantity measures how strong the coefficient magnitude of each original covariate is against its knockoff, hence the name Lasso-coefficient difference. Clearly, the LCD statistic satisfies the two properties stated in Definition 2.

Finally, a threshold for controlling the FDR under given level α ∈ (0, 1) is calculated:

(4) τ + = min t > 0 : 1 + #{j : W j ≤ -t} #{j : W j ≥ t} ∨ 1 ≤ α ,
and the set of selected variables is S = {j ∈ [p] :

W j ≥ τ + }.
Instability in Inference Results. Knockoff inference is a flexible method for multivariate inference in the sense that it can use different loss functions (least squares, logistic, etc.), and use different variable importance statistics. However, a major drawback of the method comes from the random nature of the knockoff variables X obtained by sampling: different draws yield different solutions (see Figure 1 in Section 5.1). This is a major issue in practical settings, where knockoff-based inference is used to prove the conditional association between features and outcome.

AGGREGATION OF MULTIPLE KNOCKOFFS

Algorithm Description

One of the key factors that lead to the extension of the original (vanilla) knockoff filter stems from the observation that knockoff inference can be formulated based on the following quantity.

Definition 3 (Intermediate p-value). Let W = {W j } j∈[p] be a knockoff statistic according to Definition 2. For j = 1, . . . , p, the intermediate p-value π j is defined as:

(5)

π j =      1 + #{k : W k ≤ -W j } p if W j > 0 1 if W j ≤ 0 .
We first compute B draws of knockoff variables, and then knockoff statistics. Using Eq. ( 5), we derive the corresponding empirical p-values π

(b) j , for all j ∈ [p] and b ∈ [B]
. Then, we aggregate them for each variable j in parallel, using the quantile aggregation procedure introduced by [MMB09]:

(6) πj = min    1, q γ {π (b) j : b ∈ [B]} γ   
where q γ (•) is the γ-quantile function. In the experiments, we fix γ = 0.3 and B = 25. The selection of these default values is explained more thoroughly in Section 5.1.

Finally, with a sequence of aggregated p-values π1 , . . . , πp , we use Benjamini-Hochberg step-up procedure (BH, [START_REF] Benjamini | Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing[END_REF]) to control the FDR.

Definition 4 (BH step-up, [START_REF] Benjamini | Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing[END_REF]). Given a list of p-values π1 , . . . , πp and predefined FDR control level α ∈ (0, 1), the Benjamini-Hochberg step-up procedure comprises three steps:

1. Order p-values such that: π(1) ≤ π(2) ≤ • • • ≤ π(p) .
2. Find:

(7) k BH = max k : π(k) ≤ kα p . 3. Select S = {j ∈ [p] : π(j) ≤ π( k BH ) }.
This procedure controls the FDR, but only under independence or positive-dependence between p-values [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF]. As a matter of fact, for a strong guarantee of FDR control, one can consider instead a threshold yielding a theoretical control of FDR under arbitrary dependence, such as the one of [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF]. We call BY step-up the resulting procedure. Yet we use BH step-up procedure in the experiments of Section 5, as we observe empirically that the aggregated p-values πj defined in Eq. (5) does not deviate significantly from independence (details in Appendix).

Definition 5 (BY step-up, [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF]). Given an ordered list of p-values as in step 1 of BH step-up π(1) ≤ π(2) ≤ • • • ≤ π(p) and predefined level α ∈ (0, 1), the Benjamini-Yekutieli step-up procedure first finds:

(8) k BY = max k ∈ [p] : π(k) ≤ kβ(p)α p , with β(p) = ( p i=1 1/i) -1
, and then selects

S = j ∈ [p] : π(j) ≤ π( k BY ) .
[BR09] later on introduced a general function form for β(p) to make BY step-up more flexible. However, because we always have β(p) ≤ 1, this procedure leads to a smaller threshold than BH step-up, thus being more conservative.

Algorithm 1 AKO -Aggregation of multiple knockoffs

Input: X ∈ R n×p , y ∈ R n , B -number of bootstraps ; α ∈ (0, 1) -target FDR level Output: SAKO -Set of selected variables index

for b = 1 to B do X(b) ← sampling_knockoff(X) W (b) ← knockoff_statistic(X, X(b) , y) π (b) ← convert_statistic(W (b) ) //
Return: SAKO ← j ∈ [p] : πj ≤ π k
The AKO procedure is summarized in Algorithm 1. We show in the next section that with the introduction of the aggregation step, the procedure offers a guarantee on FDR control under mild hypotheses. Additionally, the numerical experiments of Section 5 illustrate that aggregation of multiple knockoffs indeed improves the stability of the knockoff filter, while bringing significant statistical power gains.

Related Work

To our knowledge, up until now there have been few attempts to stabilize knockoff inference. Earlier work of [START_REF] Su | Communication-efficient false discovery rate control via knockoff aggregation[END_REF] rests on the same idea of generating multiple knockoff bootstrap as ours, but relies on the linear combination of the so-called one-bit p-values (introduced as a means to prove the FDR control in original knockoff work of [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF]). As such, the method is less flexible since it requires a specific type of knockoff statistic to work. Furthermore, it is unclear how this method would perform in high-dimensional settings, as it was only demonstrated in the case of n > p. More recently, the work of [START_REF] Holden | Multiple Model-Free Knockoffs[END_REF] incorporates directly multiple bootstraps of knockoff statistics for FDR thresholding without the need of p-value conversion. Despite its simplicity and convenience as a way of aggregating knockoffs, our simulation study in Section 5.1 demonstrates that this method somehow fails to control FDR in several settings.

In a different direction, [START_REF] Gimenez | Improving the Stability of the Knockoff Procedure: Multiple Simultaneous Knockoffs and Entropy Maximization[END_REF] and [START_REF] Emery | Controlling the FDR in variable selection via multiple knockoffs[END_REF] have introduced simultaneous knockoff procedure, with the idea of sampling several knockoff copies at the same time instead of doing the process in parallel as in our work. This, however, induces a prohibitive computational cost when the number of bootstraps increases, as opposed to the AKO algorithm that can use parallel computing to sample multiple bootstraps at the same time. In theory, on top of the fact that sampling knockoffs has cubic complexity on runtime with regards to number of variables p (requires covariance matrix inversion), simultaneous knockoff runtime is of O(B 3 p 3 ), while for AKO, runtime is only of O(Bp 3 ) and O(p 3 ) with parallel computing. Moreover, the FDR threshold of simultaneous knockoff is calculated in such a way that it loses statistical power as the number of bootstraps increases, when the sampling scheme of vanilla knockoff by [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF] is used. We have set up additional experiments in the Appendix to illustrate this phenomenon. In addition, the threshold introduced by [EK19] is only proven to have a theoretical control of FDR in the case where n > p.

THEORETICAL RESULTS

We now state our theoretical results about the AKO procedure.

Equivalence of Aggregated Knockoff with Single Bootstrap

First, when B = 1 and γ = 1, we show that AKO+BH is equivalent to vanilla knockoff.

Proposition 1 (Proof in Appendix A.1). Assume that for all j, j = 1, . . . , p,

P(W j = W j , W j = 0, W j = 0) = 0
that is, non-zero LCD statistics are distinct with probability 1. Then, single bootstrap version of aggregation of multiple knockoffs (B = 1), using γ = 1 and BH step-up procedure in Definition 4 for calculating FDR threshold, is equivalent to the original knockoff inference by [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF].

Remark 2. Although Proposition 1 relies on the assumption of distinction between non-zero W j s for all j = 1, . . . , p, the following lemma establishes that this assumption holds true with probability one for the LCD statistic up to further assumptions.

Lemma 1 (Proof in Appendix A.2). Define the equi-correlation set as:

J λ = j ∈ [p] : x j (y -X β) = λ/2
with β, λ defined in Eq. ( 2). Then we have:

(9) P W j = W j , W j = 0, W j = 0, rank(X J λ ) = | J λ | = 0
for all j, j ∈ [p] : j = j . In other words, assuming X J λ is full rank, then the event that LCD statistic defined in Eq. ( 3) is distinct for all non-zero value happens almost surely.

Validity of Intermediate P-values

Second, the fact that the π j are called "intermediate p-values" is justified by the following lemma.

Lemma 2. If Assumption 1 holds true, and if |S c | ≥ 2, then, for all j ∈ S c , the intermediate p-value π j defined by Eq. (5) satisfies:

∀t ∈ [0, 1] P(π j ≤ t) ≤ κp |S c | t where κ = √ 22 -2 7 √ 22 -32 ≤ 3.24.
Proof. The result holds when t ≥ 1 since κp ≥ p ≥ |S c | and a probability is always smaller than 1. Let us now focus on the case where t ∈ [0, 1), and define m = |S c |-1 ≥ 1 by assumption. Let F 0 denote the c.d.f. of P 0 , the common distribution of the null statistics {W k } k∈S c , which exists by Assumption 1. Let j ∈ S c be fixed. By definition of π j , when W j > 0 we have:

π j = 1 + #{k ∈ [p] : W k ≤ -W j } p = 1 + #{k ∈ S : W k ≤ -W j } p + #{k ∈ S c \ {j} : W k ≤ -W j } p (since W j > 0 > -W j ) ≥ m p F m (-W j ) + 1 p (10) where ∀u ∈ R, F m (u) ∆ = #{k ∈ S c \ {j} : W k ≤ u} m is the empirical cdf of {W k } k∈\{j} .
Therefore, for every t ∈ [0, 1),

P(π j ≤ t) = P(π j ≤ t and W j > 0) + P(π j ≤ t and W j ≤ 0) =0 since π j =1 when W j ≤0 = E P(π j ≤ t | W j )1 W j >0 ≤ E P m p F m (-W j ) + 1 p ≤ t | W j 1 W j >0 by (10) ≤ P m p F m (-W j ) + 1 p ≤ t . (11)
Notice that W j has a symmetric distribution around 0, as shown by Remark 1, that is, -W j and W j have the same distribution. Since W j and {W k } k∈S c \{j} are independent with the same distribution P 0 by Assumption 1, they have the same joint distribution as

F -1 0 (U ), F -1 0 (U 1 ), . . . , F -1 0 (U m )
where U, U 1 , . . . , U m are independent random variables with uniform distribution over [0, 1], and F -1 0 denotes the generalized inverse of F 0 . Therefore, Eq. ( 11) can be rewritten as

P(π j ≤ t) ≤ P m p F m F -1 0 (U ) + 1 p ≤ t (12) where ∀v ∈ R, F m (v) ∆ = 1 m m k=1 1 F -1 0 (U k )≤v .
Notice that for every u ∈ R,

G m (u) ∆ = 1 m m k=1 1 U k ≤u ≤ 1 m m k=1 1 F -1 0 (U k )≤F -1 0 (u) = F m F -1 0 (u) since F -1 0
is non-decreasing. Therefore, Eq. ( 12) shows that

P(π j ≤ t) ≤ P m G m (U ) ≤ tp -1 = 1 0 P m G m (u) ≤ tp -1 du . (13)
Now, we notice that for every u ∈ (0, 1), m G m (u) follows a binomial distribution with parameters (m, u). So, a standard application of Bernstein's inequality [BLM13, Eq. 2.10] shows that for every 0

≤ x ≤ u ≤ 1, P m G m (u) ≤ mx ≤ exp -m 2 (u -x) 2 2mu + m(u-x) 3 = exp -3mx u x -1 2 7u x -1 .
Note that for every λ ∈ (0, 1/7), we have

∀w ≥ 1 -λ 1 -7λ ≥ 1 , w -1 7w -1 ≥ λ hence ∀u ≥ x 1 -λ 1 -7λ , P m G m (u) ≤ mx ≤ exp -3mλx u x -1 .
As a consequence, ∀λ ∈ (0, 1/7),

1 0 P m G m (u) ≤ mx du ≤ 1 -λ 1 -7λ x + 1 1-λ 1-7λ x exp -3mλ(u -x) du ≤ 1 -λ 1 -7λ x + +∞ 6λ 1-7λ x exp(-3mλv)dv ≤ 1 -λ 1 -7λ x + 1 3mλ exp -3mλ 6λ 1 -7λ x ≤ 1 -λ 1 -7λ x + 1 3mλ
.

Taking x = (tp -1)/m, we obtain from Eq. ( 13) that ∀λ ∈ (0, 1/7)

P(π j ≤ t) ≤ 1 -λ 1 -7λ tp -1 m + 1 3mλ = 1 -λ 1 -7λ tp m + 1 3λ - 1 -λ 1 -7λ 1 m . ( 14
)
Choosing λ = (5 -√ 22)/3 ∈ (0, 1/7), we have 1 3λ = 1-λ 1-7λ hence the result with

κ = 1 -λ 1 -7λ = √ 22 -2 7 √ 22 -32 ≤ 3.24.
Remark 3. If the definition of π j is replaced by

(15) π j,c ∆ =      c + #{k : W k ≤ -W j } p if W j > 0 1 if W j ≤ 0
for some c > 0, the above proof also applies and yields an upper bound of the form

∀t ≥ 0 , P(π j,c ≤ t) ≤ κ(c)t
for some constant κ(c) > 0. It is then possible to make κ(c) as close to 1 as desired, by choosing c large enough. Lemma 2 corresponds to the case c = 1.

Note that we also prove in the Appendix that if p → +∞ with |S| p, then for every j ≥ 1 such that β * j = 0, π j is an asymptotically valid p-value, that is,

(16) ∀t ∈ [0, 1] , lim sup p→+∞ P(π j ≤ t) ≤ t .
Yet, proving our main result (Theorem 1) requires a non-asymptotic bound such that the one of Lemma 2.

FDR control for AKO

Finally, the following theorem provides a non-asymptotic guarantee about the FDR of AKO with BY step-up.

Theorem 1. If Assumption 1 holds true and |S c | ≥ 2, then for any B ≥ 1 and α ∈ (0, 1), the output S AKO+BY of aggregation of multiple knockoff (Algorithm 1), with the BY step-up procedure, has a FDR controlled as follows:

E | S AKO+BY ∩ S c | | S AKO+BY | ∨ 1 ≤ κα
where κ ≤ 3.24 is defined in Lemma 2.

Sketch of the proof. The proof of [MMB09, Theorem 3.3], which itself relies partly on [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF], can directly be adapted to upper bound the FDR of S AKO+BY in terms of quantities of the form P(π (b) j ≤ t) for j ∈ S c and several t ≥ 0. Combined with Lemma 2, this yields the result. A full proof is provided in Appendix A.5.

Note that Theorem 1 loses a factor κ compared to the nominal FDR level α. This can be solved by changing α into α/κ in the definition of S AKO+BY . Nevertheless, in our experiments, we do not use this correction and find that the FDR is still controlled at level α.

EXPERIMENTS

Compared Methods. We make benchmarks of our proposed method aggregation of multiple knockoffs (AKO) with B = 25, γ = 0.3 and vanilla knockoff (KO), along with other recent methods for controlling FDR in high-dimensional settings, mentioned in Section 3.2: simultaneous knockoff, an alternative aggregation scheme for knockoff inference introduced by [GZ19] (KO-GZ), along with its variant of [START_REF] Emery | Controlling the FDR in variable selection via multiple knockoffs[END_REF] (KO-EK); the knockoff statistics aggregation by [START_REF] Holden | Multiple Model-Free Knockoffs[END_REF] (KO-HH); and debiased Lasso (DL-BH) [START_REF] Javanmard | False discovery rate control via debiased lasso[END_REF].

Synthetic Data

Simulation Setup. Our first experiment is a simulation scenario where a design matrix X (n = 500, p = 1000) with its continuous response vector y are created following a linear model assumption. The matrix is sampled from a multivariate normal distribution of zero mean and covariance matrix Σ ∈ R p×p . We generate Σ as a symmetric Toeplitz matrix that has the structure:

Σ =        ρ 0 ρ 1 . . . ρ p-1 ρ 1 . . . . . . ρ p-2 . . . . . . . . . . . . ρ p-1 ρ p-2 . . . ρ 0       
where the ρ ∈ (0, 1) parameter controls the correlation structure of the design matrix. This means that neighboring variables are strongly correlated to each other, and the correlation decreases with the distance between indices. The true regression coefficient β * vector is picked with a sparsity parameter that controls the proportion of nonzero elements with amplitude 1. The noise is generated to follow N (µ, I n ) with its magnitude σ = Xβ * 2 /(SNR 2 ) controlled by the SNR parameter. The response vector y is then sampled according to Eq. (1). In short, the three main parameters controlling this simulation are correlation ρ, sparsity degree k and signal-to-noise ratio SNR.

Aggregation Helps Stabilizing Vanilla Knockoff.

To demonstrate the improvement in stability of the aggregated knockoffs, we first do multiple runs of AKO and KO with α = 0.05 under one simulation of X and y. In order to guarantee a fair comparison, we compare 100 runs of AKO, each with B = 25 bootstraps, with the corresponding 2500 runs of KO. We then plot the histogram of FDP and power in Figure 1. For the original knockoff, the false discovery proportion varies widely and has a small proportion of FDP above 0.2 = 4α. Besides, a fair amount of KO runs returns null power.

On the other hand, AKO not only improves the stability in the result for FDP -the FDR being controlled at the nominal level α = 0.05-but it also improves statistical power: in particular, it avoids catastrophic behavior (zero power) encountered with KO.

Inference Results on Different Simulation Settings.

To observe how each algorithm performs under various scenarii, we vary each of the three simulation parameters while keeping the others unchanged at default value. The result is shown in Figure 2. Compared with KO, AKO improves statistical power while still controlling the FDR. Noticeably, in the case of very high correlation between nearby variables (ρ > 0.7), KO suffers from a drop in average power. The loss also occurs, but is less severe for AKO.

Moreover, compared with simultaneous knockoff (KO-GZ), AKO gets better control for FDR and a higher average power in the extreme correlation (high ρ) case. Knockoff statistics aggregation (KO-HH), contrarily, is spurious: it detects numerous truly significant variables with high average statistical power, but at a cost of failure in FDR control, especially when the correlation parameter ρ gets bigger than 0.6. Debiased Lasso (DL-BH) and KO-EK control FDR well in all scenarii, but are the two most conservative procedures. Choice of B and γ for AKO. Figure 3 shows an experiment when varying γ and B. FDR and power are averaged across 30 simulations of fixed parameters: SNR=3.0, ρ = 0.7, sparsity=0.06. Notably, it seems that there is no further gain in statistical

GWAS on Flowering Phenotype of Arabidopsis thaliana

To test AKO on real datasets, we first perform a genome-wide association study (GWAS) on genomic data. The aim is to detect association of each of 174 candidate genes with a phenotype FT_GH that describes flowering time of Arabidopsis thaliana, first done by [AHV + 10]. Preprocessing is done similarly to [AGS + 13]: 166 data samples of 9938 binary SNPs located within a ±20-kilobase window of 174 candidate genes that have been selected in previous publications as most likely to be involved in flowering time traits. Furthermore, we apply the same dimension reduction by hierarchical clustering as [START_REF] Slim | kernelPSI: a Post-Selection Inference Framework for Nonlinear Selection[END_REF] to make the final design matrix of size n = 166 samples × p = 1560 features. We list the detected genes from each method in Table 1.

The three methods that rely on sampling knockoff variables detect AT2G21070. This gene, which is responsible for the mutant FIONA1, is listed by [KKY + 08] to be vital for regulating period length in the Arabidopsis circadian clock. FIONA1 also appears to be involved in photoperiod-dependent flowering and in daylength-dependent seedling growth. In particular, the time for opening of the first flower for FIONA1 mutants are shorter than the ones without under both long and short-day conditions. In addition to FIONA1 mutant, AKO also detects AT4G02780 and AT5G47640. It can be found in studies dating back to the 90s [SCS98] that AT4G02780 encodes a mutation for late flowering. Meanwhile, AT5G47640 mutant delay flowering in long-day but not in short-day experiments [CBE + 07].

Functional Magnetic Resonance Imaging (fMRI) analysis on Human Connectome Project Dataset

Human Connectome Project (HCP900) is a collection of neuroimaging and behavioral data on 900 healthy young adults, aged 22-35. Participants were asked to perform different tasks inside an MRI scanner while blood oxygenation level dependent (BOLD) signals of the brain were recorded. The analysis investigates what brain regions are predictive of the subtle variations of cognitive activity across participants, conditional to other brain regions. Similar to genomics data, the setting is high-dimensional with n = 1556 samples acquired and 156437 brain voxels. A voxel clustering step that reduces data dimension to p = 1000 clusters is done to make the problem tractable.

When decoding brain signals on HCP subjects performing a foot motion experiment (Figure 4, left), AKO recovers an anatomically correct anti-symmetric solution, in the motor cortex and the cerebellum, together with a region in a secondary sensory cortex. KO only detects a subset of those. Moreover, across seven such tasks, the results obtained independently from DL-BH are much more similar to AKO than to KO, as measured with Jaccard index of the resulting maps (Figure 4, right). The maps for the seven tasks are represented in the Appendix. Note that the sign of the effect for significant regions is readily obtained from the regression coefficients, with a voting step for bootstrap-based procedures.

DISCUSSION

In this work, we introduce a p-value to measure knockoff importance and design a knockoffs bootstrapping scheme that leverages this quantity. With this we are able to tame the instability inherent to the original knockoff procedure. Analysis shows that aggregation of multiple knockoffs retains theoretical guarantees for FDR control. However, i) the original argument of [START_REF] Barber | Controlling the false discovery rate via knockoffs[END_REF] no longer holds (see Appendix); ii) a factor κ on the FDR control is lost; this calls for tighter FDR bounds in the future, since we always observe empirically that the FDR is controlled without the factor κ. Moreover, both numerical and realistic experiments show that performing aggregation results in an increase in statistical power and also more consistent results with respect to alternative inference methods.

The quantile aggregation procedure from [START_REF] Meinshausen | p-Values for High-Dimensional Regression[END_REF] used here is actually conservative: as one can see in Figure 2, the control of FDR is actually stricter than without the aggregation step. Nevertheless, as often with aggregation-based approaches, the gain in accuracy brought by the reduction of estimator variance ultimately brings more power.

We would like to address here two potential concerns about FDR control for AKO+BH. The first one is when the {W j } j∈S c are not independent, hence violating Assumption 1. In the absence of a proof of Theorem 1 that would hold under a general dependency, we first note that several schemes for knockoff construction (for instance, the one of [START_REF] Candès | Panning for gold: 'model-X' knockoffs for high dimensional controlled variable selection[END_REF]) imply the independence of (x i -xi ) i∈ [p] , as well as their pseudo inverse. These observations do not establish the independence of W j . Yet, intuitively, the Lasso coefficient of one variable should be much more associated with its knockoff version than with other variables, so it should not be much affected by these other variables, making the Lasso-coefficient differences weakly correlated if not independent. Moreover, in the proof of Lemma 2 and Theorem 1, Assumption 1 is only used for applying Bernstein's inequality, and several dependent versions of Bernstein's inequality have been proved [Sam00, MPR09, HS17, among others]. Similarly, the proof of Eq. ( 16) only uses Assumption 1 for applying the strong law of large numbers, a result which holds true for various kinds of dependent variables (for instance, [START_REF] Abdesselam | The weakly dependent strong law of large numbers revisited[END_REF], and references therein). Therefore we conjecture that independence in Assumption 1 can be relaxed into some mixing condition. Overall, given that the unstability of KO with respect to the KO randomness is an important drawback (see Figure 1), we consider Assumption 1 as a reasonable price price to pay for correcting it, given that we expect to relax it in future works.

The second potential concern is that Theorem 1 is for AKO with k computed from the BY procedure, while BH step-up may not control the FDR when the aggregated p-values (π j ) j∈[p] are not independent. We find empirically that the (π j ) j∈ [p] do not exhibit spurious Spearman correlation (Figure B.6 in Appendix) under a setting where the W j satisfy a mixing condition. This is a mild assumption that should be satisfied, especially when each feature X j only depends on its "neighbors" (as typically observed on neuroimaging and genomics data). It is actually likely that the aggregation step contributes to reducing the statistical dependencies between the (π j ) j∈ [p] . Eventually, it should be noted that BH can be replaced by BY [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF] in case of doubt.

To conclude on these two potential concerns, let us emphasize that the FDR of AKO+BH with B > 1 is always below α (up to error bars) in all the experiments we did, including preliminary experiments not shown in this article, which makes us confident when applying AKO+BH on real data such as the ones of Sections 5.2-5.3.

A practical question of interest is to handle the cases where n p, that is, the number of features overwhelms the number of samples. Note that in our experiments, we had to resort to a clustering scheme of the brain data and to select some genes. A possible extension is to couple this step with the inference framework, in order to take into account that for instance the clustering used is not given but estimated from the data, hence with some level of uncertainty.

The proposed approach introduces two parameters: the number B of bootstrap replications and the γ parameter for quantile aggregation. The choice of B is simply driven by a compromise between accuracy (the larger B, the better) and computation power, but we consider that much of the benefit of AKO is obtained for B ≈ 25. Regarding γ, adaptive solutions have been proposed [START_REF] Meinshausen | p-Values for High-Dimensional Regression[END_REF], but we find that choosing a fixed quantile (0.3) yields a good behavior, with little variance and a good sensitivity. since #{i :

W (i) ≥ W (k) } = k because {W (j) } j∈[p]
is ordered descendingly and because of the assumption that non-zero LCD statistics are distinct. Furthermore, finding the maximum index k of the descending ordered sequence is equivalent to finding the minimum value in that sequence, or

k = min W (k) > 0 : 1 + #{i : W (i) ≤ -W (k) } #{i : W (i) ≥ W (k) } ≤ α ,
since all W (j) ≤ 0 (corresponding with π (k) = 1) have been excluded. Without loss of generality, can write:

t + = min t > 0 : 1 + #{i : W i ≤ -t} #{i : W i ≥ t} ≤ α .
This threshold t + is exactly the same as the definition of threshold τ + in Eq. ( 4) from the original KO procedure.

A.2 Proof of Lemma 1

Setting. Let X ∈ R n×q , β * ∈ R q , λ > 0, σ > 0 be fixed. Define

y = Xβ * + ∀β ∈ R q , L(β) := y -Xβ 2 2 + λ β 1 .
with ∼ N (0, σI n ) the Gaussian noise and • p the L p norm.

Classical Optimization Properties.

Since L is convex, non-negative, and tends to +∞ at infinity, its minimum over R q exists and is attained (although may not be unique). Since L is convex, its minima are characterized by a first-order condition:

β λ ∈ argmin β∈R q {L(β)} ⇔ 0 ∈ ∂L(β)
which is equivalent to

     ∃ z ∈ [-1, 1] q : X X β λ = X y - λ 2 z ∀j s.t. ( β λ ) j = 0, z j = sign(( β λ ) j ) (A.17)
As shown by [Gir14, Section 4.5.1] for instance, the fitted value f λ ∈ R n is uniquely defined:

∃! f λ ∈ R n such that ∀ β λ ∈ argmin β∈R q {L(β)} , f λ = X β λ .
As a consequence, the equicorrelation set

J λ := j ∈ {1, . . . , q} : x j (y -X β λ ) = λ/2
is uniquely defined. We also have,

(A.18) ∀ β λ ∈ argmin β∈R q {L(β)} , j : ( β λ ) j = 0 ⊂ J λ
(but these two sets are not necessarily equal, and the former set may not be uniquely defined).

Note that for every set J ⊂ {1, . . . , q} such that ∀j / ∈ J, ( β λ ) j = 0, we have X β λ = X J ( β λ ) J so that (X X β λ ) J = X J X J ( β λ ) J . As a consequence, taking J = J λ , by eq. (A.17) and (A.18), any minimizer

β λ of L over R q satisfies (A.19) X J λ X J λ ( β λ ) J λ = X J λ y - λ 2 z J λ for some z J λ ∈ {-1, 1} J λ . If the matrix X J λ X J λ is non-singular (that is, if X J λ is of rank | J λ |), then the argmin of L is unique [Gir14, Section 4.5.1].
Result A.1. For every α ∈ R q \{0}, the event rank(X J λ ) = | J λ | , α β λ = 0 and ∃ j ∈ {1, . . . , q} , α j ( β λ ) j = 0 , (A.20) where { β λ } = argmin β∈R q {L(β)} is well-defined by the first property, has probability zero.

Proof. Let Ω be the event defined by Eq. (A.20). If Ω holds true, then there exists some J ⊂ {1, . . . , q} and some z ∈ {-1, 1} q such that rank(X J ) = |J|, ( β λ ) J c = 0, and

X J X J ( β λ ) J = X J y - λ 2 z J .
Indeed, this is a consequence of Eq. (A.18) and (A.19), by taking J = J λ and z such that z J λ = sign ( β λ ) J λ . Therefore, using that X J X J is non-singular, we get

( β λ ) J = M (J) + v(J, z)
where

M (J) := (X J X J ) -1 X J and v(J, z) := (X J X J ) -1 X J Xβ * -(X J X J ) -1 λ 2 z J , hence α β λ = α J M (J) + α v(J, z)
follows a normal distribution with variance σ 2 α J M (J)M (J)

α J = σ 2 M (J) α J 2 .
Now, on Ω, we also have the existence of some j such that α j ( β λ ) j = 0. Since ( β λ ) J c = 0, we must have j ∈ J, which shows that α J = 0.

Overall, we have proved that

Ω ⊂ J∈J ,z∈{-1,1} q Ω J,z
where J := {j ∈ {1, . . . , q} : rank(X J ) = |J| and α J = 0} and Ω J,z := α J M (J) + α v(J, z) = 0 .

For every J ∈ J , M (J) α J = 0 since α J = 0 and M (J) is of rank |J|. As a consequence, for every J ∈ J and z ∈ {-1, 1} p , P(Ω J,z ) is the probability that a Gaussian variable with non-zero variance is equal to zero, so it is equal to zero. We deduce that

P(Ω) ≤ J∈J ,z∈{-1,1} q P(Ω J,z ) = 0
since the sets J and {-1, 1} q are finite.

Applying Result A.1 to the case where X concatenates the original p covariates and their knockoff counterparts (hence q = 2p), we get that, apart from the event where X J λ is not full rank, for every j ∈ {1, . . . , p}, W j takes any fixed non-zero value with probability zero (with α j = ±1, α j+p = ±1, α k = 0 otherwise).

Similarly, the above lemma shows that for every j = j ∈ {1, . . . , p}: P(X J λ is full-rank and ∃j = j , W j = W j , W j = 0, W j = 0) = 0.

As a consequence, with probability 1, all the non-zero W j are distinct if X J λ is fullrank.

Remark A.4. The proof of Result A.1 is also valid for other noise distributions: it only assumes that the support of the distribution of is not included into any hyperplane of R n .

A.3 Asymptotic Validity of Intermediate P-values

We consider in this section an asymptotic regime where p → +∞.

Assumption A.2 (Asymptotic regime p → +∞). When p grows to infinity, n, X, β * , and y all depend on p implicitly. We assume that for every integer j ≥ 1, 1 β * j =0 does not depend on p (as soon as p ≥ j), and that

|S| p = {j ∈ [p] : β * j = 0} p ----→ p→+∞ 0 .
When making Assumption 1, we also assume that P 0 does not depend on p.

Lemma A.3. If Assumptions 1 and A.2 hold true, then for all j ≥ 1 such that β * j = 0, the empirical p-value π j defined by Eq. ( 5) is a valid p-value asymptotically, that is,

∀t ∈ [0, 1] , lim p→+∞ P(π j ≤ t) ≤ t .
Note that our proof of Theorem 1 in Section A.5 relies on the use of Lemma 2 with t that can be of order 1/p. Therefore, Lemma A.3 above is not sufficient for our needs. Nevertheless, it still provides a interesting insight about the π j , and justifies (asymptotically) their name, which is why we state and prove this result here.

Proof. By definition, π j ≤ 1 almost surely, so the result holds when t = 1. Let us now focus on the case where t ∈ [0, 1). Let F 0 denote the c.d.f. of P 0 , the common distribution of the null statistics {W j } 1≤j≤p / β * j =0 , which exists by Assumption 1. Let j ≥ 1 such that β * j = 0 be fixed, and assume that p ≥ j is large enough so that |S c | ≥ 2. Let m = |S c | -1 ≥ 1 as in the proof of Lemma 2. Note that m depends on p, and m/p → 1 as p → +∞ by Assumption A.2, hence m → +∞ as p → +∞.

By definition of π j , when W j > 0 we have: ----→ p→+∞ F 0 (u) .

π j = 1 + #{k ∈ [p] : W k ≤ -W j } p (since W j > 0 > -W j ) = 1 + #{k ∈ S : W k ≤ -W j } + #{k ∈ S c \ {j} : W k ≤ -W j } p ≥ #{k ∈ S c \ {j} : W k ≤ -W j } p = F m (-W j ) α p (A.
Since we assume lim p→+∞ |S| /p = 0, lim p→+∞ α p = 1 and we get that for all u ∈ R,

1 α p F m (u)
a.s.

----→ p→+∞ F 0 (u) .

Since W j is independent from {W k } k∈S c \{j} , this result also holds true conditionally to W j , with u = -W j . Given that almost sure convergence implies convergence in distribution, we have: conditionally to W j ,

(A.22) 1 α p F m (-W j ) (d) ----→ p→+∞ F 0 (-W j ) (d) = F 0 (W j )
where the latter equality comes from the fact that W j has a symmetric distribution, as shown in Remark 1.

So, when W j > 0, for every t ∈ [0, 1), Therefore, for every t ∈ [0, 1), lim sup p→+∞ P(α p π j ≤ t) = lim sup p→+∞ P(α p π j ≤ t and W j > 0) + P(α p π j ≤ t and W j ≤ 0)

=0 since αp≥1>t and π j =1 when

W j ≤0 = lim sup p→+∞ E[P(α p π j ≤ t | W j )1 W j >0 ] ≤ E lim sup |S c |→+∞ {P(α p π j ≤ t | W j )1 W j >0 )} ≤ P F 0 (W j ) ≤ t by Eq. (A.23) ≤ t ,
which concludes the proof.

A.4 A General FDR Control with Quantile-aggregated P-values

The proof of Theorem 1 relies on an adaptation of results proved by [MMB09, Theorems 3.1 and 3.3] about aggregation of p-values. The results of [START_REF] Meinshausen | p-Values for High-Dimensional Regression[END_REF], whose proof relies on the proofs of [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF], are stated for randomized p-values obtained through sample splitting. The following lemma shows that they actually apply to any family of p-values. : 1 ≤ b ≤ B} where q γ (•) is the γ-quantile function,

h ∆ = max i ∈ [p] : Q (i) ≤ iᾱ where Q (1) ≤ • • • ≤ Q (p) ,
and

S ∆ = j ∈ [p] : Q j ≤ Q ( h) .
Then,

E | S ∩ N | | S| ∨ 1 ≤ p-1 j=1 1 j - 1 j + 1 F (j) + F (p) p (A.24) where ∀j ∈ [p] , F (j)
∆ = 1 γ 1 B B b=1 i∈N P π (b) i ≤ j ᾱγ p .
As a consequence, if some C ≥ 0 exists such that

(A.25) ∀t ≥ 0, ∀b ∈ [B] , ∀i ∈ N , P π (b) i ≤ t ≤ Ct , then we have (A.26) E | S ∩ N | | S| ∨ 1 ≤ |N |C p   p j=1 1 j   ᾱ .
Let us emphasize that Lemma A.4 can be useful in general, well beyond knockoff aggregation. To the best of our knowledge, Lemma A.4 is new. In particular, the recent preprint by [START_REF] Romano | Multiple Data Splitting for Testing[END_REF], that studies p-values aggregation procedures, focuses on FWER controlling procedures, whereas Lemma A.4 provides an FDR control for a less conservative procedure.

Proof. For every i, j, k ∈ [p], let us define

p i,j,k =    P Q i ∈ (j -1)ᾱ, j ᾱ , i ∈ S and | S| = k if j ≥ 2 P Q i ∈ [0, ᾱ] , i ∈ S and | S| = k if j = 1 . Then, | S ∩ N | | S| ∨ 1 = p k=1 1 | S|=k i∈N 1 i∈ S k = i∈N p k=1 1 k 1 | S|=k and i∈ S = i∈N p k=1 1 k 1 | S|=k and i∈ S and 0≤Q i ≤k ᾱ since i ∈ S and | S| = k implies that Q i ≤ Q ( h) ≤ hᾱ = k ᾱ.
Taking an expectation and writing that

1 0≤Q i ≤k ᾱ = 1 Q i ∈[0, ᾱ] + k j=2 1 Q i ∈((j-1) ᾱ,j ᾱ] ,
we get -following the computations of [MMB09, proof of Theorems 3.3], which themselves rely on the ones of [BY01]-, Since the above upper bound is equal to

E | S ∩ N | | S| ∨ 1 ≤ i∈N p k=1 1 k k j=1 p i,j,k = i∈N p j=1 p k=j 1 k p i,j,k
F (1) + p j=2 1 j F (j) -F (j -1) = p j=1 1 j - 1 j + 1 F (j) + F (p) p , we get that E | S ∩ N | | S| ∨ 1 ≤ p j=1 1 j - 1 j + 1 F (j) + F (p) p . (A.27) Notice also that F (j) = i∈N P(Q i ≤ j ᾱ and i ∈ S) ≤ i∈N P(Q i ≤ j ᾱ) .

Now, as done by [MMB09, proof of Theorems 3.1], we remark that

Q i ≤ j ᾱ is equivalent to 1 B b ∈ [B] : pπ (b) i γ ≤ j ᾱ = 1 B B b=1 1 pπ (b) i ≤j ᾱγ ≥ γ so that P(Q i ≤ j ᾱ) = P 1 B B b=1 1 pπ (b) i ≤j ᾱγ ≥ γ ≤ 1 γ E 1 B B b=1 1 pπ (b) i ≤j ᾱγ by Markov inequality = 1 γ 1 B B b=1 P pπ (b) i ≤ j ᾱγ .
Therefore,

F (j) ≤ i∈N 1 γ 1 B B b=1 P pπ (b) i ≤ j ᾱγ = F (j) ,
so that Eq. (A.27) implies Eq. (A.24).

If condition (A.25) holds true, then, for every j ∈ [p],

F (j) ≤ |N |C γ j ᾱγ p = |N |C ᾱ p j , hence Eq. (A.24) shows that E | S ∩ N | | S| ∨ 1 ≤ p-1 j=1 F (j) j(j + 1) + F (p) p ≤ |N |C ᾱ p p-1 j=1 1 j + 1 + |N |C ᾱ p = |N |C ᾱ p   p j=1 1 j   ,
which is the desired result.

A.5 Proof of Theorem 1

We can now prove Theorem 1. We apply Lemma A.4 with ᾱ = β(p)α, N = S c , so that S = S AKO+BY . Since the π

APPENDIX B: ADDITIONAL EXPERIMENTAL RESULTS

B.1 Demonstration of Aggregated Multiple Knockoff vs. Simultaneous Knockoff

Using the same simulation settings as in Section 5.1 with n = 500, p = 1000 and varying simulation parameters to generate Figure 2 in the main text, we benchmark only aggregation of multiple knockoffs (AKO) with 5 and 10 bootstraps (B = 5 and B = 10) and compare with simultaneous knockoffs with 2 and 3 bootstraps (κ = 2 and κ = 3). Results in Figure B.5 show that while increasing the number of knockoff bootstraps makes simultaneous knockoffs more conservative, doing so with AKO makes the algorithm more powerful (and in the worst case retains the same power with fewer bootstraps).

B.2 Empirical Evidence on the Independence of Aggregated P-values π

Using the same simulation settings as in Section 5.1 with n = 500, p = 1000, ρ = 0.6, snr = 3.0, sparsity = 0.06 we generate 100 observations of aggregated p-values π. Then, we compute the Spearman rank-order correlation coefficient of the Null πj for these 100 observations along with their two-sided p-value (for a hypothesis test whose null hypothesis is that two sets of data are uncorrelated).

The results are illustrated in Figure B.6: the Spearman correlation values are concentrated around zero, while the distribution of associated p-values seems to be a mixture between a uniform distribution and a small mixture component consisting of mostly non-null p-values. This indicates near independence between the aggregated p-values using quantile-aggregation [START_REF] Meinshausen | p-Values for High-Dimensional Regression[END_REF], hence justifies our use of BH step-up procedure for selecting FDR controlling threshold in the AKO algorithm.

Remark B.5. Again, it is worth noticing that the empirical evidence we have shown is only done in a setting with a Toeplitz structure for the covariance matrix. However, as explained in the main text, this correlation setting is usually found in neuroimaging and genomics data. Hence, we believe that assuming short-distance correlations between the (X j ) j∈[p] is a mild assumption, which should be satisfied in the practical scenarios where we want to apply our algorithm.

The decoding maps returned by the KO, AKO and DL inference procedures are presented in Figure B.7. As quantified by the Jaccard index in the main text, we observe that the AKO solution is typically closer to an external method based on the desparsified lasso (DL). Moreover, AKO is also typically more sensitive than KO alone.

The seven classification problems are the following:

• Emotion: predict whether the participant watches an angry face or a geometric shape. • Gambling: predict whether the participant gains or loses gambles.

• Motor foot: predict whether the participant moves the left or right foot.

• Motor hand: predict whether the participant moves the left or right hand.

• Relational: predict whether the participant matches figures or identified feature similarities. • Social: predict whether the participant watches a movie with social behavior or not. • Working Memory: predict whether the participant does a 0-back or a 2-back task. Emotion, gambling, motor foot, motor hand, relational, social and working memory refer to 7 binary tasks that were considered based on the HCP900 dataset. We observe that AKO is typically more sensitive than KO, and yields solution closer to that of an independent solution based on a desparsified-Lasso (DL) estimator.

k

  ← fdr_threshold(α, (π1, π2, . . . , πp)) // Using either Eq. (7) or Eq. (8)

Figure 1 :

 1 Figure 1: Histogram of FDP and power for 2500 runs of KO (blue, top row) vs. 100 runs of AKO with B = 25 (teal, bottom row) under the same simulation. Simulation parameter: SNR = 3.0, ρ = 0.5, sparsity = 0.06. FDR is controlled at level α = 0.05.

Figure 2 :

 2 Figure2: FDR (left) and average power (right) of several methods for 100 runs with varying simulation parameters. For each varying parameter, we keep the other ones at default value: SNR = 3.0, ρ = 0.5, sparsity = 0.06. FDR is controlled at level α = 0.1. The benchmarked methods are: aggregation of multiple knockoffs (AKO -ours); vanilla knockoff (KO); simultaneous knockoff by[START_REF] Gimenez | Improving the Stability of the Knockoff Procedure: Multiple Simultaneous Knockoffs and Entropy Maximization[END_REF] (KO-GZ) and by[START_REF] Emery | Controlling the FDR in variable selection via multiple knockoffs[END_REF] (KO-EK); knockoff statistics aggregation (KO-HH); debiased-Lasso (DL-BH).

Figure 3 :

 3 Figure 3: FDR and average power for 30 simulations of fixed parameters: SNR=3.0, ρ = 0.7, sparsity=0.06. There is virtually no gain in statistical power when B > 25 and when γ ≥ 0.1.

Figure 4 :

 4 Figure 4: Detection of significant brain regions for HCP data (900 subjects). (left) Selected regions in a left or right foot movement task. Orange: brain areas with positive sign activation. Blue: brain areas with negative sign activation. Here the AKO solution recovers an anatomically correct pattern, part of which is missed by KO. (right) Jaccard index measuring the Jaccard similarity between the KO/AKO solutions on the one hand, and the DL solution on the other hand, over 7 tasks: AKO is significantly more consistent with the DL-BH solution than KO.

  ∈ S c \ {j} : W k ≤ u} m is the empirical cdf of {W k } k∈S c \{j} .Now, since {W k } k∈S c \{j} are i.i.d. with distribution P 0 by Assumption 1 , the law of large numbers implies that, for all u ∈ R, F m (u) a.s.

  j ≤ t | W j ) ≤ lim sup p→+∞ P F m (-W j ) α p ≤ tW j by Eq. (A.21) ≤ 1 F 0 (W j )≤t (A.23) by Eq. (A.22) combined with the Portmanteau theorem.

Lemma A. 4 .

 4 Let (π (b) j ) 1≤j≤p,1≤b≤B be a family of random variables with values in [0, 1]. Let γ ∈ (0, 1], ᾱ ∈ [0, 1] and N ⊂ [p] be fixed. Let us define ∀j ∈ [p] , Q j

j

  , b = 1, . . . , B, have the same distribution as π j , by Lemma 2, condition (A.25) holds true with C = κp/|S c |, and Eq. (A.26) yields the desired result. Note that an FDR control for AKO such as Theorem 1 cannot be obtained straightforwardly from the arguments of [BC15] and [CFJL18]. One key reason for this is that their proof relies on a reordering of the features according to the values of (|W j |) j∈[p] [BC15, Section 5.2], such a reordering being permitted since the signs of the W j are iid coin flips conditionally to the (|W j |) j∈[p] [CFJL18, Lemma 2]. In the case of AKO, we must handle the (W (b) j ) j∈[p] simultaneously for all b ∈ [B], and conditioning with respect to (|W (b) j |) j∈[p],b∈[B] may reveal some information about the signs of the (W (b) j ) j∈[p] as soon as B > 1. At least, it does not seem obvious to us that the key result of [CFJL18, Lemma 2] can be proved conditionally to the (|W (b) j |) j∈[p],b∈[B] when B > 1, so that the proof strategy of [CFJL18] breaks down in the case of AKO with B > 1.

Figure B. 5 :Figure B. 6 :Figure B. 7 :

 567 Figure B.5: Aggregation of multiple knockoffs (B = 5 and B = 10) vs. simultaneous knockoffs (κ = 2 and κ = 3).A clear loss in statistical power is demonstrated in the latter method when increasing the number of bootstraps κ, while the former (AKO) shows the opposite: with B = 10 bootstraps there are small, yet consistent gains in the number of true detections compared to using only B = 5 bootstraps.

Table 1 List of detected genes associated with phenotype FT_GH.

 1 Empty line (-) signifies no detection. Detected genes are listed in well-known studies dated up to 20 years ago.

	Method Detected Genes
	AKO	AT2G21070, AT4G02780, AT5G47640
	KO	AT2G21070
	KO-GZ	AT2G21070
	DL-BH	-

ACKNOWLEDGEMENTS

The authors thank anonymous reviewers for their helpful comments and suggestions. We are grateful for the help of Lotfi Slim and Chloé-Agathe Azencott on the Arabidopsis thaliana dataset, and the people of Human Connectome Project on HCP900 dataset. This research is supported under funding of French ANR project FastBig (ANR-17-CE23-0011), the KARAIB AI chair and Labex DigiCosme (ANR-11-LABEX-0045-DIGICOSME).

APPENDIX

The Appendix is organized as follows. First, the main theoretical results of the article are proved:

• Proof of Proposition 1: AKO+BH with B = 1 and γ = 1 is equivalent to vanilla KO. • Proof of Lemma 1: for Lasso-coefficient differences, the non-zero W j are distinct. Second, the results of some additional experiments are reported:

• Additional experiments to show that the KO-GZ alternative aggregation procedure by [START_REF] Gimenez | Improving the Stability of the Knockoff Procedure: Multiple Simultaneous Knockoffs and Entropy Maximization[END_REF] has decreasing power when the number κ of knockoff vectors x considered simultaneously increases (we compare κ = 2 with κ = 3). We show empirically that this is not the case for AKO with respect to B.

• Empirical evidence for the near independence of p-values π j .

• Additional figures for HCP 900 experiments.

APPENDIX A: DETAILED PROOFS

A.1 Proof of Proposition 1

We begin by noticing that the function f : R

This means the first step of both FDR control step-up procedures, that involves ordering the intermediate p-values ascendingly, is the same as arranging the knockoff statistic in descending order:

. Therefore from Eq. ( 7) and the definition of π j we have:

(note that we can exclude all the π (k) = 1 due to the fact that ∀ k ∈ [p], α ∈ (0, 1) : kα/p < 1).

This can be written as: