Wind Corrections in Flight Path Planning

Abstract When operating autonomous unmanned aerial vehicles (UAVs) in real environments it is necessary to deal with the effects of wind that causes the aircraft to drift in a certain direction. In such conditions it is hard or even impossible for UAVs with a bounded turning rate to follow certain trajectories. We designed a method based on an Accelerated A* algorithm that allows the trajectory planner to take the wind effects into account and to generate states that are reachable by UAV. This method was tested on hardware UAV and the reachability of its generated trajectories was compared to the trajectories computed by the original Accelerated A*.

Introduction
This paper addresses the problem of three-dimensional trajectory path planning for UAVs with a bounded turning rate in the presence of constant wind, obstacles and other aircraft. To represent the motion dynamics of airplane-like UAVs, we apply constraints on flight manoeuvres and restrictions on the smoothness of the trajectory (trajectory and its first derivation are continuous).
In real environments UAVs are significantly susceptible to being drifted by wind. Particularly, aircraft with a bounded turning rate require different turning radii depending on the flight direction relative to the wind vector. For a successful and precise execution of a planned trajectory, the trajectory has to reflect the physical capabilities of the aircraft and the assumed impact of wind.
In this paper we present an approach based on an Accelerated A* (AA*) algorithm [11] for collision-free path planning in an environment with no-flight zones and other aircraft, which is capable of generating physically satisfiable plans in the presence of constant wind. The AA* algorithm was chosen because it is both complete and optimal and thanks to adaptive sampling it is fast even for large environments with a high number of obstacles.
A*-like algorithms are primarily intended to be used in fully observable deterministic environments. Real environments are non-deterministic, however we assume that the UAV has a GPS receiver to correct its position estimation and minimize the effects of non-deterministic actions. We also assume a prior knowledge of a map of the environment with obstacles defined as no-flight zones, thus considering the environment fully observable.
Our method has been successfully deployed and tested on a real fixed-wing UAV using Kestrel autopilot with significantly better results than those achieved by regular AA* algorithms in terms of plan execution precision.
The paper is organized as follows. In the following section we summarize related work dealing with path planning and the effects of wind. In Section 3 we describe the AA* algorithm our work is based on. Our modified version of the AA* algorithm is introduced in Section 4. In Section 5 we present the results we achieved and Section 6 concludes the paper.

Related work
The area of UAVs' trajectory planning has been the subject of many research works over the past several years. Numerous methods for flight trajectory planning in no-wind cases have been proposed. There are algorithms based on randomness, e.g., random-walk [3], rapidly exploring random trees (RRT) [8] or randomized potential fields [2]. These methods are fast but they are still based on random sampling and generally cannot find optimal solutions with respect to given criteria.
There are also algorithms computing optimal solutions that exploit pre-computed structures for a given environment definition -potential field path planning [7], vector field methods [9] or the 3D Field D* algorithm [4]. However, these approaches need to recalculate their inner data structures every time the environment changes.
The impact of wind effects on trajectory planning was dealt with in [10] by iteratively solving a no wind case problem with a moving virtual target, or in [13] by generating candidates of smooth (continuous curvature) time-optimal paths between initial and final points in the plane. These approaches search only for a horizontal plan between the start and final configuration and do not take possible obstacles into account. Cobano et al. [5] used a genetic algorithm and the Monte-Carlo method to deal with wind effects and to solve collisions noncooperatively. However this is also a randomized approach that does not have to find optimal solution and is not suitable for cooperative collision avoidance methods.
Our approach is based on an AA* algorithm, which can, thanks to adaptive sampling, efficiently (in a short time) find optimal paths in three dimensional environments with obstacles. The algorithm was modified to take the effects of wind and the limitations of the aircraft into account and is suitable for cooperative collision avoidance.
3. Accelerated A* Accelerated A* (AA*) is a modified version of the A* algorithm presented in [11] that overcomes the trade-off between search speed (efficiency) and search precision by using an adaptive sampling mechanism. The algorithm generates samples of state space using elementary flight manoeuvres where the length of the manoeuvres is given by the aircraft's distance to the nearest obstacle -the further the obstacle the sparser the sampling and vice versa (see Figure 1). Search precision is defined by a minimal sampling step. The AA* algorithm finds optimal trajectories only with respect to the sampling used. However, thanks to the use of an adaptive sampling mechanism (dependent on the distance to the closest obstacle), AA* is in this domain typically faster than the original A* algorithm with uniform sampling, which allows the use of a much smaller minimal sampling step while preserving the low computational requirements.
Manoeuvres that define the sampling and correspond to the aircraft's elementary actions are straight elements that are given by the start configuration and length, horizontal turn and vertical turn that are defined by the turning radius, arc angle and turn direction -left or right, respectively, and up or down.
The algorithm works in the same way as regular A*. It removes the best state from the OPEN list and if it can be directly connected with the goal state, it smoothes the found path by application of Dubins manoeuvres [6] when possible and returns the plan. Otherwise, it expands the current state by applying actions represented by the flight manoeuvres, estimates the new position and new sampling, checks the presence of this new state in the CLOSED list and checks the validity of the element, counts heuristics using Dubins curves and inserts the new state into the OPEN list. The sampling of continuous space can generate an infinite number of states where none are exactly the same. To prevent this, two states are considered equal if their Euclidean distance and their direction vectors differ less than a specified value given by the sampling parameterization used.

Modified AA* for planning in wind
Wind influences the minimal turn radius and ground speed of UAVs depending on the flight direction relative to the wind. The proposed AA* planning algorithm does not take wind into consideration. It does not allow dynamic changes of turn radius and does not take the ground-and air speeds into account during the planning process. Thus, in windy conditions the planned trajectories are hard or impossible to follow.
The impact of wind on UAV flight is most significant in the horizontal plane because of (i) better manoeuvrability of a fixed wing UAV in the vertical plane and (ii) slower vertical winds in the lower layers of the atmosphere. That is why we modelled wind as a two-dimensional vector and modified only the horizontal properties of the AA* algorithm.

Geometry of flight in wind
Aeroplane-like UAVs can turn to both sides with a limited turn radius and in the presence of wind they drift in its direction. That is why the geometrical model of turn manoeuvres was changed from a turn along circle with a constant radius to a turn along trochoid and straight elements that changed according to wind speed and relative wind direction.
A trochoid is a locus of a point (x, y) orbiting at a constant rate � with radius r around an axis located at �� � � � � � with parameterization: with the axis being translated in the xy-plane at a constant rate in a straight line.
When applying turn manoeuvres the aircraft still moves in circles in the air mass, however the presence of wind makes its movement trochoidal according to the ground. As shown, for example, by McGee et al. in [10], the shape of the flight curve depends only on the aircraft's turn radius and on the wind speed to aircraft's air speed ratio. Figure 2 demonstrates the effect of wind on a trajectory of an aircraft flying in circles that drifts in the direction of the wind. � � is the airspeed vector and � � is the wind vector. In [1] it was shown that if the wind speed to airspeed ratio (labelled as β) is less than one, then any two points with arbitrary orientations are reachable for a vehicle with Dubins type kinematics operating in the presence of wind. Thus, our approach should be valid for wind speeds up to the maximal airspeeds of the UAV.

Changes to AA* algorithm
To implement the geometry of the flight manoeuvres into the AA* algorithm we needed to change the way a new manoeuvre is applied to the current state and the way two state configurations are connected by turn and straight elements. Also, we needed to modify the heuristic function to reflect the effects of wind on the optimality of Dubins manoeuvres. The rest of the algorithm remained the same.
The startAngle parameter determines the initial position in the trochoid (this corresponds to an angle in a circle, see Figure 5) where the curve is tangent to the state's airspeed direction vector. The angle parameter corresponds to the manoeuvre's position shift along the trochoid measured from a given state's position (see Figure 5). Parameter r is an aircraft's turn radius and � ��� and � ��� are corresponding components of the wind vector. � � and � � specify the position of the trochoid in space. The angleXY is the aircraft's heading used to estimate the startAngle according to the following equation: for left and right turns respectively.
We know the wind vector that we consider constant during the whole planning process. From the known track's heading (ground speed direction vector), which needs to be continuous and the known aircraft's airspeed (direction and aircraft's heading is unknown) we can estimate angleXY according to Figure 3 and the following equations. We substitute angleXY to Equation (5) and then to Equations (3) and (4) to correctly place the trochoids (see Figure 4). Obviously: where Vg is ground speed. Then we can obtain the ground speed magnitude as: By projecting the wind vector to the ground speed vector we get: and from the right-angled triangles we get: By substituting (8) and (9) to (7) we estimate the ground speed magnitude and calculate the required airspeed direction (angleXY) as:

New straight flight manoeuvre
The effect of applying the straight flight manoeuvre in wind is straightforward and can be expressed by Equation (6). New straight flight and turn manoeuvres are shown in Figure 4 for different values of an aircraft's heading.

Connecting two states
When planning a path between two state configurations without obstacles, we use a modification of Dubins curves where the turn elements are trochoidal curves connected with a straight element (see Figure 5(b)). This connection of three manoeuvres is called a Smooth manoeuvre and is used for the estimation of heuristics and for smoothing the planned trajectory.
The initial trochoid is placed so that it passes through the initial state and is tangent to its direction vector. The startAngle is estimated accordingly. The second trochoid is placed similarly with respect to the second state. Next, we need to find a tangent to both trochoids representing the straight manoeuvre between the turns. If there are more solutions to these manoeuvres we select the shortest one. The mechanism of finding the tangents differs according to the turn sides in the three consecutive manoeuvres. We distinguish two cases: i) turn -straight segment -same side turn and ii) turn -straight segmentopposite side turn (see Figure 6).

Same side turn
As can be seen in Figure 6(a) in this case the trochoids are only translated one to the other where the direction of the straight segments connecting these two curves is parallel to a translation vector � � � � � � ��. Because a trochoid is a periodic function, there are an infinite number of translations vectors. The directions on the ends of straight segments are the same (� � � � � ), so we can determine the translation vectors by the substitution � � � � � � � � ��.
Then we need to find the value of angle �, which specifies the first turn length, so that one of the translation vectors is tangential to a point on the first curve given by �. From that point we can apply a particular translation and find the straight manoeuvre that tangentially connects the two trochoids. This method is described in detail in [16].

Opposite side turn
In this case the trochoids are centrally symmetric according to a point � � �� � � � � and again there are an infinite number of these points (see Figure 6(b)). The directions on the ends of the straight segment are the same, but either the first or the second turn radius is negative (� � � �� � ). We can determine the centres of symmetry by the substitution � � � �� � � � ��� � ���.
Straight manoeuvres need to pass through one of these points. Then we need to find the value of parameter � so that the following straight element passes one of the centres of symmetry.

Modifications of heuristic function
As described in Section 3, the AA* algorithm uses Dubins curves to estimate the heuristics from a given state to the goal. However, McGee et al. [10] proved that for certain configurations in the presence of wind, non-Dubins manoeuvres need to be applied for time-optimal planning.
To be able to use the manoeuvres presented above and not to reduce the optimality of the algorithm, we defined a so-called irregular zone around a current state. If the target state lies outside the zone, we can be sure that our manoeuvres are time-optimal. For states inside this zone we are not sure whether they can be reached by Dubins manoeuvres and thus we use other less informative heuristics, such as the Cartesian distance. Another effect of this approach is that outside the irregular zone, only the curve-straight-curve elements can be time-optimal and thus we do not need to consider the curve-curvecurve Dubins elements. The size of the irregular zone and its description requires complex derivation, which is given in Appendix A.1.

Experiments
To verify the proposed methods, we deployed the modified AA* algorithm on a real fixed wing UAV (see Figure 7). We compared plan execution precision in cases with plans found by a regular AA* algorithm and with plans found by the AA* modified for planning in wind. Tests were conducted in a the presence of approximately 5m/s wind in different directions (measured by Kestrel autopilot on the basis of airspeed and GPS comparison) and the UAV's airspeed was 15m/s, giving � � ����. The scenario consisted of three waypoints placed in a triangle (see Figures 9 to 12). The blue line in Figures 9-12 represents the planned trajectory that starts at the first waypoint (WP1) then proceeds via WP2 and WP3 and ends again in WP1. The green line shows the real recorded trajectory of the UAV and black arrows emphasize different turn radii in individual cases. Red crosses are plan samples uploaded to the UAV's autopilot. The sampling period is proportional to the planned flight time in straight segments and the turning rate in turn manoeuvres. The largest deviations from the plan occurred near the third waypoint where the UAV turns to the direction of the wind. The speed of the wind and the aircraft add together and thus, a much larger turn radius is needed. This makes the trajectories created with Dubins curves with a constant turn radius hard to follow. Trajectories created with trochoidal curves, adapted to the wind speed and direction, are followed much more precisely. On the other hand, when the radius of the scheduled turns is larger than necessary, the trajectory gets shifted inwards of the turns (noticeable in Figure 11). This behaviour is caused by a navigational feature of the autopilot that switches the current waypoint target to the next one whenever the UAV reaches a predefined boundary around the waypoint (in our case 50m). This feature makes waypoint passing easier but causes the UAV's track to be shifted inwards of the turns. A modified planner can partially deal with this problem as well, because it schedules sharper curves whenever possible. Other errors are caused mainly by dynamic changes of wind during the flight. The reason why we show only one flight over the same area in Figures 9 and  10 is that wind conditions during the experiment were not stable for a longer period of time. Figure 8 shows the probability that the track following error is greater than the given distance in both East wind situations.

Conclusion and Future Work
In this paper we presented a novel approach for trajectory planning for Unmanned Aerial Vehicles (UAVs) in the presence of wind, based on a modified Accelerated A* algorithm. In contrast to any-time methods, the AA* approach optimally plans the whole UAV's trajectory. This modified algorithm was deployed and tested on a real UAV at a wind speed of 5m/s and achieved significantly better results than the original AA* in terms of the precision of plan execution.
In future research we would like to give evidence that in contrast to any-time methods, the AA* approach is more suitable for purposes of distributed cooperative trajectory conflict resolution, as stated in [12]. Furthermore, we would like to address the following topics: (i) modelling the uncertainty of plan execution for safer collision avoidance and (ii) modelling wind changes in different locations and altitudes. Currently, we are working on temporal planning in wind on more sophisticated missions.

Acknowledgements
The research described in this paper was supported by the Czech Ministry of Defence grant OVCVUT2010001, the U.S. Army grant W911NF-11-1-0252 and the Grant Agency of the Czech Technical University in Prague grant SGS13/078/OHK3/1T/13.

A.1 Optimality of modified manoeuvres for heuristics
In [10] the authors transform the problem of finding optimal paths in wind by iteratively solving a no wind case Dubins car problem (with all six manoeuvres -RLR, LRL, RSR, LSL, RSL, LSR) to intercept a moving virtual target, which moves in a straight line by the speed of the wind in the opposite direction. They describe points on this line by variable d, which is the distance to the initial target position. They define two functions � � ��� and � �� ��� where � � is the time it takes the aircraft to fly the optimal air path from its initial position and orientation to a final orientation at point d and � �� is the time it takes the virtual target to move to point d. The authors further prove that in some cases (when there is a discontinuity in G(d)) it is necessary to use longer than optimal manoeuvres (in a no-wind case) or special non-Dubins manoeuvres to reach the optimal interception point on the path of the virtual target.
In our approach we do not use three turn elements (RLR, LRL -in the next text labelled as CCC, meaning curvecurve-curve) or non-Dubins manoeuvres but instead we define a zone (irregular zone) around the UAV, outside of which, we are sure that the optimal path will consist only of Dubins curves with two turns (RSR, LSL, RSL, LSR LRL -in the next text labelled as CSC, meaning curvestraight-curve). Inside that zone we use other heuristic functions than the length of the optimal Dubins curve, e.g., the Cartesian distance.

A.1.1 Construction of the irregular zone
To prevent cases where non-Dubins manoeuvres are necessary (in case of discontinuities in G(d)) the irregular zone should be large enough to cover cases where the RSL and LSR manoeuvres do not exist. These discontinuities occur when the outer circles have centres more than �� ��� apart, which in the worst case is 4� ��� distance from the aircraft. Furthermore, according to [15] CCC manoeuvres can be optimal ones only for points closer than 4� ��� .
In the presence of wind, the target moves and thus a specific part of its trajectory needs to be checked for its presence in the irregular zone. Movement of the virtual target is depicted in Figure 13, where S is the initial state of the UAV and VT is the initial position of the virtual target. � � is the closest interception point reachable simultaneously by the UAV and the virtual target from their initial positions in the most convenient configuration (green). Length � ��� is then the Cartesian distance |�� � |. � � is an interception point reachable simultaneously by the UAV and the virtual target from their initial positions in configuration that requires the longest optimal Dubins curve (red). � ��� (resp.� ��� ) is a distance, which the virtual target travels in the time that is necessary for the UAV to fly distance � ��� (resp.� ��� ). To be able to use CSC manoeuvres and be sure that they are optimal, we need the whole line segment P � P � of the virtual target trajectory to be outside of the circle with radius 4� ��� (see Figure 14). To find the points � � and � � , we need to estimate the distances � ��� and � ��� .
We define parameter � as the wind speed to air speed ratio. We can express the relation between distance d and length � as: Because � ��� is the shortest distance to intercept the virtual target in its distance � ��� from the initial point and because the point � � at � ��� must be outside of the circle with radius 4� ��� , it must hold that: For the edge of the irregular zone: Similarly |�� � | � �� ��� (equality holds on the edge of the irregular zone). However, in this case the path � ��� between S and � � corresponds to the longest optimal Dubins curve and from the Figure 13 we can see that: � ��� � �� ��� � ��� ��� The irregular zone is then the corridor between the two translations of the �� ��� circle. The translations are specified by vectors of size � ��� and � ��� in the direction of the wind (see the grey zone in Figure 14).
We are aware of the fact, that the irregular zone subsumes the worst-case situations and could be made smaller. We are planning to closely specify it in future work.


My First Deep Learning System of 1991 + Deep Learning Timeline 1962-2013

Deep Learning has attracted significant attention in recent years. Here I present a brief overview of my first Deep Learner of 1991, and its historic context, with a timeline of Deep Learning highlights.

: My first Deep Learning system of 1991 used a deep stack of recurrent neural networks (a Neural Hierarchical Temporal Memory) pre-trained in unsupervised fashion to accelerate subsequent supervised learning [79,81,82].
In 2009, our Deep Learning Artificial Neural Networks became the first Deep Learners to win official international pattern recognition competitions [40,83] (with secret test sets known only to the organisers); by 2012 they had won eight of them (TL 1.13), including the first contests on object detection in large images (ICPR 2012) [2,16] and image segmentation (ISBI 2012) [3,15]. In 2011, they achieved the world's first superhuman visual pattern recognition results [20,19]. Others have implemented very similar techniques, e.g., [51], and won additional contests or set benchmark records since 2012, e.g., (TL 1.13, TL 1.14). The field of Deep Learning research is far older though-compare the timeline (TL) further down.
My first Deep Learner dates back to 1991 [79,81,82] (TL 1.7). It can perform credit assignment across hundreds of nonlinear operators or neural layers, by using unsupervised pre-training for a stack of recurrent neural networks (RNN) (deep by nature) as in Figure 1. Such RNN are general computers more powerful than normal feedforward NN, and can encode entire sequences of input vectors.
The basic idea is still relevant today. Each RNN is trained for a while in unsupervised fashion to predict its next input. From then on, only unexpected inputs (errors) convey new information and get fed to the next higher RNN which thus ticks on a slower, self-organising time scale. It can easily be shown that no information gets lost. It just gets compressed (much of machine learning is essentially about compression). We get less and less redundant input sequence encodings in deeper and deeper levels of this hierarchical temporal memory, which compresses data in both space (like feedforward NN) and time. There also is a continuous variant [84].
One ancient illustrative Deep Learning experiment of 1993 [82] required credit assignment across 1200 time steps, or through 1200 subsequent nonlinear virtual layers. The top level code of the initially unsupervised RNN stack, however, got so compact that (previously infeasible) sequence classification through additional supervised learning became possible.
There is a way of compressing higher levels down into lower levels, thus partially collapsing the hierarchical temporal memory. The trick is to retrain lower-level RNN to continually imitate (predict) the hidden units of already trained, slower, higher-level RNN, through additional predictive output neurons [81,79,82]. This helps the lower RNN to develop appropriate, rarely changing memories that may bridge very long time lags.
The Deep Learner of 1991 was a first way of overcoming the Fundamental Deep Learning Problem identified and analysed in 1991 by my very first student (now professor) Sepp Hochreiter (TL 1.6): the problem of vanishing or exploding gradients [46,47,10]. The latter motivated all our subsequent Deep Learning research of the 1990s and 2000s.
Particularly successful are stacks of LSTM RNN [40] trained by Connectionist Temporal Classification (CTC) [36]. On faster computers of 2009, this became the first RNN system ever to win an official international pattern recognition competition [40,83], through the work of my PhD student and postdoc Alex Graves, e.g., [40]. To my knowledge, this also was the first Deep Learning system ever (recurrent or not) to win such a contest (TL 1.10). (In fact, it won three different ICDAR 2009 contests on connected handwriting in three different languages, e.g., [83,40], TL 1.10.) A while ago, Alex moved on to Geoffrey Hinton's lab (Univ. Toronto), where a stack [40] of our bidirectional LSTM RNN [39] also broke a famous TIMIT speech recognition record [38] (TL 1.14), despite thousands of man years previously spent on HMM-based speech recognition research. CTC-LSTM also helped to score first at NIST's OpenHaRT2013 evaluation [11].
Our (multi-column [19]) GPU-MPCNN [18] (TL 1.12) were adopted by the groups of Univ. Toronto / Stanford / Google, e.g., [51,25] (TL 1.13, TL 1.14). Apple Inc., the most profitable smartphone maker, hired Ueli Meier, member of our Deep Learning team that won the ICDAR 2011 Chinese handwriting contest [22]. ArcelorMittal, the world's top steel producer, is using our methods for material defect detection and classification, e.g., [63]. Other users include a leading automotive supplier, recent start-ups such as deepmind (which hired four of my former PhD students/postdocs), and many other companies and leading research labs. One of the most important applications of our techniques is biomedical imaging [16] (TL 1.13, TL 1.14), e.g., for cancer prognosis or plaque detection in CT heart scans.
In particular, most competition-winning or benchmark record-setting Deep Learners (TL 1.10 -TL 1.14) use one of two supervised techniques developed in my lab: (1) recurrent LSTM (1997) (TL 1.8) trained by CTC (2006) [36], or (2) feedforward GPU-MPCNN (2011) [18] (TL 1.12) (building on earlier work since the 1960s mentioned in the text above). Nevertheless, in many applications it can still be advantageous to combine the best of both worlds -supervised learning and unsupervised pre-training, like in my 1991 system described above [79,81,82].

1991: Fundamental Deep Learning Problem
By the early 1990s, experiments had shown that deep feedforward or recurrent networks are hard to train by backpropagation (TL 1.2). My student Hochreiter [46] discovered and analyzed the reason, namely, the Fundamental Deep Learning Problem due to vanishing or exploding gradients. Compare [47].

1991: Deep Hierarchy of Recurrent NN
My first recurrent Deep Learning system (present paper) partially overcame the fundamental problem (TL 1.6) through a deep RNN stack pre-trained in unsupervised fashion [79,81,82] to accelerate subsequent supervised learning. This was a working Deep Learner in the modern post-2000 sense, and also the first Neural Hierarchical Temporal Memory.

2010: Plain Backpropagation on GPUs Yields Excellent Results
In 2010, a new MNIST record was set by good old backpropagation (TL 1.2) in deep but otherwise standard NN, without unsupervised pre-training, and without convolution (but with training pattern deformations). This was made possible mainly by boosting computing power through a fast GPU implementation [17]. (A year later, first human-competitive performance on MNIST was achieved by a deep MCMPCNN (TL 1.12) [22].)

2011: MPCNN on GPU / First Superhuman Visual Pattern Recognition
In 2011, Ciresan et al. introduced supervised GPU-based Max-Pooling CNN or Convnets (MPCNN) [18], today used by most if not all feedforward competition-winning deep NN (TL 1.13, TL 1.14). The first superhuman visual pattern recognition in a controlled competition (traffic signs [87]) was achieved [20,19] (twice better than humans, three times better than the closest artificial NN competitor, six times better than the best non-neural method), through deep and wide Multi-Column (MC) GPU-MPCNN [18,19], the current gold standard for deep feedforward NN.

2012: First Contests Won on Object Detection and Image Segmentation
2012 saw the first Deep learning system (a GPU-MCMPCNN [18,19], TL 1.12) to win a contest on visual object detection in large images (as opposed to mere recognition/classification): the ICPR 2012 Contest on Mitosis Detection in Breast Cancer Histological Images [2,74,16]. An MC (TL 1.12) variant of a GPU-MPCNN also achieved best results on the ImageNet classification benchmark [51]. 2012 also saw the first pure image segmentation contest won by Deep Learning (again through a GPU-MCMPCNN), namely, the ISBI 2012 Challenge on segmenting neuronal structures [3,15]. This was the 8th international pattern recognition contest won by my team since 2009 [1].

2013: More Contests and Benchmark Records
In 2013, a new TIMIT phoneme recognition record was set by deep LSTM RNN [38] (TL 1.8, TL 1.10). A new record [24] on the ICDAR Chinese handwriting recognition benchmark (over 3700 classes) was set on a desktop machine by a GPU-MCMPCNN with almost human performance. The MICCAI 2013 Grand Challenge on Mitosis Detection was won by a GPU-MCMPCNN [88,16]. Deep GPU-MPCNN [18] also helped to achieve new best results on ImageNet classification [95] and PASCAL object detection [34].
Additional contests are mentioned in the web pages of the Swiss AI Lab IDSIA, the University of Toronto, NY University, and the University of Montreal.


extfieldsofstudy: Mathematics Medicine Physics

text:
Minimum vertex cover problems on random hypergraphs: replica symmetric solution and a leaf removal algorithm

We study minimum vertex cover problems on random \alpha-uniform hypergraphs using two different approaches, a replica method in statistical mechanics of random systems and a leaf removal algorithm. It is found that there exists a phase transition at the critical average degree e/(\alpha-1). Below the critical degree, a replica symmetric ansatz in the statistical-mechanical method holdsand the algorithm estimates a solution of the problem which coincide with that by the replica method. In contrast, above the critical degree, the replica symmetric solution becomes unstable and these methods fail to estimate the exact solution.These results strongly suggest a close relation between the replica symmetry and the performance of approximation algorithm.

The more crucial part of everyday life computers bear, the more significance computer science and information theory seem to have. In particular, the computational complexity theory shows the difficulty, the limit of improving algorithms, to solve theoretical computational problems. It has revealed that the problems belong to several classes such as P and NP and there are many inclusion relations between these classes. For example, 2-satisfiability problems (2-SAT) belong to a class of P guaranteed to be solved in polynomial time. 3-SAT and the vertex cover problems belong to a class of NPcomplete [1]. These problems are deeply related to the well-known P versus NP problem plaguing the theoretical computer scientists, who have studied the worst-case performance to solve the computational problems. Among many types of combinatorial optimization problems, the minimum vertex cover problem (min-VC) belongs to a class of NP-hard. The approximation algorithm for the min-VC and its performance have been studied [2]. The application of the problem is to search a file on a file storage [3] and to improve the group testing [4].
In addition to the worst-case analysis, an important alternative is the study of typical-case behavior on a class of random instances of the computational problems. Recently, statistical-mechanical methods of random spin systems have been applied to the problems such as K-SAT and constraint-satisfaction problems [5]. These methods, developed in the spin-glass theory [6], enable us to study the typical properties of the randomized problems. For example, the statistical-mechanical approaches find a SAT/UNSAT transition of K-SAT [7], p-XOR-SAT [8], q-coloring [9] and min-VC [10][11][12][13]. These results clarify that there is a so-called replica symmetric (RS) phase where a replica symmetry ansatz provides correct estimates of the typical properties, and a replica symmetry breaking (RSB) phase where those estimates become unstable. Together with these approaches, a typical-case performance of some approximation algo-rithms has been also studied [14][15][16], suggesting that there is a non-trivial relation between the replica symmetry and the performance of approximation algorithms.
In this Letter, we study the minimum vertex cover problem on a random hypergraph. The random graph is defined by two distributions, the degree distribution and the edge size distribution. The degree means the number of edges connecting to a vertex and the edge size represents the number of vertices connected to an edge. As the former distribution, the Poisson distribution and the delta function are often used and they are called an Erödos-Rényi random graph and a regular random graph, respectively [5]. As the latter distribution, one uses the delta function with a mean α, which yields a random graph with the same edge size as α called a random αuniform hypergraph. In general, a statistical-mechanical model defined on a hypergraph has multi-body interactions determined by its edge size. In contrast to a conventional two-body interaction, the higher-order multibody interactions often change a type of phase transition and a breaking pattern of the replica symmetry as shown in the p-body spin glass model [17]. From this viewpoint, influence of an edge size on the typical estimates of random computational problems is investigated by statistical-mechanical approaches. In fact, it has been revealed that the edge size changes the properties of some problems such as K-SAT [5,7], q-coloring [18] and min-VCs on K-uniform regular random hypergraphs [19]. It is also found that there exists a P/NP transition between 2-SAT and 3-SAT [20]. Here we study the typical case of the size of the min-VC, explained later, on random αuniform hypergraphs and focus on the relation between the replica symmetry and the performance of an approximation algorithm called a leaf removal algorithm.
Let us suppose that an α-uniform hypergraph G = (HV, HE) consists of N vertices i ∈ HV = {1, · · · , N } and (hyper)edges (i 1 , · · · , i α ) ∈ HE ⊂ HV α (i 1 < · · · < i α ). We define covered vertices as a subset HV ′ ⊂ HV and covered edges as a subset of edges connecting to at least a covered vertex. The vertex cover problem on the hypergraph G is to find a set of the covered vertices HV ′ by which all edges are covered. We define the cover ratio on G as |HV ′ |/N with |HV ′ | being the size of the vertex cover problem. The min-VC on G is to search a set of the covered vertices with the minimum cover ratio. In the random α-uniform hypergraph all the edges are set independently from all α-tuples of vertices with probability p. The degree distribution of the graph converges to the Poisson distribution with the average degree c, which is given as c = pN α−1 /(α − 1)! for large N . In this Letter, we focus on an average of the minimum cover ratio x c over the sparse random hypergraphs with the average degree c being O(1).
The vertex cover problems are mapped on the lattice gas model [10,11,21] on the random hypergraphs. We define a variable ν i on each vertex, representing the existence of a gas particle, which takes 0 if a vertex i is covered and 1 if uncovered. An covered edge has at least a vertex with ν i = 0 in its connecting vertices. Thus, an indicator function for a given particle configuration which takes 1 if ν is a solutions of the vertex cover problem on the hypergraph, and 0 otherwise. Using the indicator function, the grand canonical partition function of the model reads where µ is a chemical potential and the sum is over all configurations of ν. In this formulation, only the solutions of the vertex cover problem contribute the partition function and its ground states in a large µ limit are given by the solutions of the min-VC. To study the typical case of min-VCs we need to take the average over the random hypergraphs and the limit as N → ∞. Then, the average minimum-cover ratio is represented as where · · · µ is the grand canonical average and E is the average over the random hypergraph ensemble. Our aim is to obtain the theoretical estimate of the average minimum-cover ratio as a function of the average degree c.
The average minimum-cover ratio is derived from the averaged grand potential density −(µN ) −1 E ln Ξ, which is obtained by using the replica method for finite connectivity graphs [22]. Following the standard procedure of the replica method, the original problem is reduced to solving a saddle-point equation of a replicated order parameter functional. To proceed the calculation, we assume the RS ansatz that the solution of the saddle-point equation has a replica symmetry. Introducing a local field on a vertex associated to the order parameter and its distribution function, we obtain the saddle-point equation of the distribution. Finally, under the RS ansatz, the average minimum-cover ratio is obtained as a function of the average degree c, We call this estimate the RS solution of min-VCs. This solution is also obtained by an alternative cavity method [12]. Although the instability of the RS solution such as the de Almeida-Thouless instability [23] must be examined to validate the solution, we here naively study an instability condition of the saddlepoint equation against a perturbation of the local field distribution within the RS sector. The analysis leads to a critical value of the average degree c * = e/(α−1) above which the RS solution becomes unstable. These results, x c and c * , include the case of α = 2 [10]. The obtained x c gives a correct value below the critical average degree, while a RSB solution for x c is required above it.
Here we turn our attention to the estimate of x c by using an approximation algorithm. The leaf removal algorithm has been proposed as an approximation algorithm to solve a min-VC on a graph with α = 2 [24] and has also been applied to search for a k-core [25] and a 3-XOR-SAT solution [15]. For a min-VC on a given graph, this algorithm consists of iterative steps, where vertices called a leaf, as well as the edges connecting to the leaves, are removed from the graph with covered vertices appropriately assigned to those vertices. This removal step makes new leaves and the algorithm continues in an iterative way until the leaf is empty. By this procedure, the minimum cover ratio is estimated correctly at least for the removed part of the graph. We consider the global leaf removal (GLR) algorithm [14], which removes simultaneously all the leaves found in a recursive step. We focus on the expansion of this algorithm for the min-VC on a hypergraph with α = 3, while it is straightforward to extend it to that on a hypergraph with α ≥ 4. A crucial point in our algorithm is in definition of leaf, where a leaf {i, j, k} ∈ HV 3 (i < j < k) is defined as a 3-tuple of vertices connecting to an edge (i, j, k), at least two of which the degree is one. The definition of the GLR algorithm is as follows: Step 1: The initial graph G is named G (0) . Set k = 0.
Step 2: Search all leaves from the graph G (k) . If there is no leaf, go to Step 6.
Step 3: Remove all the leaves except for the vertices which belong to more than two leaves, named bunch of leaves [14], and remove only one of leaves in each bunch.
Step 4: Assign covered vertices to the one with the maximal degree in each removed leaf from G (k) .
Step 5: The left graph is named G (k+1) , and return to Step 2 with k incresed by one.
Step 6: If there exist connected vertices in the left graph, assign all of them to covered vertices. Stop the algorithm.
It is proven that the result of the algorithm is independent of order of removal and a selection of a leaf out of a bunch of leaves in the removal process. When the recursive steps stop, the left graph consists of isolated vertices and a core, which is defined as a set of vertices connecting to edges without leaves. Vertices in a bunch of leaves which are not selected for the removal in Step 3 become isolated and the core of the order O(N ) exists in large c. We note that Step 4 can be omitted if one is interested only in the minimum cover ratio, not the covered vertices. Because the algorithm covers all vertices in the core without searching the solution of the min-VC as shown in Step 6, the existence of the core of the order O(N ) leads to overestimation of the average minimum-cover ratio. We study the core size at the end of the GLR algorithm by numerically performing the above-mentioned procedure for finite-size random hypergraphs with α = 3. While the computational time for the GLR algorithm is proportional to the number of vertices, it takes time of the order O(N 3 ) for generating a random graph. To avoid it, we use the microcanonical ensemble [14] with fixing the number of edges to the expectation number of edges cN/3, ignoring fluctuation of the average degree. We expect that such fluctuation is irrelevant in a large size N limit. In Fig. 1, the core size density obtained by numerical simulations is presented as a function of the average degree c up to the size N = 10 5 . The data averaged over 10 4 random graphs converges well for large sizes and a giant core with O(N ) emerges above a certain value of c. We discuss the asymptotic behavior of the recursive procedure in the GLR algorithm. We introduce the average fraction of the core c n and the isolated vertices i n over random hypergraphs after n-th step of the algorithm, and find i n = e 2n+1 + 2e 2n + 2ce 2n e 2 2n−1 − 2, c n = e 2n − e 2n+1 − 2ce 2n e 2 2n−1 + 2ce 3 2n−1 , where a parameter e n obeys a recursion relation e n = exp(−ce 2 n−1 ) with the initial condition e −1 = 0. A detailed derivation of the formulas will be reported in a separate paper [26]. By definition, the average fraction Below the critical value, the sequence {e n } converges to the unique value [W (2c)/(2c)] 1/2 and consequently the core size c ∞ is zero. Above the critical value, however, a bifurcation occurs in the recursion relation and the sequence has a cycle with period two. This type of the transition would occur above α = 3 at the critical average degree c * = e/(α − 1). Because e −1 = 0, an even term e 2n is larger than that at one-step later, that is e 2n+1 . We compute the limiting values lim n→∞ e 2n+1 and lim n→∞ e 2n numerically as a function of c. The difference between them yields emergence of the core of the order of O(N ). We present the core size density obtained from the asymptotic analysis of the recursion relation by the solid line in Fig. 1, which coincides with the data by numerical simulations. Thus, we confirm that a core percolation occurs at the critical average degree in the GLR algorithm, which coincides with that of the RS instability. From the analysis near the critical degree, it is found that the size of the core emerges linearly near above the critical average degree. These findings, the bifurcation in the recursion relation and the core percolation, are common in the min-VCs on random graphs with α = 2. As mentioned above, the GLR algorithm estimates the minimum cover ratio by the size of the removed part in the graph during the recursive procedure, which is given as r ∞ = 1 − i ∞ − c ∞ . Taking one-third of r ∞ and adding c ∞ to the value, we obtain the estimate of the average minimum-cover ratio by the algorithm. Thus, we find that below the critical average degree e/2 the estimate r ∞ /3 coincides with the RS solution Eq. (4) estimated by the replica method. In contrast, the sequence {e n } of the algorithm does not converge to a unique value above the critical value and the GLR algorithm could not give a precise estimate of x c there.
In order to confirm whether these analyses estimate the average minimum-cover ratio x c correctly, we also evaluate the min-VCs by the Markov chain Monte Carlo method. We use the replica exchange Monte Carlo method (EMC) [27], for accelerating the dynamics of the system, with 50 replicas in the range of the chemical potential from −2 to 10. In our Monte Carlo simulations, the smallest cover ratio found in typically 2 17 Monte Carlo steps is used as the estimate of x c for each random graph, which is averaged over 800 hypergraphs randomly generated. The number of vertices of the graph is up to N = 512. The average minimum-cover ratio is extrapolated from these numerical results for finite N . Fig. 2 shows the obtained minimum cover ratio as a function of the average degree c. Below the critical average degree e/2 where the RS solution is considered to be correct, we observe that the MC result is consistent with those by the two approaches, the replica method and the GLR algorithm. Above the critical value, on the other hand, the MC estimate stays slightly above that by the replica method and considerably deviates from that by the GLR algorithm. The former is due to the instability of the RS solution and the latter is the existence of the core of the order O(N ).
To summarize, we consider the minimum vertex cover problems on random α-uniform hypergraphs, and analyze them by the statistical-mechanical method and the approximation algorithm. The replica method estimates the average minimum-cover ratio x c as a function of the average degree c under the replica symmetric assumption. We find that there is an RS/RSB phase transition at the critical average degree c * = e/(α−1), which is well above a percolation threshold c = 1/(α−1) in the random graph. We also perform the global leaf removal algorithm and study the asymptotic behavior of the recursive procedure of the algorithm, particularly in the case of α = 3. If the average degree is below the critical value which coincides with that in the replica theory, there is a core of the order O(1) in the remaining part of the graph, which does not affect the estimate of the minimum cover ratio. In contrast, above the critical value, the core of the order O(N ) emerges, leading to a wrong estimation of the minimum cover ratio. Comparing the results obtained by MC simulations, we confirm that these estimates are correct below the critical average degree, but this is not the case above the critical degree. These results strongly suggest that there is a close relation between the replica symmetry in statistical physics and the performance of the leaf removal algorithm even when the edge size α is larger than two. It is noted that this relation is not always true for all types of random graphs. For instance, the GLR algorithm removes no vertex on regular random graphs with c ≥ 2 because no leaf is found there while, from the point of the statistical-mechanical view, the min-VCs on regular random 2-uniform graphs with degree 2 is described by the RS solution [19]. Thus, the relation depends on a type of random graphs and approximation algorithms. In addition to the leaf removal algorithm, a recent work for the min-VC problem with α = 2 [16] suggests that linear programming algorithms, which are one of the most commonly used tools for solving optimization problems, have the relation discussed in the present work. Further study will need to establish the relation between the replica symmetry and the performance of numerous algorithms. This research was supported by a Grants-in-Aid for Scientific Research from the MEXT, Japan, No. 22340109.


Fast computation of the performance evaluation of biometric systems: application to multibiometric

The performance evaluation of biometric systems is a crucial step when designing and evaluating such systems. The evaluation process uses the Equal Error Rate (EER) metric proposed by the International Organization for Standardization (ISO/IEC). The EER metric is a powerful metric which allows easily comparing and evaluating biometric systems. However, the computation time of the EER is, most of the time, very intensive. In this paper, we propose a fast method which computes an approximated value of the EER. We illustrate the benefit of the proposed method on two applications: the computing of non parametric confidence intervals and the use of genetic algorithms to compute the parameters of fusion functions. Experimental results show the superiority of the proposed EER approximation method in term of computing time, and the interest of its use to reduce the learning of parameters with genetic algorithms. The proposed method opens new perspectives for the development of secure multibiometrics systems by speeding up their computation time.

Introduction
Biometrics [1] is a technology allowing to recognize people through various personal factors. It is an active research field which design new biometric traits from time to time (like finger knuckle recognition [2]). We can classify the various biometric modalities among three main families: • Biological : the recognition is based on the analysis of biological data linked to an individual (e.g, DNA, EEG analysis, ...).
• Behavioural : the recognition is based on the analysis of the behaviour of an individual while he is performing a specific task (e.g, signature dynamics, gait, ...).
• Morphological : the recognition is based on the recognition of different physical patterns, which are, in general, permanent and unique (e.g, fingerprint, face recognition, ...).
It is mandatory to evaluate these biometric systems in order to quantify their performance and compare them.
These biometric systems must be evaluated in order to compare them, or to quantify their performance. To evaluate a biometric system, a database must be acquired (or a common public dataset must be used). This database must contain as many users as possible to provide a large number of captures of their biometric data. These data are separated into two different sets: • the learning set which serves to compute the biometric reference of each user • the validating set which serves to compute their performance.
When comparing test samples to biometric references, we obtain two different kinds of scores: • the intrascores represent comparison scores between the biometric reference (computed thanks to the learning set) of an individual and biometric query samples (contained in the validating set) • the interscores represent comparison scores between the biometric reference of an individual and the biometric query samples of the other individuals.
From these two sets of scores, we can compute various error rates, from which the EER is one functioning point which represents a very interesting error rate often used to compare biometric systems. In order to have reliable results, it is necessary to evaluate the performance of biometric system with huge datasets. These huge datasets produce numbers of scores. As the time to evaluate the performance of a biometric system depends on the quantity of available scores, we can see that evaluation may become very long on these large datasets. In this paper, we present a very fast way to compute this error rate, as well as its confidence interval in a non parametric way, on different datasets of the literature.
Nevertheless, there will always be users for which one modality (or method applied to this modality) will give bad results. These low performances can be implied by different facts: the quality of the capture, the acquisition conditions, or the individual itself. Biometric multi-modality (or multibiometrics) allows to compensate this problem while obtaining better biometric performances (i.e., better security by accepting less impostors, and better usability by rejecting less genuine users) by expecting that the errors of the different modalities are not correlated. So, the aim of multibiometrics is to protect logical or physical access to a resource by using different biometric captures. We can find different types of multibiometrics systems. Most of them are listed in [3], they use: 1. Different sensors of the same modality (i.e., capacitive or resistive sensors for fingerprint acquisition); 2. Different representations of the same capture (i.e., use of points of interest or texture); 3. Different biometric modalities (i.e., face and fingerprint); 4. Several instances of the same modality (i.e., left and right eye for iris recognition); 5. Multiple captures (i.e., 25 images per second in a video used for face recognition); 6. An hybrid system composed of the association of the previous ones.
In the proposed study, we are interested in the first four kinds of multimodality. We also present in this paper, a new multibiometrics approach using various fusion functions parametrized by genetic algorithms using a fast EER (Equal Error Rate) computing method to speed up the fitness evaluation.
This paper is related to high performance computing, because algorithms are designed to work in an infrastructure managing the biometric authentication of millions of individuals (i.e., border access control, logical acces control to webservices). To improve the recognition rate of biometric systems, it is necessary to regularly update the biometric reference to take into account intra class variability. With the proposed approach, the time taken to update the biometric reference would be lowered. The faster is the proposed method, the more we can launch the updating process (or the more users we can add to the process). We also propose an adaptation of the proposed EER computing method which gives confidence intervals in a non parametric way (i.e., by computing the EER several times through a bootstraping method). The confidence intervals are computed on a single CPU, on several CPUs on the same machine and on several machines.
The main hints of the papers are: • the proposition of a new method to approximate the EER and its confidence interval in a fast way • the proposition of two original functions for multibiometrics fusion.
The plan is organized as following. Section 2 presents the background of the proposed work. Section 3 presents the proposed method for computing the approximated value of the EER and its confidence interval. Section 4 validates them. Section 5 presents the proposed multibiometrics fusion functions and their performance in term of biometric recognition and computation time against the baseline. Section 6 gives perspectives and conclusions of this paper.

Evaluation of Biometric Systems
Despite the obvious advantages of this technology in enhancing and facilitating the authentication process, its proliferation is still not as much as attended [4]. As argued in the previous section, biometric systems present several drawbacks in terms of precision, acceptability, quality and security. Hence, evaluating biometric systems is considered as a challenge in this research field. Nowadays, several works have been done in the literature to evaluate such systems. Evaluating biometric systems is generally realized within three aspects as illustrated in figure 1: usability, data quality and security.

Usability
According to the International Organization for Standardization ISO 13407:1999 [5], usability is defined as "The extent to which a product can be used by specified users to achieve specified goals with effectiveness, efficiency, and satisfaction in a specified context of use".
• Efficiency which means that users must be able to accomplish the tasks easily and in a timely manner. It is generally measured as task time; • Effectiveness which means that users are able to complete the desired tasks without too much effort. It is generally measured by common metrics include completion rate and number of errors such failure-to-enroll rate (FTE) [6]; • User satisfaction which measures users' acceptance and satisfaction regarding the system. It is generally measured by studying several properties such as easiness to use, trust in the system, etc. The acceptability of biometric systems is affected by several factors. According to [7], some members of the human-computer interaction (HCI) community believe that interfaces of security systems do not reflect good thinking in terms of creating a system that is easy to use, while maintaining an acceptable level of security. Existing works [8,9] show also that there is a potential concern about the misuse of personal data (i.e., templates) which is seen as violating users' privacy and civil liberties. Moreover, one of our previous work [10] shows the necessity of taking into account users' acceptance and satisfaction when designing and evaluating biometric systems. More generally speaking, even if the performance of a biometric system outperformed another one, this will not necessarily mean that it will be more operational or acceptable;

Data quality
It measures the quality of the biometric raw data [11,12]. Low quality samples increase the enrollment failure rate, and decrease system performance. Therefore, quality assessment is considered as a crucial factor required in both the enrollment and verification phases. Using quality information, the bad quality samples can be removed during enrollment or rejected during verification. Such information could also be used in soft biometrics or multimodal approaches [13]. Such type of assessment is generally used to quantify biometric sensors, and could be also used to enhance system performance;

Security
It measures the robustness of a biometric system (algorithms, architectures and devices) against attacks. Many works in the literature [14,15,16] show the vulnerabilities of biometric systems which can considerably decrease their security. Hence, the evaluation of biometric systems in terms of security is considered as an important factor to ensure its functionality. The International Organization for Standardization ISO/IEC FCD 19792 [17] addresses the aspects of security evaluation of such systems. The report presents an overview of biometric systems vulnerabilities and provide some recommendations to be taking into account during the evaluation process. Nowadays, only few partial security analysis studies with relation to biometric authentication systems exist. According to ISO/IEC FCD 19792 [17], the security evaluation of biometric systems is generally divided into two complementary assessments: 1) assessment of the biometric system (devices and algorithms) and 2) assessment of the environmental (for example, is the system is used indoor or outdoor?) and operational conditions (for example, tasks done by system administrators to ensure that the claimed identities during enrolment of the users are valid). A type-1 security assessment method is presented in a personal previous work [18]. The proposed method has shown its efficiency in evaluating and comparing biometric systems.

Performance Evaluation of Biometric Systems
The performance evaluation of biometric systems is now carefully considered in biometric research area. We need a reliable evaluation methodology in order to put into obviousness the benefit of a new biometric system. Nowadays, many efforts have been done to achieve this objective. We present in section 2.2.1 an overview of the performance metrics, followed by the research benchmarks in biometrics as an illustration of the evaluation methodologies used in the literature for the comparison of biometric systems.

Performance metrics
By contrast to traditional methods, biometric systems do not provide a cent per cent reliable answer, and it is quite impossible to obtain such a response. The comparison result between the acquired biometric sample and its corresponding stored template is illustrated by a distance score. If the score is lower than the predefined decision threshold, then the system accepts the claimant, otherwise he is rejected. This threshold is defined according to the security level required by the application. Figure 3 illustrates the theoretical distribution of the genuine and impostor scores. This figure shows that errors depend from the used threshold. Hence, it is important to quantify the performance of biometric systems. The International Organization for Standardization ISO/IEC 19795-1 [6] proposes several statistical metrics to characterize the performance of a biometric system such as: • Failure-to-enroll rate (FTE): proportion of the user population for whom the biometric system fails to capture or extract usable information from biometric sample; • Failure-to-acquire rate (FTA): proportion of verification or identification attempts for which a biometric system is unable to capture a sample or locate an image or signal of sufficient quality; • False Acceptation Rate (FAR) and False Rejection Rate (FRR): FAR is the proportion of impostors that are accepted by the biometric system, while the FRR is the proportion of authentic users that are incorrectly denied. The computation of these error rates is based on the comparison of the scores against a threshold (the direction of the comparison is reversed if the scores represent similarities instead of distances). FRR and FAR are respectively computed (in the case of a distance score) as in (1) and (2), where intra i (respectively inter i ) means the intra score at position i in the set of intra score (respectively inter score at position i) and Card(set) is the cardinal of the set in argument, thr is the decision threshold, and 1 is the indicator function.
• Receiver operating characteristic (ROC) curve: the ROC curve is obtained by computing the couple of (FAR, FRR) for each tested threshold. It plots the FRR versus the FAR. The aim of this curve is to present the tradeoff between FAR and FRR and to have a quick overview of the system performance and security for all the parameters configurations.
• Equal error rate (EER): it is the value where both errors rates, FAR and FRR, are equals (i.e., FAR = FRR). It constitutes a good indicator, and the most used, to evaluate and compare biometric systems. In other words, lower the EER value is, higher the accuracy of the system. Using the ROC curve, the EER is computed by selecting the couple of (FAR, FRR) having the smallest absolute difference (3) at the given threshold τ : and returning their average (3): By this way, we have obtained the best approaching EER value with the smallest precision error. The classical EER computing algorithm is presented in the Figure 2 1 . From Figure 2, we can see that the complexity is in O(n * m) with n the number of thresholds held in the computation, and, m the number of scores in the dataset. As it is impossible to reduce m, we have to find a better method which reduces n. We did not find, in the literature, methods allowing to reduce computation time in order to obtain this EER.

Biometrics Datasets
A public dataset allows researchers to test their algorithm and compare them with those from the state of the art. It takes a lot of time and energy to build a large and significant dataset. It is very convenient to download one for research purposes. We present in this section an overview of the used datasets in this paper. Table 1 presents a summary of these datasets.  • Biometric Scores Set -Release 1 (BSSR1) The BSSR1 [19] database is an ensemble of scores sets from different biometric systems. In this study, we are interested in the subset containing the scores of two facial recognition systems and the two scores of a fingerprint recognition system applied to two different fingers for 512 users. This database has been used many times in the literature [20,21].
• BANCA The second database is a subset of scores produced from the BANCA database [22]. The selected scores correspond to the following one labelled: The database as two subsets G1 and G2. G1 set is used as the learning set, while G2 set is used as the validation set.

• PRIVATE
The last database is a chimeric one we have created for this purpose by combining two public biometric template databases: the AR [23] for the facial recognition and the GREYC keystroke [24] for keystroke dynamics [25,26]. The AR database is composed of frontal facial images of 126 individuals under different facial expressions, illumination conditions or occlusions. These images have been taken during two different sessions with 13 captures per session. The GREYC keystroke contains the captures on several sessions on two months of 133 individuals. Users were asked to type the password "greyc laboratory" 6 times on a laptop and 6 times on an USB keyboard by interlacing the typings (one time on a keyboard, one time on another). We have selected the first 100 individual of the AR database and we have associated each of these individuals to another one in a subset of the GREYC keystroke database having 5 sessions of captures. We then used the 10 first captures to create the biometric reference of each user and the 16 others to compute the intra and inter scores. These scores have been computed by using two different methods for the face recognition and two other ones for the keystroke dynamics.

Multibiometrics
We focus in this part on the state of the art on multimodal systems involving biometric modalities usable for all computers (keystroke, face, voice...). The scores fusion is the main process in multimodal systems. It can be operated on the scores provided by algorithms or in the templates themselves [27]. In the first case, it is necessary to normalize the different scores as they may not evolve in the same range. Different methods can be used for doing this, and the most efficient methods are zscore, tanh and minmax [28]. Different kinds of fusion methods have been applied on biometric systems. The fusion can be done with multiple algorithms of the same modality. For example, in [29], three different keystroke dynamics implementations are fused with an improvement of the EER, but less than 40 users are involved in the database. In [30], two keystroke dynamics systems are fused together by using weighted sums for 50 users, but no information on the weight computing is provided. The fusion can also be done within different modalities in order to improve the authentication process. In [31], authors use both face and fingerprint recognition, the impact of error rate reduction is used to reduce the error when adapting the user's biometric reference. There is only one paper (to our knowledge) on keystroke dynamics fusion with another kind of biometric modality (voice recognition): it is presented in [32], but only 10 users are involved in the experiment. In [33], multi-modality is done on fingerprints, speech, and face images on 50 individuals. Fusion has been done with SVM [34] with good improvements, especially, when using user specific classifiers.
Very few multimodal systems have been proposed for being used in classical computers and the published ones have been validated on small databases. In order to contribute to solve this problem, we propose a new approach in the following section.

Fast EER and Confidence Intervals Estimation
We propose a kind of dichotomic EER computing function, in order to quickly approximate its value. Thanks to this computing speed up, we can use it in time consuming applications. Finally, we present a confidence interval computing method based on our approximated EER calculation associated to parallel and distributed computing.

EER Estimation
Computation time to get the EER can be quite important. When the EER value needs to be computed a lot of time, it is necessary to use a faster way than the standard one. In the biometric community, the shape of the ROC curve always follows the same pattern: it is a monotonically decreasing function (when we present FRR against FAR, or increasing when we present 1-FRR against FAR) and the EER value is the curve's point having x ROC = y ROC (or F AR = F RR). Thanks to this fact, the curve symbolising the difference of y ROC against x ROC is also a monotonically decreasing function from 1 to −1, where the point at y DIF F = 0 represents the EER (and its value is x DIF F because x ROC = y ROC or F AR = F RR). With these information, we know that to get the EERs, we need to find the x DIF F for which y DIF F is the closest as possible to zero. An analogy with the classical EER computing, would be to incrementally compute y DIF F for each threshold by increasing order and stop when y DIF F changes of sign. By this way, we can expect to do half thresholds comparisons than with the classical way if scores are correctly distributed. A clever way is to use something equivalent to a divide and conquer algorithm like the binary search and obtain a mean complexity closer to O(log(n)). That is why we have implemented a polytomous version of EER computing: 1. We chose i thresholds linearly distributed on the scores sets 2. For each threshold t among the i thresholds, we compute the FAR and FRR values (F AR t , F RR t ) 3. We take the two following thresholds t1 and t2 having sign(F RR t1 − F AR t1 ) different of sign(F RR t2 − F AR t2 ) 4. We repeat step 2 with selecting i thresholds between t1 and t2 included while F RR t1 − F AR t1 ) does not reach the attended precision.
By this way, the number of threshold comparisons is far smaller than in the classical way. Its complexity analysis is not an easy task because it depends both on the attended precision and the choice of i. It can be estimated as O(log(N )). Figure 4 presents the algorithm while Figure 5 illustrates it by showing the different iterations until getting the EER value with a real world dataset. We have chosen i = 5 points to compute during each iteration. The EER is obtained in five iterations. Circle symbols present the points computed at the actual iteration, triangle symbols present the points computed at the previous iterations, and the dotted curve presents the ROC curve if all the points are computed. Very few points are computed to obtained the EER value. Figure 5f presents the real ROC curve and the ROC curve obtained with the proposed method. We can see that even if we obtain an approximated version of the real ROC curve, it is really similar around the EER value (cross with the dotted lined).   In this case, the EER value is found in five iterations. Each image represents an iteration, with: the real ROC curve, the points computed at the iteration and the points computed at the previous iteration (different thresholds may produce the same points).

Confidence Intervals Estimation
We also provide a method to compute the confidence interval of the EER value. It is based on a bootstrapping method and can be used in a parallelized way.

Bootstrapping
It is interesting to give a confidence interval of an EER value, because we are not totally sure of its value. One way is to obtain this confidence interval parametrically, but it requires to have strong hypothesis on the function of the EER value (the scores come from independant and identically distributed variables, even for the scores of the same user). As such assumption is too strict (and probably false), it is possible to use non parametric confidence intervals. One of these non parametric methods is called "bootstrap" [35]. Such method is often used when the distribution is unknown or the number of samples is too low to correctly estimates the interval. The main aim is to re-sample the scores several times, and, compute the EER value for each of these re-sampling. The boostraping method works as following: 1. Use the intra and inter scores to compute the EERχ. 2. Resample K times the intra and inter scores and store them in intra i and inter i (0 < i <= K).
• Generate the resampled intra i scores by sampling Card(intra) scores with replacement from intra. • Generate the inter i scores by sampling Card(inter) scores with replacement from inter. 3. Compute the K EERs (χ i ) for each couple of intra i and inter i (0 < i <= K). 4. Store the K residuals e i =χ − χ i . 5. The 100(1 − α)% confidence interval of the EER is formed by taking the interval fromχ−e U pper toχ−e Lower , with e Lower and e U pper which respectively represent the α/2th and the 1−α/2th percentiles of the distribution of e. Figure 6 presents the residuals of one run of the bootstraped method on a real world dataset with the lower and upper limits used to compute the confidence interval. We can see that the consuming part of this algorithm is the fact to compute K times the EER. As all the EER computing is totally independent from each other, the K computations can be done in a parallel way.

Parallelization
We propose three different ways to compute the confidence interval: • Single. The single version consists in doing all the computations in a sequential manner (cf. figure 7a). The EER with the original scores is computed. The whole set of resampled scores is created. The EER of each new distribution is computed. The confidence interval is computed from the results. Figure 6: Histogram of residuals and α/2 and 1−α/2 percentiles for a confidence interval at 90%.
• Parallel. The parallel version consists in using the several cores or processors on the computer used for the computation (cf. figure 7b). In this case, several EERs may be computed at the same time (in the better case, with a computer having n processing units, we can compute n different results at the same time). The procedure is the following: the EER with the original scores is computed. The whole set of resampled scores is created. The main program distributes the K EERs computations on the different processing units. Each processing unit computes an EER and returns the results, until the main program stops to send it new data. The results are merged together. The confidence interval is computed from the results.
• Distributed. The distributed version consists in using several computers to improve the computation (cf. figure 7c). The computation is done in a parallelized way on each computer. In this case, much more EERs can be computed at the same time. The main program generates a set S = {S 1 , ...S T } of T values symbolising T subworks, were the subwork T i must compute S i EERs (thus S = K). The procedure is the following: the EER with the original scores is computed. The main program sends the intra and inter scores on each worker (i.e., a computer). The main program distributes the T numbers to each worker. Each time a work receive such number (S i ), it computes the S i resamples sets. Then, it computes (using the Parallellized way) the S i EERs by distributing them on its processing units. It merges the S i EERs together and send them to the main program. The main program merges all the results together. The confidence interval is computed from the results.
We can see that the three schemes are totally different. The parallelized version may be used on all recent computers which have several processing units, while the distributed version needs to use several computers connected through a network.

Databases Sets
In order to do these evaluations, we have used three different biometric databases presented in section 2.2.2.

Evaluation of the EER Computing
The two different algorithms for EER computing have been run on five different sets of scores (three of keystroke dynamics and two of face recognition, generated with the PRIVATE database) with various parameters. We call classic the classical way of computing the EER and polyto the proposed version of the algorithm. The classic way is tested by using 50, 100, 500 and 1000 steps to compute the EER. The polytomous way is tested by using between 3 and 7 steps and a precision of 0.01, 0.005 and 0.003. The aim of these tests is to compare how the proposed method performs better than the classical one, and what are its best parameters.

Utility for Non-parametric Confidence Interval
Confidence intervals are also an interesting information on the performance of a system. The properties of the score distribution may forbid the use of parametric methods to compute it. This is where the boostrap method helps us by computing several times the EER with resampling methods. We have tested the computation time of confidence interval computing for two systems of each database, which give us six different systems. We have computed the EER value under three different ways: the polytomous one, the classical one (with 1000 steps), and another we called whole. The whole method is similar to the classical one, except that it uses all the possible scores present in the intra inter scores arrays as thresholds, instead of artificially generating them in a predefined interval (thus, there may be far more thresholds than in the classic method, or far less, depending on the number of available scores in the database). The results are then validated with confidence intervals. The test scripts were written with the Python language. The EER computing methods and the resampling methods have been compiled in machine code thanks to the Cython [36] framework (which speed up computation time). The parallelization is done with the joblib [37] library. The distributed version is done by using the task management provided by Ipython [38]. The standard and parallelized versions have been launched on a recent Linux machine with an Intel R Core TM i5 with 4 cores at 3.20GHz and 4Gb of RAM. Four processes are launched at the same time. For the distributed machine, the orchestrator is the same machine. It distributes the jobs on three multicores machines (itself, another similar machine and an another machine with an Intel R Xeon R with 8 cores at 2.27GHz and 4Gb of RAM. The controller sends two more jobs last machine machine). Table 2 presents the results obtained within the first tested biometric system. We present the name of the method, the error of precision while computing the EER, the computation time in milliseconds and the number of comparisons involved (each comparison corresponds to the comparison of a threshold against the whole set of intra and inter scores). The real computation time taken by a comparison is given in (4), where n is the number of thresholds to compare, A is the timing to do a comparison and B and C depends on the algorithm.

Evaluation of the EER Computing
We can see that the computation time is highly related to the number of comparisons and the size of the score set. Using the Kruskal-Wallis test at a confidence degree equals to 95%, the proposed method significantly outperformed the classical method in terms of errors (with a p value = 0.0305) and computation time (with a p value = 0.002562). The obtained results are slightly similar for the five tested biometrics modalities. We can observe that, in the classic method, using 50 steps gives not enough precise results, while using 1000 gives a very good precision, but is really time consuming; depending on the dataset, 500 steps seems to be a good compromise between precision error and computation time.
In all the polytomous configurations, the computation time is far better than the fastest classic method (50 steps) while having a greatest precision. This precision is always better than the classic method with 100 steps and approach or is better than the precision in 1000 steps. This gain of time is due to the lowest number of involved comparisons. In an n steps classical computing, we need to check n thresholds, while in the polytomous way this number depends both on the dataset and the required precision: with our dataset, it can vary from 8 to 35 which is always lower than 50. As the computation time depends only on this value, we can say that the fastest algorithms are the one having the smallest number of tests to complete.
Based on the number of comparisons (and the timing computation), Table 3 presents the best results for each modality (when several methods return the same number of iterations, the most precise is chosen).
We can argue that the proposed method is better, both in terms of speed and precision error, than the classical way of computing. Based on the results of our dataset, the configuration using 3 steps and a precision of 0.010 seems to be the best compromise between speed and precision. We can now argue that the proposed EER computation will speed up genetic algorithms using the EER as fitness function.

Utility for Non-parametric Confidence Interval
All the methods under all the implementations give similar confidence intervals. We do not discuss on this point, because we are only interested in computation time. Table 4 gives for each method, under each implementation  the mean value of the computation time for all the six different biometric systems. We can observe that, in average, the polytomous version seems far more faster than the other methods (classic and whole). The distributed implementation seems also more faster than the other implementations (Parallel, Single). Using the Kruskal-Wallis test at a confidence degree equals to 95%, the computation time of the proposed method is significantly faster than both classical (p value = 0.00651) and whole (p value = 0.004407) schemes. There was no significant difference of computation time between classical and whole schemes (p value = 0.8002).

Discussion
We have demonstrated the superiority of our EER estimation method against the classical method concerning the computation time. However, the method stops when the required precision is obtained. As the method is iterative, it is not parallelizable when we want only a simple EER. However using a grid computing method greatly improves the computation of confidence intervals.

Application to Multibiometrics Fusion Function Configuration
We propose a biometric fusion system based on the generation of a fusion function parametrized by a genetic algorithm and a fast method to compute the EER (which is used as fitness function) in order to speed the computing time of the genetic algorithm.

Method
We have tested three different kinds of score fusion methods which parameters are automatically set by genetic algorithms [39]. These functions are presented in (5), (6) and (7) where n is the number of available scores (i.e., the number of biometric systems involved in the fusion process), w i the multiplication weight of score i, s i the score i and x i the weight of exponent of score i. (5) is the commonly used weighted sum (note that in this version, the sum of the weights is not equal to 1), while the two others, to our knowledge, have never been used in multibiometrics. We have empirically designed them in order to give more weights to higher scores.
The aim of the genetic algorithm is to optimize the parameters of each function in order to obtain the best fusion function. Each parameter (the w i and x i ) is stored in a chromosome of real numbers. The fitness function is the same for the three genetic algorithms. It is processed in two steps: • fusion: The generated function (5), (6) or (7) are evaluated on the whole set of scores; • error computing: The EER is computed on the result of the fusion. We use the polytomous version of computing in order to highly speed up the total computation time. Table 5 presents the parameters of the genetic algorithms. The genetic algorithms have been trained on a learning set composed of half of the intrascores and half of the interscores of a database and they have been verified with a validation set composed of the others scores. The three databases have been used separately.

Design of Fusion Functions
The generated functions are compared to three methods of the state of the art: sum, mul and min, they have been explored in [28,40]. Table 6 presents, for each database, the EER value of each of its biometric method (noted sn for method n), as well as the performance of the fusion functions of the state of the art. We can see that biometric methods from PRIVATE have more biometric verification errors than the ones of the other databases. Using the Kruskal-Wallis test, the sum (p value = 0.0038) and mul (p value = 0.0038) operators outperformed the min operator. There was no significant difference (p value = 0.935) between both operators sum and mul operators.

Magnitude Of the Gain in Computation Time
We also want to prove that using the proposed EER computation method improves the computation time of the genetic algorithm run. To do that, the previously described process has been repeated two times: • Using the proposed EER computing method with the following configuration: 5 steps and stop at a precision of 0.01.
• Using the classical EER computing method with 100 steps.
The total computation time is saved in order to compare the speed of the two systems. These tests have been done on a Pentium IV machine with 512 Mo of RAM with the Matlab programming language.

Design of Fusion Functions
The EER of each generated function of each database is presented in Table 8 for the learning and validation sets, while Figure 10 presents there ROC  curve on the validation set. We can see that the proposed generated functions are all globally better than the ones from the state of the art (by comparison with Table 6) and the obtained EER is always better than the ones of the sum and mul. The two new fusion functions ( (6) and (7)) give similar or better results than the weighted sum (5). We do not observe over-fitting problems: the results are promising both on the learning and validation sets. We also could expect to obtain even better performance by using more individuals or more generations in the genetic algorithm process, but, in this case, timing computation would become too much important. Their will always be a tradeoff between security (biometric performance) and computation speed (genetic algorithm performance). By the way, the best individuals were provided in the As a conclusion of this part, we increased the performance of multibiometrics systems given the state of the art by reducing errors of 58% for BANCA, 45% for BSSR1 and 22% for PRIVATE. Table 8 presents a summary of the performance of the generated methods both in term of EER and timing computing improvement. The column gain presents the improvement of timing computation between the proposed EER polytomous computation time and the classical one in 100 steps.

Magnitude Of the Gain in Computation Time
We can observe that, in all the cases, the proposed computation methods outperform the classical one (which is not its slowest version). We can see that this improvement depends both on the cardinal of the set of scores and the function to evaluate: there are better improvements for (5). The best gain is about 78% while the smallest is about 28%.

Discussion
Once again, we can observe the interest of our EER estimation method which allows to obtain results far quickly. We can note that the generated function all generate a monotonically decreasing ROC curve which allows to use our method. If the ROC curve does not present this shape, we would be unable to Figure 9: ROC Curve Of The Generated Multibiometrics Fusion Functions on the Validation Set of the BSSR1 dataset obtain the estimated EER (such drawbacks as been experimented using genetic programming [41] instead of genetic algorithms).

Conclusion and Perspectives
The contribution of this paper is twofold: a fast approximated EER computing method (associated to its confidence interval), and two score fusion functions having to be parametrized thanks to genetic algorithms. Using these two contributions together allows to speed up the computation time of the genetic algorithm because its fitness function consists on computing the EER (thus, allows to use a bigger population).
The fast EER computing method has been validated on five different biometric systems and compared to the classical way. Experimental results showed the benefit of the proposed method, in terms of precision of the EER value and timing computation.
The score fusion functions have been validated on three significant multibiometrics databases (two reals and one chimerical) having a different number of scores. The fusion functions parametrized by genetic algorithm always outperform simple state of the art simple functions (sum, min, mul), and, the two new fusion functions have given better are equal results than the simple weighted sum. Using the proposed fast EER computing method also considerably speed up the timing computation of the genetic algorithms. These better results imply Figure 10: ROC Curve Of The Generated Multibiometrics Fusion Functions on the Validation Set of the Private dataset that the multibiometrics system has a better security (fewer impostors can be accepted) and is more pleasant to be used (fewer genuine users can be rejected).
One limitation of the proposed method is related to the shape of the ROC curve and the atended precision wanted. In some cases, the method is unable to get the EER at the wanted precision, and, is not able to return the result (we did not encounter this case in these experiments).
Our next research will focus on the use of different evolutionary algorithms in order to generate other kind of complex functions allowing to get better results.