Marco Gramaglia

Oscar Trullols-Cruces

Diala Naboulsi

Marco Fiore

Maria Calderon

Vehicular Networks on Two Madrid Highways

Keywords: vehicle-to-vehicle communications, vehicular networks, mobility modeling, connectivity analysis communications entre vehicules, réseaux vehiculaires, modélisation de mobilité, analyse de connectivité

There is a growing need for vehicular mobility datasets that can be employed in the simulative evaluation of protocols and architectures designed for upcoming vehicular networks. Such datasets should be realistic, publicly available, and heterogeneous, i.e., they should capture varied traffic conditions. In this paper, we contribute to the ongoing effort to define such mobility scenarios by introducing a novel set of traces for vehicular network simulation. Our traces are derived from high-resolution real-world traffic counts, and describe the road traffic on two highways around Madrid, Spain, at several hours of different working days. We provide a thorough discussion of the real-world data underlying our study, and of the synthetic trace generation process. Finally, we assess the potential impact of our dataset on networking studies, by characterizing the connectivity of vehicular networks built on the different traces. Our results underscore the dramatic impact that relatively small communication range variations have on the network. Also, they unveil previously unknown temporal dynamics of the topology of highway vehicular networks, and identify their causes.

Introduction

Vehicular communications are regarded as a key enabling technology within upcoming Intelligent Transportation Systems (ITS), which are in turn expected to significantly improve road safety and traffic management. Indeed, communication-enabled vehicles will be able to exchange data in a direct way, via the so-called vehicle-to-vehicle (V2V) communication, so as to build spontaneous, distributed, self-organizing vehicular networks. The latter are anticipated to complement the existing cellular network, by providing short-range, low-latency, inexpensive transfers that well integrate today's pervasive but strongly centralized radio access network architecture. The same features also make vehicular networks especially suitable to implement ITS services such as cooperative awareness, collision avoidance, and danger warning.

The moment when V2V technologies and the associated services will hit the market is in fact approaching. Standards especially designed for V2V communication, such as IEEE 802.11p, IEEE 1609 OSI CALM-M5 and ETSI ITS are being finalized, and the automotive industry is actively participating in the first field tests, such as those carried out via the sim TD platform in Frankfurt, Germany, or through the Ann Arbor Safety Pilot in Michigan, USA.

However, experimental trials remain an exception, due to their costs and complexity. The vast majority of applications, protocols and architectures for upcoming vehicular networks is still evaluated via simulative studies. Within such a context, the level of realism of the simulation is a paramount aspect to account for, and the way the mobility of individual vehicles is represented is often the single feature that introduces the largest bias in the results [START_REF] Fiore | The networking shape of vehicular mobility[END_REF].

During the past decade, there have been efforts aiming at gathering real-world road traffic data [START_REF] Huang | Performance evaluation of suvnet with real-time traffic data[END_REF][3][4], developing tools for the simulation of vehicular movement [5,[START_REF] Härri | Vehicular mobility simulation with VanetMo-biSim[END_REF], and generating synthetic mobility datasets [START_REF] Raney | An agent-based microsimulation model of swiss travel: first results[END_REF][START_REF] Uppoor | Generation and analysis of a large-scale urban vehicular mobility dataset[END_REF][START_REF] Pigné | A vehicular mobility model based on real traffic counting data[END_REF]. This notwithstanding, there is still a noticeable scarcity of road traffic datasets featuring the level of realism and of spatiotemporal granularity required for network simulation [START_REF] Uppoor | Synthetic mobility traces for vehicular networking[END_REF]. The result of the lack of a set of realistic, publicly shared, heterogeneous scenarios is that simulations of vehicular networks are often unreliable and non-reproducible [START_REF] Joerer | Toward reproducibility and comparability of IVC simulation studies: a literature survey[END_REF].

In this paper, we contribute to the endeavor of enlarging the set of realistic datasets of vehicular mobility that are freely available to the research community. To that end, we present 16 original traces describing the road traffic on two highways around Madrid, Spain, at different times of multiple weekdays. The traces are derived from high-detail real-world traffic counts, and describe unidirectional free-flow traffic in quasi-stationary conditions on a road spanning over 10 km. More precisely, our contributions are as follows:

• We collect real-world traffic count data that feature an unprecedented resolution, and study the road traffic dynamics they capture. Our analysis pinpoints the heterogeneity of the empirical dataset, and confirms that widely adopted assumptions on exponential inter-arrivals of vehicles do not necessarily hold in real scenarios.

• We present a process to generate synthetic mobility traces from road traffic counts, by leveraging microscopic models and properties of the real-world data. Our analysis is relevant to the calibration of car-following and lane-changing models, and results in original, publicly available datasets of highway mobility for network simulation.

• We provide a characterization of the vehicular network connectivity in the synthetic traces, investigating the impact of time (i.e., hour of the day, day of the week), highway diversity, and vehicular communication range. Our results highlight that: (i) the level of global vehicular connectivity can be primarily ascribed to the V2V communication range; (ii) for short communication ranges, the vehicular connectivity is driven by a mixture of slow and fast temporal dynamics that are regulated by the road traffic density and by V2V contact durations, respectively.

The paper is organized as follows. After a discussion of related works, in Sec. 2, we present and analyze the real-world traffic count data we use in our study, in Sec. 3. Then, Sec. 4 describes the generation process of the synthetic vehicular mobility traces, whose vehicular network connectivity is characterized in Sec. 5. Finally, Sec. 6 concludes the paper.

Related Work

The impact of realistic mobility models on the simulation of communication protocols for vehicular networks has been emphasized in many works [START_REF] Fiore | The networking shape of vehicular mobility[END_REF][START_REF] Uppoor | Generation and analysis of a large-scale urban vehicular mobility dataset[END_REF][START_REF] Joerer | Toward reproducibility and comparability of IVC simulation studies: a literature survey[END_REF][START_REF] Bai | Spatio-temporal variations of vehicle traffic in VANETs: facts and implications[END_REF]. This fact has pushed the research community to seek for an ever-increasing realism in road traffic traces used to feed network simulators.

A first approach is to use real-world mobility traces that have been recorded by directly logging the position of vehicles during their movements. Unfortunately, the traces of this kind that are currently available yield a partial description of the road traffic, as they only report the positions of, e.g., buses [START_REF] Huang | Performance evaluation of suvnet with real-time traffic data[END_REF]3] or taxis [4]. Also, they are typically affected by limited temporal accuracy, with positions logged with periodicities in the order of tens of seconds [START_REF] Uppoor | Synthetic mobility traces for vehicular networking[END_REF].

Therefore, most works have focused instead on the generation of synthetic vehicular traces. In particular, a number of studies have considered urban mobility scenarios, created by feeding real-world road topologies of different cities to microscopic traffic simulators such as SUMO [5] or VanetMobiSim [START_REF] Härri | Vehicular mobility simulation with VanetMo-biSim[END_REF]. In order to characterize the number, origin, destination and time of trips, these works usually made use of macroscopic data (i.e., origin-destination matrices) collected from user surveys [START_REF] Raney | An agent-based microsimulation model of swiss travel: first results[END_REF][START_REF] Uppoor | Generation and analysis of a large-scale urban vehicular mobility dataset[END_REF][START_REF] Caliskan | Collection of Dedicated Information in Vehicular Ad Hoc Networks[END_REF] or from roadside detectors such as induction loops, cameras and infrared counters [START_REF] Pigné | A vehicular mobility model based on real traffic counting data[END_REF]. However, the dynamics of traffic in urban regions -characterized by vehicles traveling at low or medium speed, and frequent intersections regulated by traffic lights or roundabouts -are not comparable to those on highways -featuring instead high speeds and frequent overtakings. Moreover, none of the aforementioned works considers the problem of fine-tuning the microscopic mobility model, as we do in this study.

Closer to our approach, the study in [START_REF] Bai | Spatio-temporal variations of vehicle traffic in VANETs: facts and implications[END_REF] uses two empirical road traffic datasets to generate synthetic highway mobility traces that can be used in network simulation. The two 24-hour datasets were collected using detectors underneath the road: the first dataset is from the I-80 highway nearby Berkeley, CA, USA, and contains information from dual-loop detectors (that log data on individual vehicles passing through); the second dataset is from the Gardiner Expressway in Toronto, Canada, and includes records from metal detectors (that measure aggregated characteristics on each lane every 20 seconds). In order to generate the synthetic traces, inter-arrivals are assumed to follow exponential distributions and vehicle speeds are modeled as Gaussian distributions: the real-world data is then employed to calculate the parameters of such distributions for different hours of a day. Then, a stochastic vehicle mobility generator is developed according these models. Our study improves that in [START_REF] Bai | Spatio-temporal variations of vehicle traffic in VANETs: facts and implications[END_REF] from several viewpoints: (i) we consider more detailed and heterogeneous traffic count data; (ii) we use the empirical traffic counts as they are, and not to parameterize distributions; (iii) we employ validated microscopic car-following and lane-changing models instead of resorting to stochastic modeling; (iv) the synthetic mobility traces we generate are publicly available, which is not the case for those of I-80 and Gardiner Expressway.

Highway scenarios are considered also in [START_REF] Akhtar | VANET topology characteristics under realistic mobility and channel models[END_REF], where empirical data from the Freeway Performance Measurement System (PeMS)1 is fed to the SUMO microscopic mobility simulator to generate synthetic traces. In particular, the real-world data, collected from road sensors in Alameda County, CA, USA, is used to determine the assignment of vehicular traffic flows and the average speed values over the road. However, the empirical traffic count dataset used in [START_REF] Akhtar | VANET topology characteristics under realistic mobility and channel models[END_REF] is affected by a coarse time granularity (i.e., flow and speed sampling rates between 30 seconds and 5 minutes), which prevents using it to model the precise arrival time of each vehicle, as we do instead. Also relevant to our work are vehicular network connectivity studies. There, several research teams have focused on urban areas [4,[START_REF] Naboulsi | On the instantaneous topology of a large-scale urban vehicular network: the Cologne case[END_REF][START_REF] Thakur | Modeling and characterization of urban vehicular mobility using web cameras[END_REF][START_REF] Viriyasitavat | Dynamics of network connectivity in urban vehicular networks[END_REF]. Although interesting, their findings are hardly applicable to the highway scenarios we are interested in, due again to their significant differences with respect to urban environments. As far as highway traffic is concerned, a number of works dealt with the analytical characterization of vehicular connectivity in unidimensional scenarios [START_REF] Yousefi | Analytical model for connectivity in vehicular ad hoc networks[END_REF][START_REF] Mohimani | Mobility Modeling, Spatial Traffic Distribution, and Probability of Connectivity for Sparse and Dense Vehicular Ad Hoc Networks[END_REF][START_REF] Gramaglia | Overhearing-Assisted Optimization of Address Autoconfiguration in Position-Aware VANETs[END_REF], which however implies adopting simplifying assumptions in order to make the problem mathematically tractable.

Fewer papers rely instead on realistic mobility traces. Specifically, in [START_REF] Artimy | Connectivity in inter-vehicle ad hoc networks[END_REF] synthetic data generated via a microscopic simulator is used to confirm the intuition that increasing the traffic density or the number of lanes improves connectivity. Nevertheless, the study employs a simple highway traffic model with deterministic vehicle inter-arrivals, while our analysis relies on much more realistic mobility traces built from real-world data. More realistic mobility is accounted for in [START_REF] Füssler | Studying vehicle movements on highways and their impact on ad-hoc connectivity[END_REF], which uses (publicly unavailable) reality-audited highway movements from the FleetNet project, and in [START_REF] Akhtar | VANET topology characteristics under realistic mobility and channel models[END_REF], which employs the aforementioned Alameda County trace. In both cases, however, the connectivity assessment is limited to statistical analyses of node degree and link duration distributions. In this paper, we provide a more thorough characterization of the instantaneous vehicular network connectivity, including original quantitative results on the correlation between the radio communication range, the traffic density, and the number and size of network components.

Real-world traffic count data

The empirical data used throughout this work comes from measurement locations in the region of Madrid, Spain. The data, kindly provided to us by the Spanish office for the traffic management (Dirección General de Tráfico, DGT), detail road traffic conditions on the following two arterial highways.

• M40: the M40 is part of the middle layer of the Madrid city beltway system, which also comprises M-30 (innermost) and M-50 (outermost). This beltway, which has an average distance of 10.7 km from the city center, traverses both the most peripheral areas of the municipality and surrounding cities. The measurement point is placed at the 12.7-km milepost, where the M40 traverses the suburb of San Blas and the town of Coslada. The measures cover the internal carriageway (southbound) that includes 3 lanes with a speed limit of 100 km/h. The M40 is indicated as A in Fig. 1(a), while a photographic view of the exact measurement location is provided in Fig. 1(b).

• A6: the Autovía A6 is a motorway that connects the city of A Coruña to the city of Madrid. This road enters into the urban area from the northwest and collects the traffic demand of the conurbation that was built along it. The data collection point is placed around the 11-km milepost (Madrid direction), where the A6 features three lanes with a speed limit of 120 km/h. We remark, however, that this location is right after a popular entrance ramp that joins the rightmost lane and significantly slows down the road traffic.

The highway maps to B in Fig. 1(a). A photographic view of the collection point is in Fig. 1(c).

Collecting high-resolution traffic counts

The sensors that DGT deployed on M40 and A6 are induction loops. They consist of loops of wires buried under the concrete layer that create a magnetic field. When a vehicle passes on the vertical axis of the loop, it creates a variation in the magnetic field, which is considered as a new transit. If two loops are placed close to each other, also the vehicular speed (and possibly other metrics, such as the vehicle length) is obtained. A picture of the buried induction loops used to collect measurements on the M40 is provided in Fig. 1(d).

Usually, these devices are programmed to supply coarse-grained data, as public transportation authorities are generally interested in aggregate measures2 (e.g., the number of vehicles transiting on a road, their average speed, or the percentage of heavy vehicles) so as to detect exceptional alterations in the road capacity [START_REF] Bai | Spatio-temporal variations of vehicle traffic in VANETs: facts and implications[END_REF][START_REF] Akhtar | VANET topology characteristics under realistic mobility and channel models[END_REF]. The loops installed by DGT were configured to supply data averaged over 60-second intervals, but their configuration was changed specifically for our study: the real-world traffic counts we collected provide fine-grained time, speed and lane information on each single transiting vehicle. To the best of our knowledge, this is the most precise traffic count dataset recorded to date, and represents an ideal input to our microscopic simulation of highway traffic.

The timing of the data collection is a very important variable to account for, since vehicular traffic presents, as many other human activities, important temporal variability: for instance, rush hours yield diverse traffic conditions than off-peak hours, especially on main arterial roads like those we consider. The dataset analyzed in this work was collected on multiple days of May 2010, namely Friday the 7th, Monday the 10th, Tuesday the 11th, and Wednesday the 12th. For each measurement location and each day, two traces are available: one capturing the traffic peak (from 8:30 a.m. to 9:00 a.m.), and one during off-peak hours (from 11:30 a.m. to 12:00 p.m.). As a result, our empirical dataset provides a heterogeneous view of traffic conditions, allowing us to generate synthetic mobility traces that are representative of different highway congestion levels.

Understanding the data

Each entry of the dataset refers to one individual vehicle transiting over the induction loops, and includes the following:

• Timestamp: the time at which the vehicle passage was recorded by the induction loop.

The precision of the time reference is 100 milliseconds. • Speed: the vehicle speed, in km/h.

• Lane: the lane on which the vehicle transited. Both M40 and A6 feature three lanes at the measurement points. They are referred to as the right, center and left lanes (in the heading direction) in the remainder of the paper.

In [START_REF] Gramaglia | New insights from the analysis of free flow vehicular traffic in highways[END_REF], we performed an extensive analysis of vehicle inter-arrival times (IAT) in the traffic count dataset. We found that, although many works assume exponential distribution or even uniform distribution to model the time headway between subsequent vehicles, a mixture gaussianexponential model yields a much better approximation of the empirical data. The mixture model also provides valuable information on drivers' behavior: the gaussian part of the distribution captures bursty arrivals of vehicles that travel close to each other at similar speeds, whereas the exponential part models isolated vehicles whose movement is less constrained by that of other cars. Fig. 2 shows the excellent match between the mixture model and the experimental data on multiple combinations of highway, day and hour for each lane. An overview of two traffic count data samples is provided in Fig. 3. Each plot depicts the evolution over time of the speed on the three lanes of the highway under consideration: consistently with Spanish road regulations, the speed of the right lane is typically the lowest, while that of the left lane is normally the highest. As anticipated, the dataset used in this work presents different traffic characteristics in terms of traffic measures (i.e., density, flow and speed). However, we remark that, throughout all the dataset, the road traffic is always in a so-called free flow state, i.e., there is a quasi-linear relationship between the vehicular flow and density, and no significant congestion is observed [START_REF] Kerner | The physics of traffic[END_REF]. Free flow traffic conditions result in neatly separated speeds on different lanes, as shown, e.g., in Fig. 3(a). Only in rare occasions, we record mild congestion situations that result in a pinch effect, i.e., travel speeds on different lanes are leveled down to a same value [START_REF] Kerner | The physics of traffic[END_REF]. An example of such a phenomenon is observed in Fig. 3(b), at around 8:50 a.m.

RT n° 4 Road traffic dataset generation

Our objective is to generate road traffic traces that are representative of unidirectional free flow highway traffic in a quasi-stationary state, i.e., such that comparable traffic conditions are present in between the in-flow and out-flow boundaries of the simulated highways. The datasets should cover a geographical distance that is sufficient for networking studies (e.g., for warning message dissemination or floating car data upload). In the following, we simulate a 10-km road segment, but we remark that our approach can easily accommodate shorter or longer highway stretches.

In order to achieve the goal above, we fed the real-world traffic count data presented in Sec. 3 to a microscopic vehicular mobility simulator. Next, we first present the state-of-art car-following and lane-changing models employed by the simulator. Then, we discuss the calibration of the microscopic models to attain quasi-stationary road traffic conditions. Finally, we outline the features of the resulting dataset.

Microscopic vehicular mobility models

The car-following and lane-changing microscopic mobility models implemented by our simulator are IDM and MOBIL, both of which have been validated by the transportation research community and are among those most commonly used in vehicular networking research.

The Intelligent Driver Model (IDM) [START_REF] Treiber | Congested traffic states in empirical observations and microscopic simulations[END_REF] characterizes the behavior of the driver of a car i through her instantaneous acceleration dv i (t)/dt, which is calculated as

dv i (t) dt = a 1 - v i (t) v max i 4 - ∆x des (t) ∆x i (t) 2 (1)
∆x des (t) = ∆x saf e + v i (t)∆t saf e - v i (t)∆v i (t) 2 √ ab (2)
In (1), v i (t) is the current speed of vehicle i, v max i is the maximum speed its driver would like to travel at, and ∆x des (t) is the so-called desired dynamical distance, representing the distance that the driver should keep from her leading vehicle. The latter is computed in (2) as a function of several measures taken with respect to the car in front of vehicle i: the minimum bumperto-bumper distance ∆x saf e , the minimum safe time headway ∆t saf e , and the speed difference ∆v i (t). In both equations, a and b denote the maximum absolute acceleration and deceleration, respectively. When combined, these formulae return the instantaneous acceleration of the car, as a combination of the desired acceleration [1 -(v i (t)/v max i

) 4 ] on a free road, and the braking deceleration induced by the preceding vehicle (∆x des (t)/∆x i (t)) 2 .

The Minimizing Overall Braking Induced by Lane-changes model (MOBIL) [START_REF] Treiber | Realistische mikrosimulation von strassenverkehr mit einem einfachen modell[END_REF], follows a game theoretical approach, and allows the driver of a vehicle i to move to an adjacent lane if her advantage in doing so is greater than the disadvantage of the trailing car j in the new lane. The (dis)advantage is measured in terms of acceleration, which translates the lane movement condition above into the inequality

dv i (t) dt L - dv i (t) dt + a L ≥ p dv j (t) dt - dv j (t) dt L + ka, (3)
where the notation | • | L denotes accelerations computed as if vehicle i were traveling on the lane to its left rather than in the current one. In (3), p ∈ (0, 1] is a politeness factor that models the selfishness of the driver with respect to the new back vehicle j, ka is a hysteresis threshold (computed as a fraction k of the maximum acceleration a) that prevents lane hopping, and a L is Inria a bias acceleration that can be used to favor or limit movements to left. An identical formulation can be used for right-hand-side lane changes, and the respective advantages can be compared to determine the final lane movement, if any. Note that, in Spain, road traffic regulations enforce drivers to travel on the rightmost lane whenever possible: we thus expect a R > a L and a R > 0, i.e., right-hand-side movements to be favored over left or no movement, if equivalent conditions are present on all lanes.

Parameter calibration

In order to obtain quasi-stationary traffic conditions over the simulated highway segment, some calibration of the IDM and MOBIL parameters proved necessary. Specifically, for the acceleration a, deceleration b, politeness factor p and minimum safety distance ∆x saf e we could use the default values suggested in [START_REF] Treiber | Congested traffic states in empirical observations and microscopic simulations[END_REF][START_REF] Treiber | Realistische mikrosimulation von strassenverkehr mit einem einfachen modell[END_REF]. The other parameters had instead to be adapted to the specificities of each road scenario, as detailed next. We remark that we employed a discrete time step of 100 ms in all mobility simulations. The final IDM and MOBIL settings are summarized in Tab. 1. Maximum desired speed. Vehicles enter the simulated highway segment at the time and with the speed defined by the real-world traffic count data. However, it is then necessary to determine the maximum desired speed v max i of each driver, i.e., the cruise velocity that she would keep if she were alone on the highway. To that end, we exploit our observation in Sec. 3.2 that all traffic count datasets yield free flow conditions. We thus assume that the experimental distribution of real-world ingress speeds can be used as a base for inferring the actual distribution of desired speeds.

We thus computed the Probability Density Functions (PDF) of traffic count ingress speeds for each combination of road, day, and measurement hour. We found that the empirical distributions for the same highway at the same time were overlapping, and thus aggregated such data for different days. This left us with four combinations of highway and hour. For each such combination, we separated the results on a per-lane basis, since each lane yields significantly diverse speeds.

The resulting PDFs are shown in Fig. 4. We note that the empirical ingress velocities tend to follow gaussian distributions in all cases. Nonlinear least squares fittings performed via the Trust-Region algorithm confirm that intuition, with adjusted R-square values always above 0.97 and root mean square errors below 0.02. The fitted gaussian distributions are portrayed as solid lines in the plots of Fig. 4.

Despite being measured in free flow conditions, ingress speeds cannot be directly mapped to desired velocities. Indeed, the presence of surrounding vehicles implicitly limits the ingress speeds, even if by a limited extent. Therefore, when the values of v max i are directly extracted from the gaussian ingress speed distributions, vehicles tend to have an egress velocity from the simulated road segment that is on average 2.8 m/s (around 10 km/h) lower than expected. By adding such a speed offset to the mean of the theoretical PDF, we obtain, for each lane of every highway and hour combination, a distribution of maximum desired speed f V (v) that grants quasi-stationary speed traffic over the whole highway stretch. Interestingly, we highlight that the shapes of such f V (v) are similar for the two measurement hours on a same highway, but remarkably different for the two highways3 : this matches the intuition that the desired speed of a driver is influenced by the type of road she is traveling on, rather than by daytime. Safety time. The safe time headway ∆t saf e was reduced from the default value of 1.5 s to 0.65 s. As a matter of fact, we observed that in the real-world traffic count data the inter-arrival time between two cars on a same lane may get much smaller than 1.5s. Maintaining the default value of ∆t saf e can then bias the IDM behavior in cars that just entered the scenario. Specifically, we recorded excessive brakings that caused traffic jams at the entrance, which in turn prevented new cars to enter the highway segment as commended by the real-world data. After extensive tests, we selected ∆t saf e = 0.65 s, as that is the largest safety time for which the effect described above disappears, allowing vehicles to enter the simulated scenario according to the traffic count datasets. Lane change bias and hysteresis threshold. The default MOBIL settings resulted in a traffic that was highly skewed towards the left lane, which became then unrealistically congested. We run comprehensive trials to find the combination of right (a R ) and left (a L ) lane change bias, and lane change hysteresis threshold factor (k) that granted quasi-stationary traffic over the different lanes. The most consistent ingress and egress per-lane properties were obtained for a R = 0.2 m/s 2 , a L = 0 m/s 2 , and k = 0.3. Those values also respect the intuitive observation that traffic should be skewed to the right so as to abide Spanish road regulations, and we thus chose to use them for the generation of the traces.

Final dataset

The final road traffic dataset 4 is composed of 16 different mobility traces, resulting from all possible combinations of two highways, four days and two measurement times. By default, the position of each vehicle is recorded over a 10-km road stretch at every second. Each mobility trace has a duration of 30 minutes, dictated by the measurement intervals of the real-world traffic For instance, Fig. 5 presents heatmaps of the average vehicular density recorded on two sample traces: we can remark that traffic conditions are similar all over the road stretch (along the x axis of the plots) for each lane (y axis of the plots). Finally, Tab. 2 summarizes the main properties of the 16 traces, in terms of average values of speed, vehicular density, and headway time recorded over the 30-minute duration of each trace. As already anticipated, the speeds measured on each highway at the same time of different days tend to be similar. However, different highways and hours yield significant diversity in terms of average speeds, which range between 68 km/h (A6 at 8:00 a.m.) and 87 km/h (M40 at 11:30 a.m.). Also, we remark in all traces the correlation between the speed and the other metrics, which is consistent with well-known road traffic theory results [START_REF] Kerner | The physics of traffic[END_REF]. However, M40 and A6 show non-comparable behaviors in that sense, since lower densities and longer headway times on A6 can lead to a traffic that -counter-intuitively -is much slower than that recorded on M40 in presence of higher density and shorter headway times. This is an effect induced by the different nature of the two highways, and a proof of the heterogeneous nature of the traces in our vehicular mobility dataset.

Connectivity analysis

In this section, we study the connectivity properties of the vehicular networks built on the mobility traces derived in the previous sections. We start by presenting the model that we apply and the metrics that we consider, followed by the results we obtained. We conclude the section by discussing vehicle-to-vehicle contact duration properties.

Network model and metrics

In our analysis, we focus on the instantaneous connectivity properties of the network. For each time instant, we represent the vehicular network as an undirected graph G(V(t), E(t)), where V(t) = {v i (t)} is the set of vertices v i (t), representing each a vehicle i appearing in the network at time t, and E(t) = {e ij (t)} is the set of edges e ij (t). We consider that an edge e ij (t) is present between v i (t) and v j (t), if a communication link exists between vehicles i and j at time t. We adopt the unit disc model as a representation of the RF signal propagation. Thus, we assume that two vehicles can communicate directly if they are separated by a distance of at most R meters, where R is the communication range. We apply this simplistic model due to the fact that deterministic propagation models do not scale well to large mobile scenarios and imply expensive instantaneous computations. Moreover, the unit disc model is representative of the average behavior of the system obtained under stochastic propagation models. Specifically, in the following we present results for different values of R, from 50 to 300 meters. Despite standards claiming up to 1-km communication ranges, these values represent the lower and upper bounds for acceptable packet delivery ratios, according to different experimental studies on Dedicated Short Range Communication (DSRC) transfers [START_REF] Martelli | A measurement-based study of beaconing performance in IEEE 802.11p vehicular networks[END_REF][START_REF] Bai | Toward understanding characteristics of dedicated short range communications (DSRC) from a perspective of vehicular network engineers[END_REF].

RT n°

We use the graph model to define several metrics of interest to the instantaneous connectivity analysis. Namely, a component

C k (t) = G(V k (t), E k (t)) is a subgraph of G(V(t), E(t)), where V k (t) is a subset of V(t)
and includes a group k of vehicles that can communicate together via direct or multi-hop links at time t, and E k (t) is a subset of E(t) formed by all the communication links among the vehicles in k.

S k (t) = V k (t) is the size of the component C k (t). C(t) = {C k (t)}
refers to the set of components appearing in the network at time t, and C(t) = C(t) denotes their number. We use C max (t) to indicate the largest cluster appearing in the network at time t and S max (t) to represent its size. We denote the number of nodes appearing on the highway at t as N (t) = |V(t)|. Finally, the average size of clusters appearing at time t is referred to as S avg (t) = N (t)/C(t).

In the rest of the paper, for the sake of simplicity, we drop the time notation and we refer to all metrics at a generic time instant: N represents then the number of vehicles on the highway, C the number of components, S avg the average size of component, and S max the largest component size.

Instantaneous network connectivity

We start by studying the global connectivity properties of the network at each time instant. Thus, we focus on the distributions of the number of components C and of the size of the largest component S max . Indeed, C is a measure of how fractioned the network is, while S max is the maximum number of nodes that can be reached via multi-hop communication at a given time instant. Clearly, the lower is C and the larger is S max , the better connected is the vehicular network.

In Fig. 6, we present the distributions of C (left plots) and S max (right plots), for various values of the communication range R (on different rows). In every plot, each candlestick refers to one mobility trace, and is obtained by aggregating the metrics computed in all instantaneous graphs observed at every second during the 30-minute timespan of that trace. The lowest and highest values in the candlesticks (in grey) are the minimum and maximum values that the metric attains. The inner error bars (in red) depict the first and ninth deciles, while the box highlights the first and third quartiles (in red) around the median value (in black). Also, the step function in the S max plots is the maximum value N max of N , observed in all instances of the respective trace. It thus represents the upper bound to S max : the closer is S max to N max , the closer the vehicular network is to a fully connected single component.

Let us now focus on the case of R = 50 m, in Fig. 6(a) and Fig. 6(b). We observe that the connectivity metrics are very similar on different days of the week, forming four clear clusters in both plots. This implies that we can expect the network to have consistent topological properties from Monday to Friday. The aforementioned clusters are instead separated by highway (A6 or M40) and hour (8:30 or 11:30). Specifically, the network appears more fragmented on A6 than on M40, and, in both cases, the connectivity is worse at 11:30 a.m. than at 8:30 a.m. Intuitively, we can ascribe such a variability to the different road traffic conditions encountered on the two highways at different times of the day. In fact, we will further investigate the issue and provide a more rigorous explanation later in this section.

Most of the observations above still hold when considering different communication ranges separately. Nonetheless, when comparing the plots for varying R, it is evident that increasing the communication range can dramatically improve the connectivity of the network. Indeed, for R = 50 m, there exist on average 30-50 disconnected components, the largest of which only comprises around one tenth of the vehicles. As R grows, however, the network fragmentation is reduced, and more nodes join the largest component. When R = 300 m almost all vehicles belong all the time to one single component, as in Fig. 6(g) and Fig. 6(h). We conclude that the communication range is the foremost parameter controlling the vehicular network connectivity, as it can induce variations in C and S max that are much more significant than those imputable to road traffic conditions. Indeed, R = 300 m guarantees a well connected network independently of the traffic scenario, but reducing that value causes the topology to rapidly break apart. In this latter case, and precisely when R falls below 200 m, differences among the traces emerge. In order to better investigate the causes behind such a diversity, we analyze time series rather than temporal aggregates. More precisely, we observe the evolution of N , S max , S avg and C over time, for each trace and varying R. In Fig. 7, we show two sample plots for one of the A6 traces at 8:30 with R = 50 m and R = 100 m. We omit results for R = 200 m and R = 300 m because they lead, as seen before, to a stable network with constant behaviors over time (i.e., C ∼ 1, and S max ∼ S avg ∼ N ). The key remark here is that the time series of C, S max and S avg are combinations of two dynamics: slow variations somehow similar to shadowing in a RF signal at the receiver, and rapid variations comparable to fast fading in the RF signal analogy. Next, we analyze these two dynamics in detail. Slow connectivity dynamics. In Fig. 7, we can glimpse that the slow dynamics of C and S max seem correlated to that of N . To explore that possibility, in Fig. 8 we draw scatterplots The correlation is evident in all cases. The Pearson product-moment correlation coefficients (PPMCC) [START_REF] Karl | Notes on the history of correlation[END_REF] computed on the different scatterplots, in Tab. 3, confirm that6 . Specifically, we remark the strong negative correlation between C and N for 50-m and 100-m communication ranges. Such a correlation is reduced when R = 200 m and R = 300 m, since the number of clusters becomes then very small and nearly constant. For S max , we note a high positive correlation under all communication ranges. We conclude that the slow topological dynamics of the network are guided solely by the vehicular density, with more intense road traffic leading to a smaller number of larger components of connected vehicles. This is a key finding, since it unveils how long-term variations of the connectivity are not affected by the type of highway, the weekday or the hour considered, but only by the volume of traffic they yield. It is also interesting to note that this result assimilates highway and metropolitan scenarios, as a akin conclusion was reached for urban vehicular networks [START_REF] Uppoor | Generation and analysis of a large-scale urban vehicular mobility dataset[END_REF].

Fast connectivity dynamics The fast dynamics observed in Fig. 7 lead to changes in the connectivity metrics at every few seconds, and are also the cause of the dispersion of values around the linear correlation models in Fig. 8. We speculate that such variations are due to fastspeed vehicular movements that cause links to be established and tore down with high frequency.

To investigate this prospect, we compute the distributions of the durations of vehicle-to-vehicle (V2V) contacts in the mobility traces. Due to space limitations, we show the results for one trace only, as all other traces displayed very similar contact duration properties. Fig. 9 portrays the Probability Density Function (PDF), the Cumulative Distribution Function (CDF), and the Complementary CDF (CCDF) of the contact duration between any pair of communicating vehicles over M40, on Mon 11 at 11:30 a.m., with varying R.

As presumed, contacts are relatively short, even if we are considering unidirectional road traffic that should favor long-lived connections 7 . The CDF shows that 80% of the contacts last less than 5 minutes when R = 300 m, and less than 1 minute when R = 50 m. More importantly to us, there are multiple elements that link V2V contact durations to the fast connectivity dynamics we previously observed. First, the probability peak in the PDF for R = 50 m is attained at 30 s, i.e., it is in the same order of magnitude of the fast temporal dynamics in Fig. 7(a). Second, the relationship above holds also as R grows, since the fast temporal dynamics become slower in Fig. 7(b) exactly as the V2V contacts become more stable. Third, the effect of higher values of R is milder than one could expect, and the majority of the links last well below one minute at R = 100 m. As a result, increasing R generates a much larger number of slightly more stable links, i.e., offers larger space for connectivity disruptions: this is in agreement with the fact that fast connectivity dynamics are stronger in amplitude in Fig. 7(b) than in Fig. 7(a). Overall, these elements let us conclude that the V2V contact duration is the cause behind the short-term variations of the vehicular network connectivity.

Conclusions

In this paper, we analyzed real-world high-resolution road traffic count measurements, and employed them as the cornerstone for the generation of realistic synthetic traces of vehicular mobility on highways. The resulting datasets are representative of quasi-stationary unidirectional traffic in heterogeneous conditions, including different highways, weekdays and measurement hours. They are the current state of art in publicly available highway traffic datasets for network study.

The topological analysis carried out on the mobility datasets reveals that the factor that primarily controls the connectivity of highway vehicular networks is the communication range. Notably, a 300-m communication range appears to grant full connectivity under any condition; when instead the range falls below 200 m, the network connectivity undergoes temporal variations that are a combination of slow dynamics, controlled by the road traffic density, and fast dynamics, due to the predominant presence of short-lived contacts. Interestingly, we found such properties to be invariant throughout multiple weekdays and on different highways.

Figure 1 :

 1 Figure 1: (a): Geographical location of the two highways near Madrid,Spain. (b): Close-by view of measurement point on M40. (c): Close-by view of measurement point on A6. (d): Close-by view of the induction loop on M40.

Figure 2 :

 2 Figure 2: Inter-arrival times CDF measured on May 12, 2010. Each plot refers to a different lane of M40 at 8:30 a.m. (a, b, c) and A6 at 11:30 a.m. (d, e, f). Solid black lines show the mixture model for each distribution. The fitting is good under any combination of highway, day and hour for each lane.

  Mon 10,[START_REF] Uppoor | Generation and analysis of a large-scale urban vehicular mobility dataset[END_REF]:30 a.m.

Figure 3 :

 3 Figure 3: Per-lane speed versus time in two dataset samples. The plots are representative of free flow traffic (a) and of a mild perturbation (b) on M40.

Figure 4 :

 4 Figure 4: Empirical per-lane ingress speed distributions and fitted normal distributions for different highway and measurement hour combinations.

Figure 5 :

 5 Figure 5: Heatmaps of the average vehicular density in two sample traces, on Wednesday the 12th: M40 at 11:30 a.m. (top) and A6 at 8:30 a.m. (bottom)

Inria

Figure 6 :

 6 Figure 6: Distributions of the number of components C (left), and of the largest component size S max (right), for each trace under different values of R.

Figure 7 :

 7 Figure 7: Time series of N , S max , S avg , and C on the A6 highway at 8:30 a.m. of Friday the 7th, for R = 50 m (a) and R = 100 m (b).

Figure 8 :

 8 Figure 8: Scatterplots of C and S max versus N , aggregated over all traces.

Table 3 :Figure 9 :

 39 Figure 9: Vehicular contact duration: PDF (left), CDF (middle), and CCDF (right) for M40, on Tue 11 at 11:30 a.m. Curves map to different R values.of such measures 5 by aggregating data from all of the 16 traces, when R = 50 m and R = 100 m. There, each point represents one instantaneous connectivity graph in one of the traces. The correlation is evident in all cases. The Pearson product-moment correlation coefficients (PPMCC)[START_REF] Karl | Notes on the history of correlation[END_REF] computed on the different scatterplots, in Tab. 3, confirm that6 . Specifically, we remark the strong negative correlation between C and N for 50-m and 100-m communication ranges. Such a correlation is reduced when R = 200 m and R = 300 m, since the number of clusters becomes then very small and nearly constant. For S max , we note a high positive correlation under all communication ranges. We conclude that the slow topological dynamics of the network are guided solely by the vehicular density, with more intense road traffic leading to a smaller number of larger components of connected vehicles. This is a key finding, since it unveils how long-term variations of the connectivity are not affected by the type of highway, the weekday or the hour considered, but only by the volume of traffic they yield. It is also interesting to note that this result assimilates highway and metropolitan scenarios, as a akin conclusion was reached for urban vehicular networks[START_REF] Uppoor | Generation and analysis of a large-scale urban vehicular mobility dataset[END_REF]. Fast connectivity dynamics The fast dynamics observed in Fig.7lead to changes in the connectivity metrics at every few seconds, and are also the cause of the dispersion of values around the linear correlation models in Fig.8. We speculate that such variations are due to fastspeed vehicular movements that cause links to be established and tore down with high frequency. To investigate this prospect, we compute the distributions of the durations of vehicle-to-vehicle (V2V) contacts in the mobility traces. Due to space limitations, we show the results for one trace only, as all other traces displayed very similar contact duration properties. Fig.9portrays the Probability Density Function (PDF), the Cumulative Distribution Function (CDF), and the Complementary CDF (CCDF) of the contact duration between any pair of communicating vehicles over M40, on Mon 11 at 11:30 a.m., with varying R.As presumed, contacts are relatively short, even if we are considering unidirectional road traffic that should favor long-lived connections7 . The CDF shows that 80% of the contacts last less than 5 minutes when R = 300 m, and less than 1 minute when R = 50 m. More importantly to us, there are multiple elements that link V2V contact durations to the fast connectivity dynamics

Table 1 :

 1 IDM and MOBIL parameter settings

				Model	Parameter Meaning	Value
				IDM	a	Maximum acceleration	1 m/s 2
				IDM	b	Maximum (absolute) deceleration 2.5 m/s 2
				IDM IDM	v max i ∆x saf e	Maximum desired speed Minimum safety distance	∼ f V (v) 1 m
				IDM	∆t saf e	Minimum safe time headway	0.65 s
				MOBIL p	Politeness factor	0.5
				MOBIL a L	Bias acceleration (left)	0 m/s 2
				MOBIL a R	Bias acceleration (right)	0.2 m/s 2
				MOBIL k	Hysteresis threshold factor	0.3
	Empirical Probability Density Function	0.0 0.4 0.0 0.4 0.2 0.2			0.0 0.4 0.2
		40	80	120	160
			Speed (km/h)

Table 2 :

 2 Main average characteristics of the 16 traces

	Trace	Speed	Density	Headway
		(km/h) (vehicles/km)	(s)
	M40, Fri 7, 8:30	85.66	56	0.74
	M40, Mon 10, 8:30	84.85	59	0.72
	M40, Tue 11, 8:30	84.33	63	0.70
	M40, Wed 12, 8:30	84.98	60	0.71
	M40, Fri 7, 11:30	86.88	52	0.82
	M40, Mon 10, 11:30	87.49	49	0.84
	M40, Tue 11, 11:30	87.52	46	0.88
	M40, Wed 12, 11:30	87.71	48	0.85
	A6, Fri 7, 8:30	69.09	62	0.86
	A6, Mon 10, 8:30	68.73	62	0.85
	A6, Tue 11, 8:30	68.24	65	0.83
	A6, Wed 12, 8:30	68.65	61	0.86
	A6, Fri 7, 11:30	77.26	42	1.11
	A6, Mon 10, 11:30	77.01	40	1.14
	A6, Tue 11, 11:30	76.99	43	1.10
	A6, Wed 12, 11:30	77.19	42	1.13

RT n°

http://pems.dot.ca.gov/ Inria

As an example, Dirección General de Tráfico provides an elaboration of the traffic data via the Infocar web service at http://infocar.dgt.es, including historical data visualization for some of the observation points.

This can be observed by looking at the gaussian shapes in Fig.

4, since the final f V (v) are just shifted versions of the same over the x axis.4 Available at http://trullols.site.ac.upc.edu/downloads.Inria

Note that studying the correlation for Smax = N /C would be redundant.

We recall that the value of PPMCC varies between -1 (perfect negative linear correlation) and 1 (perfect positive linear correlation), with values around 0 indicating absence of linear correlation between the variables.

Such long-lived V2V contacts (established by vehicles traveling at similar speeds for a significant amount of time) do exist, but they occur with much lower probability than short-lived V2V contacts. This is demonstrated by the heavy tail of the distributions, shown in the log-log scale CCDF plot of Fig.9.

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-0803

grant from the Rhône-Alpes Region.


M Sadeghi
email: m.saadeghii@gmail.com

M Babaie-Zadeh

C Jutten

Learning Overcomplete Dictionaries Based on

Atom-by-Atom Updating Mostafa Sadeghi * , Massoud Babaie-Zadeh, Senior Member, IEEE, and Christian Jutten, Fellow, IEEE Abstract-A dictionary learning algorithm aims to learn a set of atoms from some training signals in such a way that each signal can be approximated as a linear combination of only a few atoms. Most dictionary learning algorithms use a two-stage iterative procedure. The first stage is to sparsely approximate the training signals over the current dictionary. The second stage is the update of the dictionary. In this paper we develop some atom-by-atom dictionary learning algorithms, which update the atoms sequentially. Specifically, we propose an efficient alternative to the well-known K-SVD algorithm, and show by various experiments that the proposed algorithm has much less execution time compared to K-SVD while its results are better. Moreover, we propose a novel algorithm that instead of alternating between the two dictionary learning stages, performs only the second stage. While in K-SVD each atom is updated along with the nonzero entries of its associated row vector in the coefficient matrix (which we name it its profile), in the new algorithm, each atom is updated along with the whole entries of its profile. As a result, contrary to K-SVD, the support of each profile can be changed while updating the dictionary.

To further accelerate the convergence of this algorithm and to have a control on the cardinality of the representations, we then propose its two-stage counterpart by adding the sparse approximation stage. We evaluate the performance of the proposed algorithms by performing two sets of experiments. The first set is the reconstruction of a true underlying dictionary, while the second set is designing a sparsifying dictionary for a certain class of signals. The results emphasize on the promising performance of the proposed algorithms.

Index Terms-Sparse approximation, dictionary learning, compressive sensing, K-SVD.

I. INTRODUCTION

A. Sparse Signal Approximation

S PARSE decomposition of signals based on some basis functions has attracted a lot of attention during the last decade [START_REF] Elad | Sparse and Redundant Representations[END_REF]. The problem consists in approximating a given signal as a linear combination of as few as possible basis functions. In this context, each basis function is called an atom and their collection as the columns of a matrix is called dictionary [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF]. The dictionary may be overcomplete, i.e., the number of atoms may be (much) more than the dimension of the atoms. Specifically, let y ∈ R n be the signal which is going to be sparsely represented in the dictionary D ∈ R n×K with K > n. This amounts to solve the following problem, P 0 : min

x ∥x∥ 0 subject to y = Dx, (1)
where ∥.∥ 0 stands for the so-called ℓ 0 pseudo-norm which counts the number of nonzero elements. Because of the combinatorial nature of P 0 [START_REF] Elad | Sparse and Redundant Representations[END_REF] that makes it hard to solve, especially in high dimensions, some relaxation methods have been introduced. A well-known relaxation approach is to use the ℓ 1 norm instead of the ℓ 0 pseudo-norm. This leads to the following convex P 1 problem which is known as Basis Pursuit (BP) [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] P 1 : min

x ∥x∥ 1 subject to y = Dx. (2)
In practical situations instead of the exact equality y = Dx, the constraint ∥y -Dx∥ 2 ≤ ϵ is used, which has a denoising and stabilizing effect [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF], [START_REF] Babaie-Zadeh | On the stable recovery of the sparsest overcomplete representations in presence of noise[END_REF]. In this case, usually the term approximation is used instead of representation.

Many algorithms have been introduced to solve the problem of finding the sparsest approximation of a signal in a given overcomplete dictionary (as a good review, see [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF]). These methods can be classified into two general categories, greedy methods such as Orthogonal Matching Pursuit (OMP) [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF], Stagewise Orthogonal Matching Pursuit (StOMP) [START_REF] Donoho | Sparse solution of underdetermined systems of linear equations by stagewise orthogonal matching pursuit[END_REF], Compressive Sampling Matching Pursuit (CoSaMP) [START_REF] Needell | Cosamp: Iterative signal recovery from in-complete and inaccurate samples[END_REF], Subspace Pursuit (SP) [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF], and relaxation methods, which replace the combinatorial P 0 problem with a tractable one, e.g., P 1 problem. Iterative Shrinkage-Thresholding (IST) [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF], [START_REF] Elad | A wide-angle view at iterated shrinkage algorithms[END_REF], [START_REF] Becker | NESTA: a fast and accurate firstorder method for sparse recovery[END_REF], Iterative Hard-Thresholding (IHT) [START_REF] Blumensath | Accelerated iterative hard threshoding[END_REF], Iteratively Reweighted Least Squares (IRLS) [START_REF] Chartrand | Iteratively reweighted algorithms for compressive sensing[END_REF], Smoothed ℓ 0 (SL0) [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF], and interior-point methods [START_REF] Kim | An interior point method for large-scale ℓ 1 -regularized least squares[END_REF] are some examples of the second category. Greedy algorithms successively choose the appropriate atoms of the dictionary that result in the greatest improvement in the quality of the approximation. These algorithms benefit from having high speed, but their accuracy is usually less than that of the second category.

Various image processing tasks (e.g., denoising, compression, inpainting, zooming) [START_REF] Elad | On the role of sparse and redundant representations in image processing[END_REF], Blind Source Separation (BSS) in underdetermined mixtures [START_REF] Bofill | Underdetermined blind source separation using sparse representations[END_REF], Compressive Sensing (CS) [START_REF] Donoho | Compressed sensing[END_REF], [START_REF] Baraniuk | Compressive sensing[END_REF], decoding real field codes [START_REF] Candès | Decoding by linear programming[END_REF], linear regression and variable selection [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] are some applications where the sparse approximation approach has already been applied.

B. Learning overcomplete dictionaries

For a given class of signals, e.g., class of natural facial images, the dictionary should have the capability of sparsely representing the signals. To this aim, atoms of the dictionary should capture the most salient features of the signals. In some applications there are predefined and universal dictionaries which are known to be well-matched to the contents of the given class of signals, for example the overcomplete DCT dictionary for the class of natural images. These non-adaptive dictionaries are appealing because of their simplicity and in some cases their fast computations. However, learning the atoms from a set of training signals would result in dictionaries with the capability of better matching the contents of the signals. In this way, the adaptive dictionaries would outperform the non-adaptive ones in many applications such as image denoising [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF], image compression [START_REF] Bryt | Compression of facial images using the K-SVD algorithm[END_REF], classification tasks [START_REF] Mairal | Task-driven dictionary learning[END_REF], and so on.

Most dictionary learning algorithms are indeed a generalization of the K-means clustering algorithm [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF]. While in K-means each training signal is forced to use only one atom (cluster center) as its approximation, in the dictionary learning problem each signal is allowed to use more than one atom provided that it uses as few as possible atoms. The approach of K-means to optimize a set of atoms (called codebook) is to iteratively perform two stages [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF]. In the first stage, known as the clustering stage, training signals are assigned to their nearest atoms (usually in the sense of the ℓ 2 norm). The second stage is the update of the atoms in which each atom is updated as the average of the signals in its cluster. This procedure is repeated several times. As a generalization of this approach, dictionary learning algorithms iteratively perform the two stages of sparse approximation and dictionary update. In the first stage, which is actually the clustering of the signals into a union of subspaces, the sparse approximations of the signals are computed using the current dictionary. The second stage is the update of the dictionary. Up to our best knowledge, most dictionary learning algorithms differ mainly in the way of updating the dictionary [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF], [START_REF] Engan | Method of optimal directions for frame design[END_REF], [START_REF] Yaghoobi | Dictionary learning for sparse approximations with the majorization method[END_REF]. Some algorithms such as K-Singular Value Decomposition (K-SVD) [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF] are based on updating the atoms one-by-one, while some others such as Method of Optimal Directions (MOD) [START_REF] Engan | Method of optimal directions for frame design[END_REF] update the whole set of atoms at once.

In another point of view, dictionary learning algorithms are divided into two groups: online [START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF], [START_REF] Skretting | Recursive least squares dictionary learning algorithm[END_REF], [START_REF] Labusch | Robust and fast learning of sparse codes with stochastic gradient descent[END_REF] and batch-based [START_REF] Engan | Method of optimal directions for frame design[END_REF], [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF], [START_REF] Yaghoobi | Dictionary learning for sparse approximations with the majorization method[END_REF] algorithms. Online algorithms continuously update the dictionary using only one (or a mini batch of) training data. Because of this, they enjoy from having the capability of handling very large sets of signals, which is common in, for example, image processing tasks (see the experiments of [START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF] on a high-resolution image). Recursive Least Squares dictionary learning algorithm (RLS-DLA) [START_REF] Skretting | Recursive least squares dictionary learning algorithm[END_REF] is an example of online algorithms, which can be considered as a generalization of the McQueen variant of K-means [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF]. One major benefit of online algorithms is to learn gradually, which is denoted as the scaling the past data in [START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF], and by using a forgetting factor in [START_REF] Skretting | Recursive least squares dictionary learning algorithm[END_REF]. Batch-based algorithms use the whole set of training data to update the dictionary. This increases the computational burden of these algorithms in high dimensions. However, an advantage of batch-based algorithms over online algorithms is their low computational loads and memory requirements in relatively low dimensions. This is due to the fact that at each sparse approximation stage, the whole set of signals are approximated in the same dictionary. So, the sparse coding algorithm can be optimized to avoid performing common operations. Batch-OMP [START_REF] Rubinstein | Efficient implementation of the K-SVD algorithm using batch orthogonal matching pursuit[END_REF] is an example of such optimized sparse coding algorithms.

C. Our contributions and the structure of the paper

The focus of this paper is on the algorithms that update atoms one-by-one. We first propose an efficient way of updating each atom along with its associated row vector in the coefficient matrix: we call this idea parallel atom-updating. From this, we propose an efficient and fast alternative for K-SVD: we call this algorithm PAU-DL, and show by various experiments that it performs better than K-SVD while its computational burden is substantially lower than it. We then propose a novel method of dictionary learning whose structure is different from the existing algorithms. Specifically, instead of alternating between the two dictionary learning stages, the proposed algorithm performs only the dictionary update stage. We call this algorithm OS-DL. In each alternation of OS-DL, each atom is updated along with the whole entries of its associated row vector in the coefficient matrix. The main differences between the new algorithm and the current atomby-atom dictionary learning algorithms, are as follows. Firstly, the new algorithm does not explicitly perform the sparse approximation stage. Secondly, in this algorithm the support of each row vector of the coefficient matrix is allowed to change.

To further accelerate the convergence of this method and have a control on the cardinality of the representations, we propose to add the sparse approximation stage to this algorithm, and hence we derive another new algorithm, called APrU-DL. OS-DL and APrU-DL use the parallel atom-updating approach. Experimental results show the promising performance of the proposed algorithms relative to K-SVD.

The paper is organized as follows. Section II is devoted to the formulation of the dictionary learning problem. In Section III we describe our proposed algorithms. This section begins with introducing the parallel atom-updating approach; a general idea for updating the atoms of the dictionary oneby-one. Using this idea in the dictionary update stage, we derive the parallel atom-updating dictionary learning (PAU-DL) algorithm. We then propose the two novel atom-by-atom dictionary learning algorithms, i.e., OS-DL and APrU-DL. Section IV presents some experimental results.

D. Notation

We use the following notations. For vector and matrix valued quantities we use small and capital bold face characters, respectively. The superscript "T " denotes vector or matrix transpose. We denote the ith column vector of a matrix X as x i , and its ith row vector as x T

[i] . ∥X∥ F = √ ∑ ij x 2 ij and ∥x∥ p = ( ∑ i |x i | p )
1 p denote the Frobenius norm of the matrix X, and the ℓ p norm of the vector x, respectively. In an iterative algorithm, we denote a parameter in the kth iteration with the iteration number in parenthesis, e.g., d (k) . For a matrix X with L columns and ω ⊆ {1, 2, . . . , L}, we define X(:, ω) as a matrix containing those columns of X that their indices are in ω. Also, x(ω) is a vector containing those entries of x that are indexed by ω. Finally, |ω| denotes the cardinality of ω, that is, its number of entries.

II. DICTIONARY LEARNING PROBLEM

Let {y i } L i=1 be a set of L training signals in R n . Putting these signals as the columns of the matrix Y, the general dictionary learning problem is then to find a sparsifying dictionary, D, by solving the following problem

min D∈D,X∈X ∥Y -DX∥ 2 F , (3)
where D and X are admissible sets of the dictionary and the coefficient matrix, respectively. D is usually defined as the set of all dictionaries with unit column-norms. Since we require that each signal has a sparse approximation, X is the set of all matrices X with sparse columns (see Fig. 1). The sparsity measure can be, for example, the non-convex ℓ 0 pseudo-norm, or the convex ℓ 1 norm. The general approach to solve (3) is to use alternating minimization over X and D, i.e., by fixing D, the objective function is minimized over X, and vice versa.

A. Sparse Approximation

With a fixed D, the minimization of (3) with respect to X is equivalent to sparsely approximating the training signals over D. Among the various sparse coding algorithms, OMP (or its variants) is very appealing for this stage. This is due to two reasons. The first reason is the high speed of OMP compared to relaxation algorithms. The second one is its capability to be efficiently implemented in the batch mode. This fact along with the use of Cholesky factorization result in the significant acceleration of OMP, which leads to Batch-OMP algorithm [START_REF] Rubinstein | Efficient implementation of the K-SVD algorithm using batch orthogonal matching pursuit[END_REF].

Another option, which we use in this paper, is the IST algorithms. Although the speed of these algorithms is lower than that of OMP, their accuracy is better.

B. Dictionary Update

As stated previously, dictionary learning algorithms differ mainly in the way they perform the second stage, i.e., the dictionary update. In this stage, some desired properties about the dictionary may be applied, such as having bounded Frobenius norm [START_REF] Yaghoobi | Dictionary learning for sparse approximations with the majorization method[END_REF], and having bounded self-coherence (i.e., pairwise similarity of the atoms) [START_REF] Sigg | Learning dictionaries with bounded self-coherence[END_REF]. However, the constraint of having unit column-norms is usually used in this stage.

In spite of its name, in the dictionary update stage some information about the coefficient matrix may be updated, too. This information may be the nonzero coefficients of X (as in K-SVD and [START_REF] Dai | Simultaneous codeword optimization (SimCO) for dictionary update and learning[END_REF]) or the whole coefficients of it (as in two of our suggested algorithms). Furthermore, as stated previously, one of the differences among various dictionary learning algorithms is that they either update the whole set of atoms at once (as in MOD) or each atom separately (as in K-SVD).

C. MOD and K-SVD

MOD [START_REF] Engan | Method of optimal directions for frame design[END_REF] is one of the simplest dictionary learning algorithms which firstly finds the unconstrained minimum of ∥Y -DX∥ 2 F and then projects the solution onto the set D. This leads to the closed-form expression

D = YX T (XX T ) -1 , (4)
followed by normalizing the columns of D. K-SVD [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF] is one of the most successful dictionary learning algorithms. In its dictionary update stage, only one atom is updated at a time. Moreover, while updating each atom, the nonzero entries in the associated row vector of X are also updated. In other words, only those signals that have used a specific atom participate in updating that atom. This is in accordance with the approach of K-means in which each atom is updated using its own cluster signals. As stated in [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF], this also prevents each row vector in X to be filled and thus violating the sparsity constraint of the coefficient matrix.

Assume that we want to update the ith atom, d i , along with the nonzero entries of x T

[i] , the ith row of X. We define

ω i = { j : x T [i] (j) ̸ = 0 } as the support of x T [i] .
Then the problem of updating d i along with x T

[i] (ω i ) amounts to solve the following minimization problem

min d,xr ∥E r i -dx T r ∥ 2 F subject to ∥d∥ 2 2 = 1, (5)
where

E r i = E i (:, ω i ), in which E i = Y - ∑ j̸ =i d j x T [j]
denotes the approximation error matrix when d i is removed, and x T r is a row vector of length |ω i |. The above problem is in fact equivalent to finding the closest rank-1 approximation to E r i , which can be easily solved via SVD of E r i . For more details refer to [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF].

III. ATOM-BY-ATOM DICTIONARY UPDATE

In this section we propose three algorithms that update the dictionary atom-by-atom. We call the ith row vector of the coefficient matrix, the profile of the ith atom, because this vector indicates which signals use this atom in their representations. We first introduce the idea of parallel atomupdating and from which we then propose the parallel atomupdating dictionary learning (PAU-DL) algorithm, which is introduced to overcome the computational burden of K-SVD. We then proceed with introducing a novel method for dictionary learning which is based on performing only the second stage of the general dictionary learning procedure, and thus we name it One-Stage Dictionary Learning (OS-DL). In OS-DL, the support of each profile is allowed to change during the update of the dictionary. To further accelerate the convergence rate of OS-DL (as will be shown in the simulations) and control the sparsity level of the representations, we propose an algorithm, which we call Atom-Profile Updating Dictionary Learning (APrU-DL), which in fact adds the first stage of dictionary learning to OS-DL. In all of these algorithms, we choose D to be the set of all unit column-norm dictionaries in R n×K .

A. Parallel Atom-Updating Dictionary Learning (PAU-DL)

The main drawback of K-SVD is its computational burden especially in high dimensions. This is due to performing SVD for atom updating. An alternative way of solving ( 5) is to use the idea of alternating minimization [START_REF] Rubinstein | Efficient implementation of the K-SVD algorithm using batch orthogonal matching pursuit[END_REF]. In other words, ( 5) is alternatively minimized over d and x r . A few (e.g., 3) alternations give a fast approximation to SVD. The resulting algorithm is known as the Approximate K-SVD (AK-SVD) [START_REF] Rubinstein | Efficient implementation of the K-SVD algorithm using batch orthogonal matching pursuit[END_REF]. Although performing more alternations gives a better approximation, the average performance will not exceed the performance of the exact solution, i.e., via SVD.

In this subsection we describe a different way of performing alternating minimization to update the atoms and their profiles. To this aim, consider the overall error matrix,

E = Y -(A 1 + A 2 + . . . + A K ), ∀i : A i = d i x T [i] . (6
) In K-SVD (or AK-SVD), in order to update (the nonzero columns of) for example A i , the updated versions of A 1 , . . . , A i-1 are used to compute E i , while A i+1 , . . . , A K have not been yet updated. Keeping this point in mind, we propose to update the atoms in parallel. In other words, instead of fully updating each A i by performing "A" alternations between d i and x r , "A" alternations are performed in such a way that in each alternation all of A i 's are partially updated (using only one alternation). In this way, in the subsequent alternations, all A i 's have been partially updated. As our experimental results in Section IV suggest, parallel updating of the atoms may result in further accelerating the convergence rate and the quality of the final results. In other words, the new algorithm outperforms K-SVD, which is based on exact solving of the rank-1 approximation problem.

To update each A i , we need to compute the error matrix E i . It can be easily seen that this matrix can be updated as follows. The overall error matrix is firstly computed as E = Y -DX using the current dictionary and coefficient matrix.

Then E i = E + A i and after updating A i to A (new) i , the error matrix E is updated as E = E i -A (new) i
. Algorithm 1 gives a description of the dictionary update based on parallel atom-updating. PAU-DL is an alternative to K-SVD that uses this atom-updating procedure. Algorithm 2 gives a complete description of PAU-DL and AK-SVD. By Batch-OMP(Y, D, τ ) we mean the sparse approximation of Y in D and with threshold τ . Depending on the application at hand, τ may be the threshold on the approximation error or the maximum allowed number of atoms in the sparse approximation of each training signal.

Considering Algorithm 2, we see that in K-SVD: A = 1, with a large B, in AK-SVD: A = 1, B = 3, and in PAU-DL:

A = 3, B = 1.
Here the reader may suggest to use A = 3, B = 3. However, as we saw in our simulations, the results are for i = 1, . . . , K do

4:

E i = E + d i x T [i] 5: Update x T [i] 6:
Update d i 7: for i = 1, . . . , K do 8:

E = E i -d i x T [i]
E i = E + d i x T [i]
9:

E r i = E i (:, ω i ) where ω i = { j : x T [i] (j) ̸ = 0 } 10:
for b = 1, . . . , B do 11:

d i = E r i x T [i] (ω i )
12:

d i = d i /∥d i ∥ 2 13:
x

T [i] (ω i ) = d T i E r i 14:
end for 15:

E = E i -d i x T [i]
16:

end for 17: end for 18: Until convergence similar to those of PAU-DL, yet with a higher computational load.

At first glance, one may think that PAU-DL differs from AK-SVD by simply changing the operation orders. This, however, is not just a simple rescheduling. As explained earlier, the main idea behind PAU-DL is to partially update each atom before moving to the next atom. In this way, we have in our disposal more reliable updates of A i 's in order to compute the error matrix associated with the atom we are going to update; see [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF]. Moreover, as will be seen in Section IV, PAU-DL has a sufficiently better performance than AK-SVD and even original K-SVD.

B. One-stage Dictionary Learning (OS-DL)

In [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF], after developing the main steps of K-SVD algorithm, the authors asserted that "Here one might be tempted to suggest skipping the step of sparse coding and using only updates of columns in D, along with their coefficients, applied in a cyclic fashion, again and again. This, however, will not work well, as the support of the representations will never be changed, and such an algorithm will necessarily fall into a local minimum trap."

In this subsection, we describe OS-DL algorithm. This algorithm ignores the sparse approximation stage and perform dictionary learning by updating atoms along with their profiles in a cyclic fashion. To prevent the above mentioned problem, a sparsity constraint on the profiles of the atoms is introduced, which at the same time allows the support of each profile to change (and probably not getting trapped into a local minimum) and prevents each profile to be filled. Allowing the support of each profile to change lets each atom adaptively find its cluster members (see Subsection III-C). Each atom and its profile are updated by solving the following minimization problem ∀i :

{ d i , x [i] } = argmin d,z 1 2 ∥E i -dz T ∥ 2 F + λ∥z∥ 1 , (7)
subject to the constraint ∥d i ∥ 2 2 = 1. Note that if we set λ = 0 and restrict E i and x T

[i] to those training data that have used d i in their approximations, then equation ( 7) is exactly the one used in K-SVD. Note also that the above problem is indeed a regularized rank-1 approximation of E i . It is worth mentioning that in OS-DL there is no control on the cardinalities of the representations, and indeed each one may have a different cardinality. However, in the next algorithm described in Subsection III-C, this problem is avoided by performing the sparse approximation stage.

To solve [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF], we use alternating minimization. At the first stage, we update z with a fixed d, and at the second stage, we update d using the previously updated z. The update formula for z is obtained by solving [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF] with d being fixed. Interestingly, this leads to a simple closed-form formula [START_REF] Elad | Why simple shrinkage is still relevant for redundant representations?[END_REF]. To see this, note that ( 7) can be de-coupled for each entry of z. Then, we should solve L scalar problems of the form

min z f (z) = 1 2 ∥e -zd∥ 2 2 + λ|z|, ( 8
)
where e is the corresponding column of E i . Since the last term in f (z) is non-differentiable, we require that zero is included in the subgradient of f (z), i.e., 0 ∈ ∂f (z). This results in

0 = -d T (e -zd) + λsgn(z). (9)
Now, considering the assumption ∥d∥ 2 2 = 1, the final solution is found using the soft-thresholding operation

z = sgn(e T d). max(0, |e T d| -λ) = S λ (e T d), (10)
where S λ (a) is the well-known soft-thresholding function [START_REF] Elad | Why simple shrinkage is still relevant for redundant representations?[END_REF] defined as follows

S λ (a) ≜    a -λ if a > λ 0 if |a| ≤ λ a + λ if a < -λ . ( 11
)
The final solution is

x [i] = S λ (E T i d i )
, where S λ (.) acts component-wise.

The update formula for d is also of the simple form

d i = E i x [i] , (12)
followed by a normalization. Like PAU-DL, in OS-DL the atoms are updated in parallel. A few iterations (A = 3 in our experiments) of the alternating minimization is sufficient to obtain x T [i] and d i . This process is repeated several iterations until a stopping condition is satisfied. The final algorithm is summarized in Algorithm 3. for i = 1, . . . , K do

6:

E i = E + d i x T [i]
7:

x

[i] = S λ (E T i d i )
8:

d i = E i x [i]
9:

d i = d i /∥d i ∥ 2 10: E = E i -d i x T [i]
11:

end for 12: end for

How does it work?

According to [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF], only those columns of E i are represented by d i that are sufficiently similar to it. After finding these signals, the cluster members of d i are actually found. Similar to K-means, each atom is updated as the (weighted) average of its cluster members, followed by a normalization. Indeed, each atom is updated as a sparse linear combination of the columns of the error matrix, see [START_REF] Elad | A wide-angle view at iterated shrinkage algorithms[END_REF]. In OS-DL, a sparse coding-like stage is implicitly performed 1 . This is because after updating the whole row vectors of the coefficient matrix, its columns are also updated. This way of sparse coding is very similar to an ordinary sparse coding algorithm that is based on updating the coefficients sequentially using a coordinate descent algorithm [START_REF] Friedman | Pathwise coordinate optimization[END_REF], [START_REF] Elad | Why simple shrinkage is still relevant for redundant representations?[END_REF].

C. Atom-Profile Updating Dictionary Learning (APrU-DL)

Recall that the idea of K-SVD to update each atom is indeed a generalization of K-means clustering algorithm. In other words, each atom is updated using only its cluster members, i.e., those signals that have used it in their representations. However, there is a main difference between these two algorithms. As stated in [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF], since in K-means each training data uses only one atom, the atom updating problems are decoupled, while this is not the case in K-SVD. In other words, each training data belongs possibly to multiple atom clusters, thus the atom update problems are indeed coupled. So, we believe that restricting the support of the profiles to be fixed during atom updates, as done in K-SVD, is not well justified.

To overcome this problem, we propose to allow the support of each profile to change during the atom updating procedure. In this way, each atom adaptively finds its cluster members. This idea is used in OS-DL algorithm described in the previous section. In this section, we describe the two-stage counterpart of this algorithm by performing the sparse approximation stage, too. Contrary to OS-DL, by this method, the cardinality of the representations can be controlled. Moreover, as will be seen in the simulations, the new algorithm has a better Algorithm 4 Batch-FISTA(Y,D,X (0) ,λ)

1: Require: Y, D, X (0) ∈ R K×L , λ 2: Initialization: Z (0) = X (0) , c ≤ 1/(2λ max (D T D)), A = D T D, B = D T Y, θ 0 = 1 3: for k = 0, 1, 2 . . . do 4: X (k+1) = S c.λ (Z (k) -cAZ (k) + cB) 5: θ k+1 = (1 + √ 1 + 4θ 2 k )/2 6: β k = (θ k -1)/θ k+1 7:
Z (k+1) = X (k+1) + β k (X (k+1) -X (k) ) 8: end for performance, both in convergence rate and quality of the results.

For the sparse approximation stage any sparse coding algorithm can be used. However, the relaxation-based methods have a better performance compared to greedy ones (to see this, we conduct an experiment in the next section). Among the relaxation-based methods, the soft-thresholding ones are probably more attractive, for both their simplicity and good performance. These methods target the following problem

min D∈D,X 1 2 ∥Y -DX∥ 2 F + λ s ∥X∥ 1 , (13)
where λ s is a regularization constant. Here, we use the Fast Iterative Shrinkage-Thresholding (FISTA) algorithm [START_REF] Beck | Fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] whose global convergence rate is much better than the ordinary IST algorithms, while preserving their simplicity. The batchmode version of FISTA is shown in Algorithm 4, in which, λ max (X) denotes the largest eigenvalue of X. For initialization of the coefficient matrix in each alternation, we use its final estimate at the previous alternation. Once the support of each column of the coefficient matrix is found, we project the associated training signal onto the subspace spanned by the corresponding atoms of the dictionary. This process is called debiasing, and it is known to improve the results [START_REF] Figueiredo | projection for sparse reconstruction: Application to compressed sensing and other inverse problems[END_REF].

In order to update the dictionary in the second stage, we follow the same approach as in OS-DL algorithm. In other words, the problem of updating each atom together with its profile amounts to solve [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF]. Algorithm 5 gives a description of APrU-DL algorithm.

IV. SIMULATIONS

We evaluate the efficiency of our proposed algorithms with two sets of experiments. The first set of experiments is on synthetic data, where we aim to evaluate the capability of our algorithms in recovery of a known dictionary. To this aim, we generate a set of training signals, each as a linear combination using a different set of atoms from an underlying dictionary. We then give these training signals (after adding a certain amount of Gaussian noise) to each algorithm and compare the output dictionary to the original one. In this way, a dictionary learning algorithm should have the ability to extract the common features of the set of signals, which are actually the generative atoms. The second set of experiments is on an autoregressive (AR) signal, where there is no underlying dictionary, and we just evaluate the capability of the algorithms Algorithm 5 APrU-DL 1: Task: Learning a dictionary for Y 2: Initialization: D = D (0) and X = X (0) 3: Repeat: 4: Sparse Approximation: X =Batch-FISTA(Y, D, X, λ s ) 5: Dictionary Update: set E = Y -DX 6: for a = 1, . . . , A do 7:

for i = 1, . . . , K do

8:

E i = E + d i x T [i]
9:

x

[i] = S λ (E T i d i )
10:

d i = E i x [i]
11:

d i = d i /∥d i ∥ 2 12: E = E i -d i x T [i]
13:

end for 14: end for 15: Until convergence in learning a good (i.e., sparsifying) dictionary, or extracting a set of good features. We consider an AR(1) signal as in [START_REF] Skretting | Recursive least squares dictionary learning algorithm[END_REF], and generate the training signals by chopping this signal into a number of blocks. We compare the performance of our proposed algorithms to those of K-SVD and AK-SVD 2,3 . As it was said in Section III-A, for PAU-DL: A = 3, B = 1, and for AK-SVD: A = 1, B = 3.

Our simulations were performed in MATLAB R2010b environment on a system with 3.8 GHz CPU and 8 GB RAM, under Microsoft Windows 7 operating system. As a rough measure of complexity, we will mention the run times of the algorithms.

A. Synthetic Data

We generated a dictionary by normalizing a random matrix of size 20 × 50, with zero-mean and unit-variance independent and identically distributed (i.i.d.) Gaussian entries. A collection of 2000 training signals {y i } 2000 i=1 were produced, each as a linear combination of s different columns of the dictionary, with zero-mean and unit-variance i.i.d. Gaussian coefficients in uniformly random and independent positions. We varied s from 3 to 6. We then added white Gaussian noise with Signal to Noise Ratio (SNR) levels of 10, 20, 30, and 100 dB. The exact value of s was given to PAU-DL, K-SVD, AK-SVD, and APrU-DL (with OMP). For OS-DL we used a fixed value of λ = 0.3. For APrU-DL (with FISTA) we used λ = 0.3 and λ s = 0.6. We applied all algorithms onto these noisy training signals, and compared the resulting recovered dictionaries to the generating one as follows. Assume that d i is a generating atom and di is the atom in the recovered dictionary that best matches d i among the others. We say that the recovery is successful if |d T i di | is above 0.99 [START_REF] Aharon | K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF]. The percentage of the correct recovery was used as the measure of successfully reconstructing the generating dictionary. We performed 100 alternations between sparse approximation and dictionary update stages for all algorithms. The initial dictionary was made by randomly choosing different columns of the training set followed by a normalization. We repeated each experiment (corresponding to a certain value of s and a certain noise level) 50 times and reported the averaged results.

The average percentage of successfully recovering the underlying atoms is shown in Table I. The average execution times of the algorithms for this experiment is shown in Fig. 2. To see the convergence behaviour of K-SVD, AK-SVD, PAU-DL, and APrU-DL, the improvement of recovery ratio along the alternation number is shown in Fig. 3. As we saw in our simulations, the results of AK-SVD were very close to those of K-SVD (as can be seen from Fig. 3). So, AK-SVD has been omitted from Table I. Also, the average execution time of AK-SVD is nearly the same as PAU-DL. The convergence behaviour of OS-DL and APrU-DL with OMP are shown in Fig. 4, where the results of APrU-DL with FISTA are also shown for comparison. With these results in mind, we deduce the following observations:

• PAU-DL has a better successful recovery results compared to both AK-SVD and K-SVD in average. This is especially observable at s = 5 and s = 6. The average execution time of PAU-DL is also much smaller than that of K-SVD. • APrU-DL has the best results in average. Also, the results of APrU-DL with FISTA are better than those with OMP. However, the average runtime of FISTA is higher than that of OMP4 . • OS-DL outperforms PAU-DL and K-SVD, both in convergence rate and the final success rate (this is especially observable at low SNRs and s = 5 and s = 6). This shows the promising performance of one-stage dictionary learning.

• The convergence rate of PAU-DL is better than K-SVD and AK-SVD, while that of APrU-DL is the best. The convergence behaviour of APrU-DL and OS-DL is approximately the same for different values of s, while those of the other three algorithms deteriorate by increasing s.

B. AR(1) signal

In this experiment, we consider an AR(1) signal (according to [START_REF] Skretting | Recursive least squares dictionary learning algorithm[END_REF]), which is generated as v(k) = 0.95v(k - training vector in the sparse approximation stage of PAU-DL, AK-SVD, K-SVD, and APrU-DL. In order to compute SNR in OS-DL, after each update of the dictionary the sparse approximations of the signals have been computed using OMP. For both APrU-DL and OS-DL a value of λ = 0.05 found to yield promising results. For all algorithms 100 alternations were done.

SNR versus alternation number (averaged over 50 trials) is plotted in Fig. 5. Again, the result of AK-SVD was very similar to that of K-SVD, and thus has been omitted from the figure. This figure shows that the SNR is improved during the learning process for all algorithms. The final values of SNR as well as the average execution times of the algorithms are reported in Table II. Based on these results, we deduce the following observations:

• PAU-DL has reached a higher SNR value compared to K-SVD, while its average execution time is much less than that of K-SVD. • OS-DL outperforms PAU-DL. This observation again indicates the promising performance of one-stage dictionary learning algorithms. • APrU-DL has the best results, in the sense of the final value of SNR and the rate of convergence.

As a conclusion, taking into account the performance and the execution time, PAU-DL seems to be the best compromise.

V. CONCLUSION

In this paper, we addressed the dictionary update stage in a general dictionary learning problem. We especially considered the atom-by-atom dictionary update procedure. We introduced the idea of parallel atom-updating and based on this, we derived PAU-DL, an efficient alternative for the well-known K-SVD algorithm which provides atom sequential updating, i.e., fully updating of one atom before moving to the next one. Conversely, the main idea of parallel atom-updating is to update the atoms in parallel. In this way, before moving to update each atom, previous atoms have been partially updated. So, we have more reliable updates of the previous atoms to be used for updating the next atom. We then proposed a novel dictionary learning algorithm which we called OS-DL. This algorithm sequentially updates each atom along with the whole entries of its corresponding row vector in the coefficient matrix. OS-DL differs from the existing dictionary learning algorithms in the sense that it does not perform the first stage, i.e., the sparse approximation stage. This is because, here, we constrain the sparsity criterion on the rows of the coefficient matrix, not on its columns. To have a better control on the sparsity level of the representations, we then proposed the two-stage counterpart of OS-DL, which is obtained by adding the sparse approximation stage to OS-DL. In this way, as it was seen in Section IV, the convergence rate of OS-DL is increased. Our simulations on recovery of a known dictionary, as well as designing a sparsifying dictionary for an AR(1)

Fig. 1 .

 1 Fig. 1. A dictionary learning problem is to factorize a training data matrix Y as Y ≃ DX with X a sparse-column matrix.

Algorithm 1

 1 Parallel Atom-Updating 1: E = Y -DX 2: for a = 1, . . . , A do 3:

2 1 ) 1 :

 211 AK-SVD (A = 1, B = 3) and PAU-DL (A = 3, B = Task: Learning a dictionary for Y 2: Initialization: D = D (0) 3: Repeat: 4: Sparse Approximation: X = Batch-OMP(Y, D, τ ) 5: Dictionary Update: set E = Y -DX 6: for a = 1, . . . , A do 7:

Algorithm 3 2 : 3 :

 323 OS-DL 1: Task: Learning a dictionary for Y Initialization: D = D (0) and X = X (0) The main loop: Set E = Y -DX 4: for a = 1, . . . , A do 5:

  1) + e(k), where e(k) is a zero-mean and unit-variance Gaussian noise. A collection of L = 2000 training signals were made by chopping this signal into vectors of length n = 20. Number of atoms was set to m = 40. As in [32], we computed SNR as SNR = 10 log ∥Y∥ 2 F ⧸∥Y -DX∥ 2 F . For the sparse approximation stage of APrU-DL we have used OMP to have an exact control on the cardinality of the representations. A number of s = 5 atoms were used to approximate each

Fig. 2 .

 2 Fig. 2. Average execution times (in second) versus number of nonzero coefficients.

Fig. 3 .

 3 Fig. 3. Convergence behaviour of K-SVD, AK-SVD, PAU-DL, and APrU-DL (with FISTA) in reconstruction of a known dictionary. Each figure corresponds to a certain amount of sparsity (number of nonzero coefficients). Noise level is SNR = 30 dB.

Fig. 4 .

 4 Fig. 4. Convergence behaviour of OS-DL and APrU-DL with (Batch) OMP and (Batch) FISTA as the sparse approximation stage. The setup is the same as in Fig. 3.

Fig. 5 .

 5 Fig. 5. SNR in dB is plotted versus the alternation number during the learning process.

TABLE II FINAL

 II SNR IN DB AND AVERAGE EXECUTION TIMES (AET) IN SECONDS, FOR VARIOUS ALGORITHMS.

	Algorithm	K-SVD	PAU-DL	OS-DL	APrU-DL
	SNR(dB)	17.03	17.16	17.59	17.67
	AET(s)	14.57	2.53	12.93	21

signal indicate that our algorithms outperform K-SVD.

Note that this differs from the usual sparse approximation stage performed in a typical dictionary learning algorithm. Actually, in[START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF], the sparsity constraint (i.e., the ℓ 1 norm penalty) is on the rows of the coefficient matrix, not on its columns. As explained in the text, we include this constraint to prevent each profile to be filled during the atom updating procedure.

As pointed out in the previous works, e.g.,[START_REF] Skretting | Recursive least squares dictionary learning algorithm[END_REF], the performances of MOD and K-SVD are very similar in these experiments. So, we omitted MOD from the simulations.

For K-SVD, AK-SVD and OMP we have used K-SVD-Box v10 and OMP-Box v10 available at http://www.cs.technion.ac.il/ ∼ ronrubin/software. html

Note that the average execution times for all algorithms have been calculated for 100 alternations. As can be seen from Fig.3, APrU-DL has converged in about 20 alternations. So, considering the convergence times, APrU-DL is fast enough compared to K-SVD.

Foundation (INSF) under Contract 91004600, and by European project 2012-ERC-AdG-320684 CHESS.