Déborah Nourrit-Lucas

Adaté Olivier Tossa

Grégory Zélic

Didier Delignières

Learning, Motor Skill and Long-Range Correlations

Keywords: Motor skills, motor learning, 1/f noise, long-range correlations, degeneracy

Long-range correlations have been evidenced in a number of experiments, generally using over-learned and over-practiced tasks. We hypothesized that long-range correlation could represent the byproduct of learning. We analyzed the series of periods produced by a group of expert and a group of novices during prolonged trials on a ski-simulator. Results showed a very low variability in expert's series, as compared to novices. Fractal analyses showed that fluctuations were significantly more structured and correlated in experts. These results suggest that learning could be conceived as the progressive installation of complexity in the system.

Introduction

Fractal fluctuations have been recently evidenced in series of performances collected in cyclical or repetitive tasks, such as serial reaction time [START_REF] Van Orden | Self-organization of cognitive performance[END_REF], finger tapping [START_REF] Gilden | 1/F Noise in Human Cognition[END_REF][START_REF] Lemoine | Testing for the presence of 1/f noise in continuation tapping data[END_REF], circle drawing [START_REF] Torre | Long-range correlation properties in motor timing are individual and task specific[END_REF], forearm oscillations [START_REF] Torre | Oscillating in Synchrony with a Metronome: Serial Dependence, Limit Cycle Dynamics, and Modeling[END_REF], reciprocal aiming [START_REF] Slifkin | Amplitude requirements, visual information, and the spatial structure of movement[END_REF], bimanual coordination [7], walking, running [START_REF] Hausdorff | Altered fractal dynamics of gait: Reduced stride-interval correlations with aging and Huntington's disease[END_REF][START_REF] Jordan | Long range correlations in the stride interval of running[END_REF], or in bimanual coordination tasks [START_REF] Torre | 1/f β fluctuations in bimanual coordination: An additional challenge for modelling[END_REF].

These fluctuations are characterized by specific properties, namely self-similarity, or scale invariance, meaning that statistical features in the series are similar whatever the scale of observation, and long-range correlation, revealed by the presence of positive serial correlations between successive values, which persist over time, often over dozens, sometimes over hundreds of observations [START_REF] Diniz | Contemporary theories of 1/f noise in motor control[END_REF][START_REF] Eke | Physiological time series: distinguishing fractal noises from motions[END_REF]. These ubiquitous and amazing statistical properties present a special interest for behavioral scientists, as they are conceived as theoretically closely linked to complexity, adaptability and health [START_REF] Goldberger | Fractal dynamics in physiology: Alterations with disease and aging[END_REF]. Long-range correlated series represent the typical output of complex and healthy organisms, characterized by essential properties of robustness and adaptability. In contrast, aging and disease seem marked by a loss of complexity, which is typically revealed by a decrease of long-range correlations in output series [START_REF] Hausdorff | Altered fractal dynamics of gait: Reduced stride-interval correlations with aging and Huntington's disease[END_REF]. These close relationships between long-range correlations, robustness and adaptability are of special importance here, as the two latter represent the main properties of the skilled behavior, and the essential by-products of learning.

Importantly, long-range correlations have been essentially evidenced in overlearned tasks, such as tapping, circle drawing, reciprocal aiming, walking or running. In all cases these tasks required the exercise of basic skills, acquired from years and extensively practiced.

Obviously, seeking for long-range correlation in performance series requires the collection of very long time series of hundreds successive data points, supposing that participants are sufficiently familiar with the task at hand. On could hypothesize, however, that the presence of such fractal fluctuation could be related to the fact that performance is underlain by a well established skill.

An interesting result supporting the present claim has been reported by Wijnants, Bosman, Hasselman, Cox and Van Orden [START_REF] Wijnants | 1/f Scaling in Movement Time Changes with Practice in Precision Aiming[END_REF]. The authors analyzed serial correlations in series of movement times in a reciprocal aiming task. The task was performed with the non-dominant hand, and the experimental design included five successive blocks of 1100 trials. Results showed an increase of serial correlations in the series with practice, with a clearer evidence for 1/f fluctuation in the last block.

It could be interesting here to clearly distinguish between learning and practice. Learning can be defined as the acquisition of a new skill, which is not initially present in the repertoire of the individual [START_REF] Nourrit | On discontinuities in motor learning: A longitudinal study of complex skill acquisition on a ski-simulator[END_REF][START_REF] Teulier | The nature of the transition between novice and skilled coordination during learning to swing[END_REF][START_REF] Teulier | The evolution of oscillatory behavior during learning on a ski simulator[END_REF]. In contrast, practice refers to the repeated exercise of a task, leading to a refinement of an existing skill, but not necessarily to the adoption of a qualitatively modified behavior. Practice is essential for learning, but extensive practice is often required for an effective learning to occur [START_REF] Nourrit | On discontinuities in motor learning: A longitudinal study of complex skill acquisition on a ski-simulator[END_REF]. From this point of view, the aforementioned results seem more related to the effects of practice and to those of effective learning.

Practice and learning, however, often produce similar and related effects, such as the decrease of performance variability, an enhancement of efficiency, a better robustness facing external perturbations, and a better adaptation to related tasks [START_REF] Schmidt | Motor Control And Learning: A Behavioral Emphasis[END_REF]. As such, long-range correlation could be conceived a logical byproduct of practice and learning.

In the present experiment, we analyzed performance series collected in novice and expert participants in a complex task. We hypothesized to evidence stronger long-range correlations in experts, suggesting that the expert behavior was characterized by a higher level of complexity than the initial, novice behavior.

Material and Methods

Participants

Nine volunteered participants (two females, seven males) participated in this study. There were separated in two experimental groups. The expert group included one female and three males (mean age: 39.2 years ±6.3, mean weight: 73.2 kg ± 8.46; mean height: 179.6 cm ± 3.5). These participants were involved in a series of experiments on the ski-simulator task.

Ten years ago, they were involved in a first longitudinal experiment including 390 1-min trials over 13 weeks, and were proven to have adopted a skilled behavior, qualitatively different than their initial behavior on the task [START_REF] Nourrit | On discontinuities in motor learning: A longitudinal study of complex skill acquisition on a ski-simulator[END_REF]. They were also involved in two retention tests, the first one five months after the completion of the learning protocol [START_REF] Deschamps | Influence of a stressing constraint on stiffness and damping functions of a ski simulator's platform motion[END_REF], and the second ten years after [START_REF] Nourrit-Lucas | Persistent coordination patterns in a complex task after 10 years delay -How validate the old saying "Once you have learned how to ride a bicycle, you never forget![END_REF]. In both cases the retention tests evidenced the persistence of the skilled behavior initially acquired by participants.

The novice group (one female, four males, mean age: 23.2 years ± 2.5, mean weight: 70.5 kg ± 4.2; mean height: 1.80 cm ± 5.8) was composed of occasional skiers (with average three days of practice per year), but none had specific training on the ski-simulator. All participants signed a consent form, and were not paid for their participation.

Experimental device

The task was performed on a ski-simulator (Skier's Edge Co., Park City, UT), which consisted of a platform on wheels, which moved back and forth on two bowed, parallel metal rails (Figure 1). We used a modified version of the simulator by replacing the two independent feet supports of the original apparatus with a 30-cm wide board, in unstable balance over a sagittal rotation axis (for more details see [START_REF] Nourrit | On discontinuities in motor learning: A longitudinal study of complex skill acquisition on a ski-simulator[END_REF]) -Insert Figure 1 -

Procedure

Participants were instructed to make cyclical sideways movements on the ski simulator, "as ample and frequent as possible". They had to keep their hands behind their back at all times, and to fix their eyes on a point located on the floor, three meters in front of the apparatus.

They performed a unique session of 10 minutes, allowing the performance of approximately 550 complete oscillations on the apparatus.

Data collection

A passive marker was fixed in the front of the simulator platform. The displacement of this marker was recorded in three dimensions by a VICON motion analyzer (Biometrics) with seven cameras (100 Hz). Analyses focused on the series of positions of the platform, along the transverse axis, computed from the collected 3-D data.

The position time series were filtered with a dual-pass Butterworth filter with a cut-off frequency of 10 Hz. A peak-finding algorithm was used to localize the left reversal points of the platform motion and the period was calculated for each oscillation as the time interval between two successive reversal points. We retained for analysis the 512 last points of the series, for each participant.

Statistical analysis

We first characterized series in terms of descriptive statistics (mean and standard deviation).

We then applied three analyses aiming at evidencing and measuring serial correlation in the series.

Autocorrelation function

Autocorrelation functions were computed up to lag 30. We extracted from these function two variables of interest: the lag-one autocorrelation [ACF(1)], and the average autocorrelation for lags comprised between 10 and 30 [<ACF> ]. As long-range correlated series are characterized by the persistence of correlations over time, <ACF>(10-30) was expected to be significantly higher in the expert group.

Detrended Fluctuation Analysis (DFA)

DFA is a widely used method that allows to quantify correlation in time series [START_REF] Peng | Longrange anticorrelations and non-Gaussian behavior of the heartbeat[END_REF]. The series x(t) is first integrated, by computing for each t the accumulated departure from the mean of the whole series:

X(t) = x(i ) -x [ ] i=1 t å (1)
This integrated series of length N is divided into k non-overlapping intervals of length n. The last N -(kn) data points are excluded from analysis. In each interval, a least squares line is fit to the data (representing the trend in the interval). The series X(t) is then locally detrended by substracting the theoretical values X n (t) given by the regression. For a given interval length n, the characteristic size of fluctuation for this integrated and detrended series is calculated by:

F(n) = 1 N -(kn) X(t) -X n (t) [ ] 2 t=1 N-kn å (2)
This computation is repeated over all possible interval lengths. Typically, F increases with interval length n. A power law is expected, as

F(n) µ n a (3)
 is expressed as the slope of the double logarithmic plot of F(n) as a function of n. The value α = 0.5 indicates the absence of correlations (white noise), α > 0.5 indicates persistent longrange correlations, meaning that large (small) values are more likely to be followed by large (small) values.

We considered interval lengths ranging from n = 10 to n = N/2. In order to avoid any bias due to the logarithmic distributions of the points in the diffusion plots, we divided the abscissa into intervals of 0.1(log 10 Δt), and computed the average points within each interval (13 points were obtained for an initial series length of 512 data points). Finally, in order to control for the effects of noisy perturbations that mainly affect short-term fluctuations and tend to flatten the diffusion plot, we focused on the long-term slope (i.e. the 6 last points [START_REF] Box | Time Series Analysis: Forecasting and Control[END_REF].

Power spectral density analysis (PSD)

This method works on the basis of the periodogram obtained by the Fast Fourier Transform algorithm. In the frequency domain, long-range correlated series are characterized by the following scaling law:

S(f)  1/f β ( 4
)
where f is the frequency and S(f) the correspondent squared amplitude. β is estimated by calculating the negative slope (-β) of the line relating log(S(f)) to log f.

We also used the improved version of PSD proposed by Eke et al. (2000) [START_REF] Eke | Physiological time series: distinguishing fractal noises from motions[END_REF], which uses a combination of preprocessing operations: First the mean of the series is subtracted from each value, and then a parabolic window is applied: each value in the series is multiplied by the following function:

)² 1 1 2 ( 1 ) (     N j j W for j = 1, 2, …, N.
(

Thirdly a bridge detrending is performed by subtracting from the data the line connecting the first and last point of the series. Finally the fitting of β excludes the high-frequency power estimates (f > 1/8 of maximal frequency). This method was proven to provide more reliable estimates of the spectral index β, and was designated as low PSD we .

Group comparisons

Considering the low sample sizes and the strong inhomogeneity of variances, we used nonparametric Mann-Whitney U tests for comparing central tendencies between groups. The significance threshold was set at 0.05.

Results

Descriptive statistics

We present in Figure 2 two example series obtained with a novice (top panel) and an expert (bottom panel). The samples of mean periods were as follows: Novices : {0.87, 1.13, 1.02, 0.96, 0.81} ; Experts : {0.82, 0.93, 0.87, 0.84}. There was no difference between the two groups (Experts: 0.87 sec ± 0.04; Novices: 0.96 sec ± 0.13; U = 6; Z = 0.98; p = 0.327; exact p = 0.413).

This figure, however, suggests evident differences in terms of variance. Indeed, the samples of standard deviations were the following: Novices : {0.17, 0.24, 0.39, 0.10, 0.25} ; Experts : {0.02, 0.04, 0.04, 0.03}. There was a statistical difference between the two groups (Experts: 0.03 sec ± 0.01; Novices: 0.23 sec ± 0.11; U = 0; Z = -2.45; p = 0.014; exact p = 0.016). As expected, experts performed the task with a very low variability, as compared with novices.

-Insert Figure 2 -

Autocorrelation functions

We present in Figure 3 the point-by-point average autocorrelation functions for the two groups. There were evident graphical differences between these two average functions. In the novice group, the autocorrelation function present just significant values for the three first lags, and then reaches very quickly values close to zero. This kind of auto-correlation function is typical of short-range correlated processes [START_REF] Box | Time Series Analysis: Forecasting and Control[END_REF]. In contrast, the average autocorrelation function of the expert group present presents a very slow decay, with significant values up to the 30 th lag. This kind of auto-correlation function corresponds to those obtained with longrange correlated series.

The following values were observed for ACF(1): Novices : {0.38, 0.03, 0.28, 0.11, 0.07} ; Experts : {0.40, 0.27, 0.38, 0.23}. There was no difference between the two groups, however, because of the large variability in the novice group (Experts: 0.32 ± 0.08; Novices: 0.17 ± 0.15; U = 4; Z = -1.47; p = 0.142; exact p = 0.190).

<ACF>(10-30) values were the following: Novices : {0.01, 0.02, 0.03, 0.00, 0.00} ; Experts : {0.31, 0.09, 0.24, 0.04}. There was a significant difference between the two groups (Experts: 0.17 ± 0.12; Novices: 0.01 ± 0.01; U = 0; Z = -2.45; p = 0.014; exact p = 0.016).

-Insert Figure 3 -

Detrended Fluctuation Analysis

We present in Figure 4 the point-by-point average diffusion plots, for the two groups. In both case, the diffusion plots present a global linear shape. A clear flattening appears for the Expert group, suggesting the influence of a white noise component in the series. This influence is less apparent for the Novice group, essentially because the global shape is close to that expected for white noise processes (α = 0.5).

The individual values were the following: Novices : {0.69, 0.63, 0.65, 0.53, 0.45} ; Experts : {1.30, 0.69, 1.16, 1.08}. There was a significant difference between the two groups (Expert:

1.06 ± 0.26; Novices: 0.59 ± 0.10; U = 0; Z = -2.45; p = 0.014; exact p = 0.016).

-Insert Figure 4 -

Power Spectral Density analysis

We present in Figure 5 the point-by-point average bi-logarithmic power spectra, for the two groups. The flattening of the spectra in the high g-frequency region revealed for both group the influence of a white noise component. The individual values of the β exponent, computed over the low-frequency region, were: Novices : {-0.16, 0.24, 0.74, 0.14, 0.17}; Experts :

{1.57, 1.04, 1.55, 1.14}. There was a significant difference between groups (Experts: 1.32 ± 0.28; Novices: 0.23 ± 0.33; Z = -2.45; p = 0.014; exact p = 0.016).

-Insert Figure 5 -

Discussion

Motor learning has been classically assumed to be characterized by the selection of the most efficient behavioral solutions, a decrease of performance variability, and an increase of smoothness in movement trajectories [START_REF] Schmidt | Motor Control And Learning: A Behavioral Emphasis[END_REF]. This point of view tends to induce the idea that learning yields a kind of simplification of the system, through the selection of proper procedures and the elimination of errors.

The present results confirm these classical assumptions, and especially the very low variability of cyclical performance in experts. The most important result, however, is the increase of serial correlations in experts, with regards to the levels observed in novices.

Expert performance seems characterized by a more complex and structured dynamics that that of novices.

This result could be interestingly related to a recent work that linked long-range correlation and degeneracy [START_REF] Delignières | Degeneracy and long-range correlations[END_REF]. Degeneracy is a design principle of complex systems which has been proposed for explaining the co-existence of the a priori paradoxical properties of robustness and evolvability [START_REF] Whitacre | Degeneracy: a link between evolvability, robustness and complexity in biological systems[END_REF]. Robustness refers to the capacity to maintain a function despite internal or external perturbations, and evolvability to the capacity to adapt to perturbations by adopting new behavior and functions. The concept of degeneracy refers to a partial overlap in the functions of the multiple components within the system. In degenerate systems, structurally different components can perform similar functions under certain conditions, but can also assume distinct roles in others conditions [START_REF] Edelman | Degeneracy and complexity in biological systems[END_REF][START_REF] Whitacre | Degeneracy: a design principle for achieving robustness and evolvability[END_REF].

Delignières and Marmelat [START_REF] Delignières | Degeneracy and long-range correlations[END_REF] propose a model of degenerate neural network composed of a chain of partially overlapping pathways. They manipulated degeneracy through the number of alternative pathways in the model. A simulation study showed that (1) such a degenerate model produces long-range correlated series, (2) the strength of correlations in the output depends on the level of degeneracy in the model, and (3) a minimal threshold in degeneracy is necessary for producing long-range correlations.

The present experiment suggests that learning could be understood as the progressive installation of degeneracy in the system. Learning is not the selection of the most appropriate solution, but the coordination of a complex network composed of multiple, alternative and overlapping pathways for producing a given outcome. Learning can then be conceived as an increase in complexity of the neural networks that underlie performance, and the overlapping between alternative pathways explains the presence of long-range correlations in output series. This enrichment of neural networks could explain the property of robustness of motor skills, essentially revealed in retention tests, but also the properties of generalizability and transfer, which are considered essential for the completeness of learning [START_REF] Schmidt | Motor Control And Learning: A Behavioral Emphasis[END_REF].

Schöllhorn, Hegen, and Davids [START_REF] Schöllhorn | The Nonlinear Nature of Learning -A Differential Learning Approach[END_REF] recently developed innovative ideas about learning that could be in resonance with the previous finding. They proposed a so-called differential learning approach, that explicitly aimed to exploit the system's complexity by its confrontation to complex and changeable environments and constraints. It is noteworthy to note, however, that this enrichment in complexity also occurs during the practice of very close and simple tasks, such as the reciprocal aiming task used by Wijnants and coll. [START_REF] Nourrit-Lucas | Persistent coordination patterns in a complex task after 10 years delay -How validate the old saying "Once you have learned how to ride a bicycle, you never forget![END_REF].

Figure 1 :Figure 2 :

 12 Figure 1: The ski simulator

Figure 3 :

 3 Figure 3: Mean autocorrelation functions. Top: Novice group (N = 5); Bottom: Expert group

Figure 4 :

 4 Figure 4: Mean DFA diffusion plots for the Novice (grey) and Expert (white) groups. Dashed

Figure 5 :

 5 Figure 5: Mean Log-log power spectra for the Novice (top) and Expert (bottom) groups.

  Figure 4


Mugurel Ionut ¸andreica
email: mugurel.andreica@cs.pub.ro

Practical Algorithmic Optimizations for Finding Maximal Matchings in Induced Subgraphs of Grids and Minimum Cost Perfect Matchings in Bipartite Graphs

Keywords: minimum cost perfect matching, maximal matching, maximum flow, grid graph. 2010 MSC: 05C21

In this paper we present practical algorithmic optimizations addressing two problems. The first one is concerned with computing a maximal matching in an induced subgraph of a grid graph. For this problem we present a faster sequential algorithm using bit operations and a way of implementing it in a parallel environment. The second problem is concerned with computing minimum cost perfect matchings in bipartite graphs. For this problem we extend the idea behind the Hopcroft-Karp maximum matching algorithm and then we consider a more general situation in which multiple minimum cost perfect matchings need to be computed in the same graph, under certain cost restrictions. We present experimental results for all the proposed optimizations.

Introduction

The problem of computing maximum or maximal matchings in bipartite graphs has been considered many times in the scientific literature. Many of the proposed algorithms use the fact that computing a maximum matching in a bipartite graph is equivalent to computing a maximum flow in a slightly modified graph. Thus, results from the theory of network flows can be used for computing maximum matchings. If only a maximal matching is needed, then simpler greedy-type algorithms can be employed. In this paper we present several practical algorithmic improvements for some of the algorithms used for computing maximal matchings in grid graphs and minimum cost perfect matchings in bipartite graphs.

The rest of this paper is structured as follows. In Section 2 we define the main terms and techniques used in this paper. In Section 3 we discuss related work. In Section 4 we present faster algorithms for computing maximal matchings in induced subgraphs of grid graphs based on algorithms which use bit operations. In Section 5 we extend an idea used for computing maximum matchings in bipartite graphs to the computation of minimum cost perfect matchings. The idea consists of using multiple edge-disjoint augmenting paths per iteration in order to reduce the number of iterations of the algorithm. In Section 6 we consider another perfect matching problem. In this problem we are interested in computing a minimum cost perfect matching in a complete bipartite graph under certain restrictions regarding the cost computation. The cost of the matching is considered to be equal to the sum of the costs of the edges from the matching except for the cost of the minimum cost edge from the matching (i.e. the minimum cost edge of the matching is considered to have cost 0 when computing the cost of the matching). In Section 7 we present experimental evaluations of all the algorithms discussed in this paper. Finally, in Section 8 we conclude.

Terms and Definitions

A bipartite graph is a graph whose vertices can be split into two sets L (left) and R (right). We consider the vertices to be numbered from 1 to |L| in the left set and from 1 to |R| in the right set (it is acceptable to have vertices with the same number in the graph, because they will be differentiated based on the set L or R to which they belong). Every edge (x, y) of the graph is between a node x ∈ L and a node y ∈ R. A matching in a bipartite graph is a set of edges such that no two edges in the set have a common vertex. A maximum matching is a matching of maximum cardinality. A maximal matching is a matching to which no more edges can be added (i.e. all the edges outside of the matching have a common vertex with at least one edge from the matching). A perfect matching is a matching in which every node of the graph is an endpoint of an edge from the matching (such a matching may exist only when |L| = |R|).

In order to reduce the maximum matching problem to a maximum flow problem we need to construct a directed graph as follows. We will have a special node S called the source and another special node T called the sink. We will also keep all the nodes from the given bipartite graph. Each edge (x, y) of the original bipartite graph will be replaced by a directed arc from x to y having capacity 1. We will also add capacity 1 arcs from S to every node x ∈ L and from every node y ∈ R to T . In case the edges of the bipartite graph have costs these costs are maintained on the directed arcs from the nodes x ∈ L to the nodes y ∈ R (we will denote by c(x, y) the cost of the edge between x ∈ L and y ∈ R). All the arcs having S or T as an endpoint will have cost 0.

One of the best known maximum flow algorithms is the Edmonds-Karp algorithm [START_REF] Edmonds | Theoretical improvements in algorithmic efficiency for network flow problems[END_REF]. This algorithm can be summarized as follows: As long as possible find a shortest path from S to T in the residual graph and augment the flow along that path. When arc costs are involved the algorithm can be adjusted in order to find a minimum cost path from S to T in the residual graph. Note that the residual graph may contain negative costs. This version of the Edmonds-Karp algorithm is known as the successive shortest path algorithm [START_REF] Todinov | Flow Networks: Analysis and Optimization of Repairable Flow Networks, Networks with Disturbed Flows, Static Flow Networks and Reliability Networks[END_REF]. A simple breadth-first search algorithm is used for finding a shortest path in the first case (i.e. when edge costs are not involved), while a minimum cost path algorithm needs to be used in the second case (i.e. when edge costs are involved), for instance, Bellman-Ford-Moore [START_REF] Papaefthymiou | Implementing parallel shortest-paths algorithms[END_REF] or even Dijkstra's algorithm [START_REF] Todinov | Flow Networks: Analysis and Optimization of Repairable Flow Networks, Networks with Disturbed Flows, Static Flow Networks and Reliability Networks[END_REF] after modifying the graph's arc costs in order to remove negative costs. Thus, the algorithm consists of multiple iterations, in each of which the flow is increased along a single path. The most time consuming part in each iteration is the traversal of the graph in order to find an augmenting path. In a graph with V vertices and E arcs finding the shortest augmenting path takes O(V + E) time when no costs are involved and O(V • E) time when costs are involved (or O(V + E • log(V)) or O(E + V • log(V)) time when Dijkstra's algorithm is used on the modified residual graph costs). Then, augmenting the flow along the found path is easy (it takes only O(V) time). In the case of bipartite graphs it is sufficient to find a path from S to an unmatched vertex in R (because this vertex is directly connected to T through an existing arc in the residual graph).

Related Work

The best algorithm for computing a maximum matching in sparse bipartite graphs is the Hopcroft-Karp algorithm [START_REF] Hopcroft | An n5/2 algorithm for maximum matchings in bipartite graphs[END_REF], which has a time complexity of O(E • √

V) where V is the number of vertices and E is the number of edges of the graph. For dense bipartite graphs the algorithm proposed in [START_REF] Alt | Computing a maximum cardinality matching in a bipartite graph in time O(n 1.5 m logn )[END_REF] has a slightly better time complexity of

O(V 1.5 E log(V)
). Both of these algorithms have a better time complexity than the Edmonds-Karp algorithm for finding a maximum flow presented in the previous section. However, due to its simplicity, the Edmonds-Karp algorithm is used in many practical implementations. Moreover, experimental evaluations showed that for some types of bipartite graphs some modified versions of the Edmonds-Karp algorithm (which use breadth-first search from all the source's neighbors for finding augmenting paths) are faster than the Hopcroft-Karp algorithm, despite having a worse theoretical time complexity [START_REF] Cherkassky | Augment or push? a computational study of bipartite matching and unit capacity maximum flow algorithms[END_REF].

Edmonds-Karp is not the only algorithm for computing maximum flows in graphs. In fact, many such algorithms were proposed in the scientific literature. Some of the most popular ones are Dinic's algorithm [START_REF] Dinic | Algorithm for solution of a problem of maximum flow in a network with power estimation[END_REF]), Karzanov's algorithm [START_REF] Karzanov | Determining the maximal flow in a network by the method of preflows[END_REF] and the push-relabel maximum flow algorithm [START_REF] Goldberg | A new approach to the maximum flow problem[END_REF].

A minimum cost perfect bipartite matching can be computed in O(V 3 ) time using the Hungarian algoritm [START_REF] Munkres | Algorithms for the assignment and transportation problems[END_REF]. The successive shortest path algorithm for minimum cost maximum flows can be implemented in O(V • (E + V • log(V))) time in order to compute a minimum cost maximum matching by using Fibonacci heaps [START_REF] Fredman | Fibonacci heaps and their uses in improved network optimization algorithms[END_REF]. The algorithm consists of O(V) iterations and each iteration runs in O(E + V • log(V)) time. Dynamic versions of the minimum cost perfect bipartite matching problem, in which edge costs can be changed, have also been considered [START_REF] Mills-Tettey | The dynamic hungarian algorithm for the assignment problem with changing costs[END_REF].

Maximum matchings can also be computed in general graphs, not just bipartite graphs (see, for instance, Gabow's algorithm [START_REF] Gabow | An efficient implementation of edmonds' algorithm for maximum matching on graphs[END_REF], having an O(V 3 ) time complexity). Minimum cost perfect matchings have also been considered in some special classes of graphs, e.g. graphs induced by points in the plane [START_REF] Varadarajan | A divide-and-conquer algorithm for min-cost perfect matching in the plane[END_REF]. Greedy algorithms for maximal matchings, including parallel versions, were presented in [START_REF] Blelloch | Greedy sequential maximal independent set and matching are parallel on average[END_REF]. The problem of maintaining maximal matchings in dynamic graphs has been addressed in [START_REF] Neiman | Simple deterministic algorithms for fully dynamic maximal matching[END_REF].

Faster Algorithm for Maximal Matchings in Induced Subgraphs of Grid Graphs using Bit Operations

We consider an M • N grid graph in which every node has a coordinate (x, y) (0 ≤ x ≤ N -1, 0 ≤ y ≤ M -1) and some nodes are missing. The graph is defined by the implicit adjacency structure of the existing nodes (i.e. two nodes at distance 1 in the grid are neighbors). We are interested in computing a maximal matching in this graph. Note that a maximal matching simply implies that no other edge of the graph can be added to the matching and not that the matching has maximum cardinality. Computing a maximum cardinality matching can be done easily, because the graph is bipartite (we can separate the nodes into two groups based on the parity of their sum of x and y coordinates) and there are many polynomial-time maximum matching algorithms in such graphs (see Section 3).

Computing a maximal matching can be achieved faster, in only O(M • N) time. Let's consider the following Greedy algorithm (Algorithm 1) which traverses the grid graph in increasing order of the y-coordinate and for each y in inceasing order of the x-coordinate.

Algorithm 1 Greedy O(M•N) Algorithm for Finding a Maximal Matching in an Induced Subgraph of a Grid Graph C = 0 {At the end of the algorithm C will be the size of the maximal matching.} for y = 0 to M -1 do for x = 0 to N -1 do if node (x, y) exists in the graph then if y > 0 and node (x, y -1) exists in the graph and node (x, y -1) is not matched then Match the nodes (x, y) and (x, y -1) C = C + 1 else if x > 0 and node (x -1, y) exists in the graph and node (x -1, y) is not matched then Match the nodes (x, y) and (x -1, y)

C = C + 1 end if end if end for end for
We can implement a faster version of the Algorithm 1 by using bit operations. Note that the presented algorithm will only compute the size of the maximal matching and not the matching itself. The speed increase is due to using bit operations and handling multiple nodes at the same time. We will split each row of the grid graph (corresponding to a y-coordinate) into blocks of B bits. Block i of each row contains bits referring to the coordinates i • B, . . . , (i + 1) • B -1. We will denote by block(y, i) the block i of the row corresponding to the coordinate y. We will have bit j of block(y, i) set to 1 if the node (i • B + j, y) exists in the graph, and set to 0 otherwise (0 ≤ j ≤ B -1). We will traverse the graph from y = 0 to y = M -1 as in Algorithm 1. During the traversal we will maintain a row of blocks corresponding to the previous row in which 1 bits will correspond to existing unmatched nodes. When considering a new row y, the first step is to perform an AND between the current row and the previous row of unmatched nodes. All the 1 bits in the result of this operation will represent nodes from the current row which are matched to nodes from the previous row. After performing this match we will clear the matched 1 bits from the current row. The next step is to match nodes from the current row which are adjacent to each other faster than O(N) time. In order to achieve this we will need to use a preprocessing step. For each sequence S of B bits we will compute MCnt(S ) that is the number of pairs of adjacent bits matched in S and MRes(S ) the B-bit sequence containing the remaining unmatched 1 bits of S . We will start with MCnt(S ) = 0 and MRes(S ) = S . Then we will traverse all the bits j of S from 1 to B -1. If MRes(S )( j) = 1 and Mres(S )( j -1) = 1 then we increase MCnt(S ) by 1 and we clear the bits j and j -1 in MRes(S ). Thus, we can compute MCnt(S ) and MRes(S ) in O(B) time, obtaining a preprocessing time of O(2 B • B). Within the same time complexity we will also compute for each B-bit sequence S the number of 1 bits in S , BCnt(S ).

With these values we can perfom the matching on the current row y. We will consider each block i from 0 to (N -1)/B and we will maintain the state of the current row as a sequence of blocks crow. First we copy block(y, i) to crow(i). Then, if i > 0 and bit B -1 of crow(i -1) is 1 and bit 0 of crow(i) is 1 we match these two bits and then we set them to zero. Afterwards we replace crow(i) by MRes(crow(i)). The detailed algorithm is presented in Algorithm 2.

C = 0 prow(i) = 0 (0 ≤ i ≤ (N -1)/B) for y = 0 to M -1 do crow(i) ← block(y, i) (0 ≤ i ≤ (N -1)/B) for i = 0 to (N -1)/B do vmatch(i) = crow(i) AND prow(i) C = C + BCnt(vmatch(i))) crow(i) = crow(i) XOR vmatch(i) end for C = C + MCnt(crow(0)) crow(0) = MRes(crow(0)) for i = 1 to (N -1)/B do if bit B -1 of crow(i -1) is 1 and bit 0 of crow(i) is 1 then C = C + 1 Clear bit B -1 of crow(i -1) and bit 0 of crow(i). end if C = C + MCnt(crow(i)) crow(i) = MRes(crow(i)) end for prow(i) = crow(i) (0 ≤ i ≤ (N -1)/B) end for The time complexity of Algorithm 2 is O(2 B • B + M • N/B). The O(2 B • B)
term is the time complexity of the preprocessing stage and the O(M • N/B) is the time complexity of the actual algorithm. The presented algorithm can even be implemented in a parallel manner. First of all the preprocessing stage is obviously parallelizable: each of the 2 B values of the tables MRes, MCnt and BCnt ca be computed independently. In order to parallelize the actual algorithm we will need to refactor it first. We will first perform all the horizontal matchings on each of the M rows. We can first perform the matching within each block of each row independently in parallel and store the result in a variable specific to each (row, block) pair (this means that we would have such a variable for each block of each row). Then we can handle the matching between bit 0 of odd-numbered blocks and bit B -1 of the preceding even-numbered block in parallel, followed by another stage in which we handle the matching between bit 0 of even-numbered blocks and bit B -1 of the preceding odd-numbered block in parallel. Then we can handle matchings between nodes on different rows. In order to parallelize this stage we will first consider all the rows corresponding to odd y coordinates being matched to the adjacent row with a smaller and even y coordinate. Obviously, each block of all of these M/2 (we consider integer division) rows can be handled independently in parallel. Then we will consider all the rows corresponding to even y coordinates being matched to the adjacent row with a smaller and odd y coordinate. Each block of these M -M/2 rows can also be handled independently, in parallel. The parallel algorithm presented here can use up to 2 B processors in the preprocessing stage and up to M • N/B processors in the maximal matching computation stage. Note that the result of the parallel version may differ from the result of the sequential algorithm (Algorithm 2) because a different maximal matching will be computed (due to the different order of performing the vertical and horizontal matches).

Using a Maximal Set of Edge-Disjoint Paths for Reducing the Number of Iterations of Minimum Cost Perfect Matching Algorithms

The best algorithm for computing a maximum matching in a sparse bipartite graph is the Hopcroft-Karp algorithm [START_REF] Hopcroft | An n5/2 algorithm for maximum matchings in bipartite graphs[END_REF] which has a time complexity of O(E • √ V). The main idea behind that algorithm is to enhance a standard augmenting path algorithm as follows.

After each BFS traversal of the graph in order to find an augmenting path, the matching will not be increased only along one path, but rather along a maximal set of edge-disjoint shortest paths (note that in this case edge-disjoint paths are also vertex-disjoint paths, because they are paths in a shortest path tree; the only common vertex is the source S ).

The same idea can be used when computing a minimum cost perfect matching. At each iteration of the successive shortest path algorithm [START_REF] Todinov | Flow Networks: Analysis and Optimization of Repairable Flow Networks, Networks with Disturbed Flows, Static Flow Networks and Reliability Networks[END_REF] we need to find a minimum cost path in the residual graph. Note that the residual graph may have arcs with negative costs, but does not have negative cycles. Thus, we either need to use a shortest path algorithm which supports negative costs (e.g. Bellman-Ford-Moore [START_REF] Papaefthymiou | Implementing parallel shortest-paths algorithms[END_REF])) or we need to modify the costs in order to obtain non-negative costs only and, thus, use Dijkstra's algorithm [START_REF] Todinov | Flow Networks: Analysis and Optimization of Repairable Flow Networks, Networks with Disturbed Flows, Static Flow Networks and Reliability Networks[END_REF].

No matter what shortest path algorithm we use, at the end of the algorithm we have the minimum cost of a path starting at the source and ending at each node x ∈ R. We can sort all these nodes in ascending order of the cost of the minimum cost path to reach them (ignoring the unreachable nodes, if any). Then, rather than only increasing the matching along the minimum cost augmenting path, we can consider these nodes in sorted order. For each node x we trace back its shortest path to the source. If the matching was already augmented at the current iteration along at least one edge of the path, then we ignore node x and we move on to the next one. If the cur-rent path does not intersect with any of the paths along which the matching was augmented at the current iteration then we can augment the matching along this path and mark its edges in order to know that no other shortest path containing (some of) these edges can be used for augmenting the matching at the current iteration.

A direct implementation of this modified matching augmentation algorithm takes O(V 2 ) time per iteration, because there may be O(V) verified paths and each verification may take O(V) time. On the other hand, we cannot guarantee that the matching will be augmented along more than one path. A scenario in which all the paths have the first edge in common (from the source to a vertex x ∈ L) and the minimum cost path has only this edge in common with the other paths is quite possible. Since O(V 2 ) may be a higher time complexity than that of computing the minimum cost paths, we may end up increasing the running time of the algorithm instead of decreasing it. Thus, we need to reduce the time complexity of the matching augmentation part. This can be achieved as follows. Let's remember that the minimum cost paths are paths in a shortest path tree (where the length of a path is its cost). We will consider the paths in the same order as before and we will consider all the edges to be initially unmarked. If the last edge of the current path is not marked then we will be able to augment the matching along the current path. After augmenting the matching along the current path, let x ∈ L be the first vertex on the path (after the source). We will traverse the whole subtree of the shortest path tree rooted at x and we will mark all the edges of this subtree. Augmenting the matching along all the possible paths takes at most O(V) time overall (because the paths are edge-disjoint). Marking the edges of the shortest path tree also takes at most O(V) time overall, because there are O(V) edges in the shortest path tree and each edge is marked at most once. Thus, the matching augmentation algorithm takes only O(V) time plus the time needed for sorting the paths in increasing order of their costs (e.g. O(V • log(V)) time).

The first version of the algorithm proposed in this section is described in pseudocode in Algorithm 3. The input to the algorithm consists of two maps: dist, containing the cost of the shortest path from S to every vertex x ∈ R (we consider dist(x) = +∞ if the vertex x is not reachable from S ), and parent, containing the parent in the shortest path tree for each vertex of the graph. Left-set vertices x are denoted as (x, L) in the algorithm and right-set vertices y are denoted as (y, R). In order to maintain the pseudocode simpler, we will mark the graph vertices instead of the edges (because, as mentioned earlier, in this case the edge-disjoint paths are also vertex-disjoint). The second version of the algorithm is described in pseudocode in Algorithm 4. The input to the Algorithm 4 also consists of the same two maps dist and parent, together with an extra map, children, which contains the children in the shortest path tree of each vertex of the graph.

This algorithm basically augments the matching along a maximal set of edge-disjoint paths (in fact, because they are paths of a shortest path tree, the paths are also vertex-disjoint except for the source vertex). It is important. however, to consider these paths in increasing order of their costs, in order to make sure that the residual graph at the next iteration does not contain negative cycles. The reason for which this optimization works is as follows. In a perfect mtaching every vertex has to be matched. When augmenting the matching along a shortest path to a node x, even if x is not the right-side node with the minimum cost path, the point is that any future minimum cost path to node x (in any future residual graph) will not have a lower cost than the current shortest path to x. So there is no reason for us not to augment the matching along that path, as long as this does not block paths with lower costs along which the matching could have been augmented.

For the problem presented in Section 6 we tested our optimization of not recomputing the perfect matching from scratch each time. The first minimum cost perfect matching algorithm that we used was the one which contained our matching augmentation optimization presented in Section 5 and tested earlier. We generated 10 bipartite graphs with n = 128 and edge costs randomly selected between 1 and 10000 (inclusive). We computed the total execution time and the total number of iterations of the successive shortest path algorithm. When the matching was computed from scratch each time (O(n 2 ) times) the total running time was 805 sec and the total number of iterations was 1305200. When we applied our optimization from Section 6, the total running time was 221 sec and the total number of iterations was 244886. When using the standard successive shortest path algorithm in order to compute a minimum cost perfect matching (and not using our optimization of not recomputing the matching from scratch each time), the total running time was 5798 sec and the total number of iterations was 20610157. When we applied our optimization from Section 6 and also used the standard minimum cost perfect matching algorithm the total running time was 232 sec and the total number of iterations was 251789. We can see that our optimization of not recomputing each minimum cost perfect matching from scratch is very effective. When combined with the optimization presented in Section 5, of augmenting the matching along multiple paths at each iteration, we obtained the best results. However, even when just the standard successive shortest path minimum cost perfect matching algorithm is used in conjunction with our optimization from Section 6 the improvements over the naive unoptimized version are significant. Nevertheless, more tests may need to be performed in the future in order to understand sufficiently well how good our proposed optimization really is.

Conclusions

In this paper we presented three practical algorithmic optimizations addressing problems like computing maximal matchings in induced subgraphs of grid graphs or computing minimum cost perfect matchings in bipartite graphs (under certain restrictions). The proposed optimizations were evaluated experimentally and compared against the unoptimized algorithms. The execution time was significantly reduced in each case, thus proving the validity and effectiveness of our optimizations.

Algorithm 2

 2 Greedy O(2 B • B+ M • N/B) Algorithm for Finding a Maximal Matching in an Induced Subgraph of a Grid Graph Compute the tables MRes, MCnt and BCnt.

Algorithm 3 Increasing the Matching Along a Maximal Set of Edge-Disjoint Paths -The O(V 2 ) Algorithm Input: dist, parent.

Set all the vertices of the graph as unmarked.

Sort the vertices x ∈ R in increasing order of dist(x). for x ∈ R in increasing order of dist(x) such that dist(x) < +∞ do y = (x, R), OK = true while y S and OK = true do if vertex y is marked then OK = f alse else y = parent(y) end if end while if OK = true then Increase the matching along the shortest path from S to (x, R) (the reverse of the path can be found by following the parent pointers starting from (x, R)). y = (x, R) while y S do Mark vertex y as marked. y = parent(y) end while end if end for Algorithm 4 Increasing the Matching Along a Maximal Set of Edge-Disjoint Paths -The O(V) Algorithm Input: dist, parent, children.

Set all the vertices of the graph as unmarked. Sort the vertices x ∈ R in increasing order of dist(x). for x ∈ R in increasing order of dist(x) such that dist(x) < +∞ do if (x, R) is not marked then Increase the matching along the shortest path from S to (x, R) (the reverse of the path can be found by following the parent pointers starting from (x, R)).

Let (y, L) be the first node on the path from S to (x, R) after S . Recursively mark all the vertices located in vertex (y, L)'s subtree of the shortest path tree. The children map will be used for retrieving for each vertex v the set children(v) that is the set of the shortest path tree children of the vertex v.

end if end for

Note that this optimization is not correct in a maximum matching algorithm. It is not correct to augment the matching along a path which does not have the globally minimum cost, because we are not sure if the right side vertex x needs to be in the optimal matching or not. And since vertices added to the matching are never removed by this algorithm, it is possible to make a mistake in this case.

This optimization does not change the theoretical time complexity of the minimum cost perfect matching algorithm, because we cannot provide any extra guarantees regarding the number of augmenting paths per iteration (and, thus, we cannot provide guarantees regarding the reduction in the number of iterations).

Minimum Cost Perfect Matching With the Minimum Cost Edge Ignored

In this section we consider the following problem. Given a complete bipartite graph with n nodes on the left side and n nodes on the right side and costs on its edges, we want to find a minimum cost perfect matching in which the cost is defined as the sum of the costs of all the edges in the matching except for the cost of the edge with the smallest cost.

A simple method for solving this problem is to iterate over all the edges (i, j) and fix them as the smallest edge in the matching. Then we would compute a (usual) minimum cost maximum matching in the bipartite graph from which left node i, right node j and all the edges (i , j ) with c(i , j ) < c(i, j) (or c(i , j ) = c(i, j) and the edge (i , j ) was considered before the edge (i, j)) are removed. If the maximum matching M has size n-1 then we found a potential solution, as follows. The potential solution consists of the n -1 edges of the found matching plus the edge (i, j). The cost of the matching (according to the definition used in this section) is equal to the sum of the costs of the n -1 edges of M. Note that the fixed edge (i, j) is the edge with the minimum cost in the perfect matching (the one whose cost is not considered towards the cost of the matching). Once the edge (i, j) was fixed we needed to minimize the total cost of the other n -1 edges of the perfect matching. Moreover, the other edges of the perfect matching needed to have costs which were larger than or equal to c(i, j). The minimum cost maximum matching M contains the n -1 edges we were looking for, in case its cardinality is n -1. If its cardinality is less than n -1 then we can conclude that there is no perfect matching containing the edge (i, j) as the minimum cost edge.

This solution requires the computation of O(n 2 ) independent minimum cost maximum matchings. The key to obtain a better solution is to notice that the O(n 2 ) matchings that we need to compute are not totally independent. We will sort all the n 2 edges first in ascending cost order and then we will consider them in this order. For the first edge we will compute the minimum cost maximum matching from scratch. Let's assume that we reached the edge (i, j). This time we will not compute the new matching from scratch. Instead, let's consider the matching M obtained for the previous edge in the sorted order. We will remove from M any edge with an endpoint at the left node i or the right node j (if any). Then we will remove from M all the edges with a cost smaller than c(i, j) (or equal to c(i, j) but for which the corresponding edge was considered before the current edge (i, j) in the sorted order). All the other edges of M will be maintained. We will start the (usual) minimum cost maximum matching for the edge (i, j) with all the remaining edges from M as part of the matching. Note that the algorithm may replace some of these edges by other edges. This can happen if the reverse of an edge (x, y) (x ∈ L and y ∈ R) from M is, at some point, part of the shortest path from S to T in the (new) residual graph. When considering the maximum matching problem as a maximum flow problem, the fact that an edge (x, y) is removed from the matching means that the flow is pushed back along that edge (in order to be redirected somewhere else).

By applying the optimizations from the previous paragraph we expect that the number of iterations required for computing each new minimum cost maximum matching will be significantly reduced.

This problem can also be viewed as a dynamic minimum cost perfect matching problem, in which the edge costs can be modified (for instance, instead of removing edges from the graph we can consider that their cost increased to +∞).

Experimental Results

We implemented the three optimizations presented in this paper and compared them against their unoptimized versions. All the tests were run on a machine running Windows 7 with an Intel Atom N450 1.66 GHz CPU and 1 GB RAM. All the algorithms were implemented in the C++ programming language and the code was compiled using the G++ compiler version 3.3.1.

First we tested our new algorithm for computing a maximal matching in an induced subgraph of a grid graph. We chose M = N = 2048 and we randomly generated the induced subgrapheach point (x, y) (0 ≤ x ≤ N -1, 0 ≤ y ≤ M -1) had an equal probability of being part of the subgraph or not. Thus, each of the tested subgraphs had approximately 50% of the nodes of the full grid graph. We generated 100 subgraphs and ran Algorithm 1 (the unoptimized version) and Algorithm 2 (the optimized version) on each of them. We computed the total running time for all the graphs. Algorithm 1 took 5.3 seconds. For Algorithm 2 we considered two values for B: B = 16 and B = 8. For B = 16 the running time was 2.74 sec and for B = 8 it was 1.13 sec. Note that in this case we computed the tables MCnt, Mres and BCnt each time (i.e. for each of the 100 tests). However, when running the algorithm on multiple tests with the same value of B, these tables only need to be computed once, in the beginning. Thus, we changed the algorithm in order to compute these tables only once in the very beginning and not for each of the 100 tests. The new running times were 0.92 sec for B = 16 and 1.13 sec for B = 8. Note that the running time is unchanged for B = 8 because the sizes of the tables are small and the time needed to compute them is negligible compared to the time needed to compute the matching. However, for B = 16, when the sizes of the tables increase significantly, it is much better to compute the tables in the beginning and reuse them for each test.

Then we repeated the tests for induced subgraphs of grid graphs containing 75% and 100% of the nodes of a full grid graph. For graphs with 75% of the nodes of a full grid graph the running times of our optimized algorithm were: 1.21 sec for B = 8 and 0.93 sec for B = 16 (note that we only considered the case when the tables are computed just once). The running time of the unoptimized version was 4.69 sec. When considering the full grid graph we obtained the following running times: 0.90 sec for B = 8 and 0.55 sec for B = 16 for the optimized version and 3.40 sec for the unoptimized version.

We did not test other values of B because the implementation would become less feasible. For B > 16 the sizes of the precomputed tables would become too large. For B = 8 and B = 16 we were able to make use of existing C/C++ data types (unsigned char and unsigned short int) in order to store a block. For B 8 and B 16 (and B ≤ 16) we cannot exactly fit a block into an existing C/C++ data type.

Second we tested the improvement brought by the use of multiple edge-disjoint paths for augmenting the matching in a minimum cost bipartite perfect matching algorithm. The expected improvement consisted in a reduction of the number of iterations. The standard algorithm would use n iterations where |L| = |R| = n. We chose n = 256 and we generated 100 complete bipartite graphs. The cost of each edge was chosen to be a random integer between 1 and 10000 (inclusive). The unoptimized algorithm we used was the standard successive shortest path algorithm with the Bellman-Ford-Moore algorithm for computing minimum cost paths at each iteration. The optimized algorithm simply included the O(n 2 ) matching augmentation along multiple edge-disjoint minimum paths described in Section 5. We measured both the total running time and the total number of iterations. The total running time (for all the 100 graphs) and the total number of iterations of the standard algorithm were: 23.17 sec and 25600 iterations. With our optimization the total running time was 2.3 sec and the total number of iterations was 1022. We notice that, even with the most basic implementation of our optimization, the running time was reduced 10 times and the number of iterations was reduced 25 times. Although the running time improvements may not translate directly when other minimum cost path computation algorithms are used (e.g. Dijkstra's algorithm) or when the O(V + V • log(V)) matching augmentation optimization is used (instead of the O(V 2 ) version), the improvement in the number of iterations does not depend on these algorithms and, thus, it is applicable to any implementation of the successive shortest path minimum cost bipartite perfect matching algorithm.

Then we considered the same testing scenario, except that the edge costs were chosen as random integers between 1 and 2 (inclusive). The total running time of our optimized algorithm was 2.98 sec and the total number of iterations was 5100. The number of iterations of the standard algorithm remained the same (as expected), but its running time dropped to 13.83 sec.

We also considered the complete bipartite graph with the following costs c(x, y) = min(x, y) (1 ≤ x, y ≤ n) and n = 256. In this case our optimization did not reduce the number of iterations at all (due to the special structure of the bipartite graph it was never able to augment the matching along more than one path per iteration). However, the running time with our optimization enabled was almost identical to the unoptimized version. We conclude that our optimization has a great potential for reducing the number of iterations of the successive shortest path minimum cost perfect matching algorithm and even in the pathological cases when it cannot reduce the number of iterations, it doesn't cause any significant overhead.

Although there is a large difference between the number of iterations obtained by our optimization in the cases of random complete bipartite graphs and in the case of the specific bipartite graph from the previous paragraph, we did not consider other types of bipartite graphs for testing. Understanding the correlation between the performance of our optimization and the specific properties of the costs of the edges of the bipartite graph is an interesting topic, but we defer its study to a later date, because we feel that this topic is more appropriate for a separate, more experimentally focused, paper.


Etienne Huot
email: etienne.huot@inria.fr

Hussein Yahia
email: hussein.yahia@inria.fr

Isaac Cohen
email: icohen@iris.usc.edu

Isabelle Herlin
email: isabelle.herlin@inria.fr

Matching Structures by Computing Minimal Paths on a Manifold

The general problem of matching structures is very pervasive in computer vision and image processing. The research presented here tackles the problem of object matching in a very general perspective. It is formulated for the matching of surfaces. It applies to objects having small or large deformation and arbitrary topological changes. The process described hinges on a geodesic distance equation for a family of curves or surfaces embedded in the graph of a cost function. This geometrical approach to object matching has the advantage that the similarity criterion can be used to define the shape of the cost function. Matching paths are computed on the cost manifolds using distance maps. These distance maps are generated by solving a general partial differential equation which is a generalization of the geodesic distance evolution scheme introduced by R. Kimmel, A. Amir, and A. F. Bruckstein (1995, IEEE Trans. Pattern Anal. Mach. Intell. 17, 635-640). An Eulerian level-set formulation is also introduced, leading to a numerical scheme used for solving partial differential equations originating from hyperbolic conservation laws, which has proven to be very robust and stable. C 2002 Elsevier Science (USA)

INTRODUCTION AND PREVIOUS WORKS

Matching structures is a major area for image processing and it involves a lot of classical applications such as stereo vision, temporal tracking, and matching of different acquisitions. The characterization of a matching function is an ill-posed problem in the sense that there is no unique solution. However, introducing structures properties such as their geometry or the underlying image representation allows the characterization of a unique matching function. Commonly used features are pixel gray level values for stereoscopic matching or optical flow [START_REF] Dhond | Structure from stereo -a review[END_REF], edges for token based approaches [START_REF] Van Doorn | Spatiotemporal integration in the detection of coherent motion[END_REF], and geometric properties of the structures [START_REF] Xu | Epipolar Geometry in Stereo, Motion, and Object Recognition: A Unified Approach[END_REF][START_REF] Zhang | Estimation of displacements from two 3D frames obtained from stereo[END_REF]. The latter properties are more robust since they can deal with situations where there is no consistency of image gray level. Relevant geometrical properties are selected on the basis of their ability to characterize a description of the structures which is invariant to the considered deformation. In the case of rigid or small elastic deformations high curvature points [START_REF] Cohen | Tracking points on deformable objects using curvature information[END_REF][START_REF] Serra | Optimal subpixel matching of contour chains and segments[END_REF] or semi-differential invariants [START_REF] Pajdla | Matching of 3D curves using semi differential invariants[END_REF] can be considered as an invariant description of the structure. Higher order geometric description such as crest and ridge lines can also be considered in order to characterize volume structures properties [START_REF] Thirion | Extremal points: definition and application to 3D image registration[END_REF][START_REF] Thompson | A surface-based technique for warping 3-dimensional images of the brain[END_REF]. Such methods perform well in the case of small deformation but cannot deal with large deformation and topological change. Such problems occur when we are interested in studying the evolution of a structure whose topology evolves in time [START_REF] Cohen | Curves matching using geodesic paths[END_REF].

In this paper we present a generalization to higher dimensions of a curve matching process described in [START_REF] Cohen | Curves matching using geodesic paths[END_REF], and we also generalize a curve evolution process introduced by Kimmel et al. [START_REF] Kimmel | Finding shortest paths on surfaces using levelset propagation[END_REF]. Curve evolution methods are becoming increasingly used in computer vision. Classical work of Kass et al. [START_REF] Kass | Snakes: Active contour models[END_REF] has been reformulated in the context of PDE-driven curves and surfaces (Malladi et al. [START_REF] Malladi | Shape modelling with front propagation: A levelset approach[END_REF], Caselles et al. [START_REF] Caselles | Geodesic avtive contours[END_REF], Gomes et al. [START_REF] Gomes | Level sets and distance functions[END_REF] etc.). In this work a new normal surface evolution is used to deduce a surface matching process.

The extension of the classical curve evolution process is presented in Section 2. The definition of the geodesic distance evolution equation is presented in Section 2.1. Sections 2.1.1 to 2.1.3 present the problem, introduce the Hodge operator, and introduce the evolution equation. Section 2.1.4 is devoted to deriving the equation of the projection, while Section 2.1.5 presents the level-set formulation. A surface matching is then deduced and described in Section 2.2. Section 2.3 presents a result related to meteorological data.

SURFACE MATCHING

This section addresses the problem of surface matching. The method proposed is based on a generalization of the algorithm introduced by Cohen et al. in [START_REF] Cohen | Curves matching using geodesic paths[END_REF] for curves. This scheme has also motivated a generalization of the geodesic curve propagation method [START_REF] Kimmel | Finding shortest paths on surfaces using levelset propagation[END_REF] into a surface evolution framework that is used for deriving a cost hypersurface W ⊂ IR 4 measuring the similarity of the surfaces to be matched. We begin by setting up a geodesic distance evolution scheme for evolving surfaces.

Geodesic Distance

2.1.1.

A geodesic distance evolution rule for propagating surfaces on a 3-manifold. Let X be a 3-manifold (or hypersurface) in IR 4 . We suppose X compact and pathwise connected. 1From these assumptions, one can derive, using an easy application of the Hopf-Rinow-De Rham Theorem (see [START_REF] Spivak | A Comprehensive Introduction to Differential Geometry[END_REF]) that given any two points M 0 and M 1 ∈ X , there exists a path γ :[0, 1] → X connecting M 0 and M 1 (γ (0) = M 0 ,γ(1) = M 1 ) whose length minimizes the lengths of all paths between the two points. The length of γ is called the geodesic distance between M 0 and M 1 and will be denoted d X (M 0 , M 1 ) in the remainder of the paper. Moreover, the path γ is necessarily a geodesic curve on X , i.e., a curve such that the second derivative d 2 γ/du 2 is always perpendicular to X (Where u is the arclength of γ ). Let Y ⊂ X be a surface (2-manifold) "drawn" on X . We consider the surfaces t ⊂ X whose points are located at geodesic distance t from Y:

t ={M ∈ X | d X (M, Y) = t}.
(

) 1
We are interested in defining a partial differential equation governing the evolution of the surface t as the parameter t evolves. For this purpose, we will need a notion of cross-product in 4-space and a method of deriving simple formulae about such a cross-product. These formulae are needed in the demonstration of intermediate propositions. One mathematical tool that can achieve these requirements is given by the Hodge * operator, a notion recalled in the next section.

Exterior algebras and the hodge * operator.

The theory is only briefly reviewed here. The reader is referred to [START_REF] Abraham | Manifolds, Tensor Analysis and Applications[END_REF] for a complete presentation of the subject. Let p (IR n )be the p-exterior power of space IR n and * the Hodge operator defined, for every λ ∈ p (IR n ) by the equality:

λ ∧ µ = * λ, µ n-p e 1 ∧ e 2 ∧•••∧e n . (2)
In Eq. ( 2) (e 1 ,...,e n ) is the standard basis of IR n , , p is the usual dot product on p (IR n ) defined on generators by

u 1 ∧•••∧u p ,w 1 ∧•••∧w p p = det( u i ,w j ), (3)
and µ ∈ n-p (IR n ). From the property * ( * λ) = (-1) p(n-p) λ, one can deduce that, if u,v, and w are three linearly independent vectors in IR 4 , the associated image * (u ∧ v ∧ w) ∈ 1 (IR 4 ) = IR 4 satisfies the following properties:

• it is a vector in IR 4 perpendicular to u,v, and w.

• The basis (u,v,w, * (u ∧ v ∧ w)) is positively oriented.

• Its squared norm:

* (u ∧ v ∧ w) 2 is equal to u, u u,v u,w v, u v, v v, w w, u w, v w, w .
Using the Hodge * operator, we will derive the geodesic distance evolution scheme for the family of surfaces α(u,v,t). This is presented in the following section.

The geodesic distance evolution equation.

We will make use of the following fact, proved in [7, p. 183]: LEMMA 2.1. Around any point on a surface one can always find a local orthogonal parameterization.

So let α(u,v,t) be a local orthogonal parameterization of a family of surfaces in X, and let τ u = α u / α u , τ v = α v / α v denote the two unitary tangent vectors determined by α and N the normal vector to X . We suppose that the family α(u,v,t) satisfies the following partial differential equation:

∂α ∂t = * (N ∧ τ u ∧ τ v ). (4)
We want to prove that, for each t,α(u,v,t) is a local parameterization of t . Let β(t)be the curve in X defined by β(t) = α(u,v,t)| u=u 0 ,v=v 0 ,fixed . Then: LEMMA 2.2. For any u 0 ,v 0 , the curve β(t) is a geodesic in X .

Proof. We prove this lemma by showing that β tt is perpendicular to τ u , τ v , and * (N ∧

τ u ∧ τ v ).
The only possibility that will remain will be that β tt is collinear to N which means that β is a geodesic (by definition). First note that

β t = dβ dt = ∂α ∂t = * (N ∧ τ u ∧ τ v ), so β t 2 = * (N ∧ τ u ∧ τ v ) 2 = 1 -N, τ u 2 -N, τ v 2 = 1
since the parameterization is orthogonal. This proves:

β tt ,β t = d 2 β dt 2 , dβ dt = d dt ( * (N ∧ τ u ∧ τ v )), * (N ∧ τ u ∧ τ v ) = 0. (5)
Hence, β tt is perpendicular to * (N ∧ τ u ∧ τ v ). Now, using methods similar to the one used in the 2D case (see proof of Lemma 1 in [START_REF] Kimmel | Finding shortest paths on surfaces using levelset propagation[END_REF]), one can show that β tt , τ u =0, and similarly β tt , τ v =0.

The curve β(t) is a geodesic in X .

LEMMA 2.3. The evolution of the family of surfaces t is given by the equation

∂α ∂t = * (N ∧ τ u ∧ τ v ). ( 7
)
The proof of this lemma is not given here, as it is a simple adaptation, using the the general Gauss lemma [START_REF] Spivak | A Comprehensive Introduction to Differential Geometry[END_REF] of the proof given in [12, Lemma 2].

The two preceding lemmas demonstrate that, using the Hodge * operator, it is possible to derive a geodesic distance evolution scheme for a family of surfaces described by local orthogonal parameterizations. In the next section, we compute the normal speed of the projection of t onto the (x, y, z) hyperplain in IR 4 . 7) cannot be solved directly in a general context. We follow the method described in [START_REF] Kimmel | Finding shortest paths on surfaces using levelset propagation[END_REF] by considering the evolution of the projection. We now make the assumption that X is a graph hypersurface W , that is to say X = W ={(x, y, z,w(x, y, z))} for a function w : IR 3 → IR . Let π : IR 4 → IR 3 be the canonical projection onto the (x, y, z) hyperplan in IR 4 and let S(t) be the projection of the image of α(u,v,t) (that is to say, t ) onto that hyperplan:

Geodesic distance computation. Equation (

C(t) = π • α.
We denote by p, q, r the following quantities: p = ∂w ∂ x , q = ∂w ∂ y , and r = ∂w ∂z . Starting from a result mentioned in [START_REF] Epstein | The curve shortening flow[END_REF], we admit that the trace of a propagating surface may be determined only by its normal velocity, as the other components of the velocity influence only the local parameterization. Our goal is then to compute the projected velocity of the evolving surface F = π • α t , n , n = (n 1 , n 2 , n 3 ) being the normal to the projected surface π • α(u,v,t). By using the same assumptions as in the bidimensional case, one can write the propagation of the projected surface C(t) along its normal by

∂C ∂t = Fn, ( 8
)
where the function F represents the propagation velocity of the surface along its normal. This velocity can be obtained by

F = n,π • α t ,
where τ u and τ v are the tangential vectors of α defined by

τ u = α u α u = (x u , y u , z u ,w u ) x 2 u + y 2 u + z 2 u + w 2 u and τ v = α v α v = (x v , y v , z v ,w v ) x 2 v + y 2 v + z 2 v + w 2 v ,
and N the normal vector to W :

N = * (W x ∧ W y ∧ W z ) * (W x ∧ W y ∧ W z ) = (-p, -q, -r, 1) 1 + p 2 + q 2 + r 2 .
The propagation equation becomes

α t = 1 K × qz u w v -qw u z v -ry u w v + ry v w u + y v z u -qz u w v + qw u z v + ry u w v + y u z v -ry v w u -y v z u py u z v -py v z u -qx u z v + rx u y v + qx v z u -rx v y u additional term with K = 1 + p 2 + q 2 + r 2 x 2 u + y 2 u + z 2 u x 2 v + y 2 v + z 2 v .
Using the fact that the parameterization is orthogonal, one finds:

F = π • α t , n = an 2 1 + bn 2 2 + cn 2 3 -dn 1 n 2 -en 1 n 3 -fn 2 n 3 (9)
with

a = 1 + q 2 + r 2 1 + p 2 + q 2 + r 2 , b = 1 + p 2 + r 2 1 + p 2 + q 2 + r 2 , c = 1 + p 2 + q 2 1 + p 2 + q 2 + r 2 , d = 2 pq 1 + p 2 + q 2 + r 2 , e = 2 pr 1 + p 2 + q 2 + r 2 , f = 2qr 1 + p 2 + q 2 + r 2 .
Having found the normal equation evolution of the projected surface, we now proceed to set up an Eulerian formulation for S(t), by writing that projected surface as a level-set ϕ -1 (0). We derive such a formulation in the next section.

Level-set formulation.

Given a function ϕ : IR3 → IR such that its zero level-set tracks the projected surface S(t) = ϕ -1 (0), one can determine that the function ϕ follows2 the propagation equation

∂ϕ ∂t = a ∂ϕ 2 ∂ x + b ∂ϕ 2 ∂ y + c ∂ϕ 2 ∂z -d ∂ϕ ∂ x ∂ϕ ∂ y -e ∂ϕ ∂ x ∂ϕ ∂z -f ∂ϕ ∂ y ∂ϕ ∂z , (10)
the coefficients a, b, c, d, e, and f having the same values as in Eq. ( 9). As mentioned in the curve evolution process described in [START_REF] Kimmel | Finding shortest paths on surfaces using levelset propagation[END_REF] such as Eulerian formulation leads to numerical resolution schemes able to handle problems caused by a time varying coordinate system (u,v,t): curvature singularities and topological changes [START_REF] Sethian | Level Set Methods: Evolving Interfaces in Geometry[END_REF]. We describe, in the following section, the numerical resolution method used to solve Eq. ( 10). Then we use this numerical algorithm to build distance maps.

The numerical implementation is a generalization of the finite difference approximation described in [START_REF] Kimmel | Finding shortest paths on surfaces using levelset propagation[END_REF][START_REF] Osher | Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF] for Hamilton-Jacobi type equations. It consists in an explicit temporal scheme where spatial derivatives are approximated by finite differences using the minmod function [START_REF] Osher | Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Sethian | Level Set Methods: Evolving Interfaces in Geometry[END_REF][START_REF] Sethian | Crystal growth and dendritic solidification[END_REF]. The derivative estimates can be bounded and the variations of the solution can therefore be controlled.

Using this numerical scheme, we can now use the computed function ϕ to generate distance maps on a manifold. In order to use Eq. [START_REF] Huot | Surface matching with large deformations and arbitrary topology: a geodesic distance evolution scheme on a 3-manifold[END_REF] for computing the geodesic distance map of the surface discribed by a parameterization α(u,v,0) on the 3-manifold W ,wehave to define an initial estimate ϕ 0 such that the initial surface is represented through a level-set of ϕ 0 . This initial estimate can be obtained in several ways according to the data. We use a Euclidean distance map [START_REF] Borgefors | Distance transformations in arbitrary dimensions[END_REF] in such a way that:

ϕ 0 (x, y, z) =      -d(x, y, z)i f ( x, y, z) is inside ϕ -1 0 (0) 0i f ( x, y, z) ∈ ϕ -1 0 (0) +d(x, y, z)i f ( x, y, z) is outside ϕ -1 0 (0). (11)
Given a graph hypersurface X = W and the initial estimate ϕ -1 0 (0) on this hypersurface, Eq. ( 10) characterizes the distance map of the area whose boundary is defined by ϕ -1 (0).

Surface Matching

We can now use the previous theory to generalize the curve matching method presented in [START_REF] Cohen | Curves matching using geodesic paths[END_REF] to a matching process between two arbitrary surfaces S and D in IR 3 . The two surfaces S and D are represented as 0-level-sets of two functions ϕ 0 and ψ 0 . The two functions ϕ 0 and ψ 0 are computed from the initial data using the rule presented in [START_REF] Kass | Snakes: Active contour models[END_REF]. Then given these two initial functions ϕ 0 and ψ 0 , the numerical process presented in Eq. ( 10) is used to generate two distance maps on the graph surface W : D S and D D .

In this section we introduce the definition of the manifold W and how we use the distance maps to solve the surface matching. It is important to understand that the similarity criterion used in the matching process relies completely on the choice of W : each choice of W leads to a different matching result. In the following, we give an example for the definition of W .

Definition of a cost hypersurface modeling geometric properties.

To take into account the geometrical properties of the surfaces to be matched, we use a graph hypersurface which takes into account the mean curvature of surfaces. This definition uses the curvature information in a small neighborhood of the source and destination areas. Within this neighborhood curvature information is relevant since only small deformation occurs and its curvature measure may be used to strengthen the matching of similar points. We characterize the local properties of the structures through the curvatures variation of the curves: κ = κ Sκ D , where κ S and κ D are respectively the mean curvatures of S and D.Inthe case of small deformation we want the graph W to reflect the similarity between points having similar curvatures. This property will be down-weight as the distance increases. For that purpose, one can define ρ as

ρ( κ, d) = 1 - x 2 1 + d 2 x 2 /σ , ( 12
)
where σ is a scale parameter and d the Euclidean distance. The use of such a function to measure the similarity of the geometric properties allows reject or down-weight large dissimilarities. A set of robust error functions having the same behavior can be found in [START_REF] Robert | Dense depth map reconstruction: A minimization and regularization approach which preserves discontinuities[END_REF]. As the dissimilarities grow, the function ρ decreases and the function ρ approaches zero.

The graph hypersurface is then given by the following equation:

W = (x, y, z, min(|ϕ 0 |ρ( κ, |ϕ 0 |), |ψ 0 |ρ( κ , |ψ 0 |))). ( 13
)
|ϕ 0 | and |ψ 0 | are the Euclidean distance from respectively S and D. The mean curvature is easily computed from the level-set representation of the surfaces:

κ = (ϕ yy + ϕ zz )ϕ 2 x + (ϕ xx + ϕ zz )ϕ 2 y + (ϕ xx + ϕ yy )ϕ 2 z -2ϕ x ϕ y ϕ xy -2ϕ x ϕ z ϕ xz -2ϕ y ϕ z ϕ yz ϕ 2 x + ϕ 2 y + ϕ 2 z 3/2 .
To smooth irregularities coming from real data or singularities in the graph W , we perform a gaussian convolution on the graphs W . This guarantees that the partial derivatives ∂w ∂ x , ∂w ∂ y , and ∂w ∂z are always defined.

Similarity measure and matching paths.

Let S and D be two surfaces drawn on W and identified with "source" and "destination" surface. We want to compute minimal paths, on the manifold W , connecting S and D. Wherever W is considered as a cost hypersurface such a path corresponds to paths of minimal cost between the two structures. Computing such a path amounts to searching for an optimal path among all the paths p X S starting at X S ∈ S and ending at a point X D ∈ D and minimizing a given cost function f . A cost function f can be easily derived using the fact, proven in [START_REF] Kimmel | Finding shortest paths on surfaces using levelset propagation[END_REF], that the minimal path between S and D minimizes the sum of the distance maps. So we consider the cost function given by f (x, y, z) = ϕ(x, y, z) + ψ(x, y, z), [START_REF] Osher | Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF] where ϕ and ψ are the bivariate functions defining the distance maps D S and D D . Note that the function f may take negative values, as we take into account the interior and exterior of the surfaces. To define the matching paths we use a property of the graph surfaces D S and D D which relates the equal distance contour to minimal paths: the minimal paths are orthogonal to equal distance contours. Since the equal distance contours are level sets of the map D S + D D , a vector tangent to the matching path is defined by the gradient of the cost function:

∇ f =∇ϕ +∇ψ.
This approach gives a reliable path construction scheme: given a point X S on the source area, move this point in the opposite direction of the gradient of f until reaching a point on the destination area. The result is the minimal cost path connecting the two areas and starting at X S . This path p(s)isdefined by the parameterized curve p(s) such that p(0) = X S , p(1) = X D and

∂ p ∂s =-∇(ϕ + ψ),
where X S ∈ ϕ -1 0 (0) is given and X D ∈ ψ -1 0 (0) is unknown. These orbits are obtained by applying a classical Runge-Kutta scheme.

Meteorological Result

Results displayed in the figures come from an experiment led on a Meteosat temporal sequence of an atmospheric depression in the tropics. The Meteosat satellite has a 30 minutes acquisition's frequency. In some situations, this temporal sampling is not sufficient to characterize the structures evolution and a model has to be considered. For example, in the tropics clouds structures aggregate and disaggregate rapidly. observed on these data, it has been proved that cloud temperature is directly linked to the elevation hence leading to a measure of cloud elevation. This example is very interesting since the structure may deform rapidly and it displays large and small deformations at specific locations and changes of topology simultaneously. The matching is performed between the elevation maps of the two cloud structures extracted from Fig. 1. Figure 2 shows the result of the matching method. The leftmost cloud structure of Fig. 1 is clearly depicted with its two components. Above is the cloud structure corresponding to the elevation map of the rightmost image in Fig. 1. The matching paths are depicted in the two views of Fig. 2.

CONCLUSION

The research presented in this work introduces a new method for answering the problem of matching two objects. It relies on a geodesic distance evolution scheme written in the form of a partial differential equation for computing distance maps. The surface evolution scheme is a generalization of the geodesic curve evolution process presented in [START_REF] Kimmel | Finding shortest paths on surfaces using levelset propagation[END_REF], leading to an algorithm for computing distance maps in 3-manifolds. The matching paths minimize a cost criterion which can incorporate various properties. The resulting matching algorithm allows us to take into consideration large deformations and topological changes between the source and destination surfaces. Use of the Hodge * operator for deriving the surface propagation rule unveils a generalization of the matching process to any number of dimensions [START_REF] Huot | Surface matching with large deformations and arbitrary topology: a geodesic distance evolution scheme on a 3-manifold[END_REF]. The level of generality of the matching process presented in this paper should give new insights to many problems in computer vision where the matching of complex structures is relevant.

FIG. 1 .

 1 FIG. 1. Two consecutive images extracted from the same METEOSAT sequence.

FIG. 2 .

 2 FIG. 2.Matching paths between the two meteorological structures.

These assumptions do not restrict the validity of the theory presented in this work, as the manifold X will appear as the graph of a cost function, which automatically satisfies these requiremenrs in practice.

It is important to note that the function ϕ depends not only on (x, y, z) ∈ IR

, but also of the parameter t.T o simplify the notations, we do not write explicity that dependence on t, but it is important to keep it in mind.