Mohamed Lamine Tazir
email: tazir.med@gmail.com

Tawsif Gokhool
email: tawsif.gokhool@u-picardie.fr

Paul Checchin
email: paul.checchin@uca.fr

Laurent Malaterre
email: laurent.malaterre@uca.fr

Laurent Trassoudaine
email: laurent.trassoudaine@uca.fr

Cluster ICP: Towards Sparse to Dense Registration

Keywords: registration, dense to sparse, selection, clustering, matching

come   L'archive ouverte pluridisciplinaire

Introduction

In a generic representation of the environment, point clouds can be viewed as a collection of 3D point entities bearing a color or intensity information depending on the acquisition sensor (LiDAR, RGBD or time of flight cameras). To obtain a more meaningful information about the semantic structure of the environment, it is more instructive to rather consider a collective set of points representing the same surface. From there, several surface indices can be extracted such as their normals, curvatures and region bounds for example. Moreover, due to limitation of the field of view of 3D sensors, coupled with the complex geometry of the scanned surrounding, registration methods are required to be more robust in order to deal with data taken from large viewpoints as well as different sensor resolution.

Extraction of surface normals is a double-edged tool, which guarantees good results if accurately exploited, but can also lead to divergence of the alignment process if badly used. Since normal features are based essentially on the estimation of neighbouring points attributed to the same surface in general, they are however subjected to sensor noise, resolution and scanning patterns. Consequently, this reverberates on the alignment process due to weak inter surface correspondences between the source and the target point cloud. Furthermore, with the advancement of 3D sensors in the market, the problematic of sparse to dense registration has emerged out [START_REF] Agamennoni | Point Clouds Registration with Probabilistic Data Association[END_REF]. In this trend, software packages such as PCL [START_REF] Rusu | 3D is here: point cloud library[END_REF] has made the identification and the treatment of the above mentioned problem more accessible. Eventually, an elegant solution provided by a successful sparse to dense point cloud alignment, results in interesting robotics applications such as the case of a monocular camera localization in a 3D model [START_REF] Wolcott | Visual localization within LIDAR maps for automated urban driving[END_REF][START_REF] Caselitz | Monocular Camera Localization in 3D LiDAR Maps[END_REF] or augmenting the environment with more consistent data as in [START_REF] Maddern | Real-time probabilistic fusion of sparse 3d lidar and dense stereo[END_REF].

In order to support the claims stated above, an illustration of sparse to dense registration is given in Fig. 1. A dense point cloud is obtained from a 3D LiDAR Leica P20 3scanner, whilst the sparser one is extracted from an HDL-32E Velodyne 4 . Figures 1 (g), (h), (i), (j) are samples of various places in a scene. The 3D points of the source and target clouds are represented in blue and green respectively, whilst their normals are in white and red. Because of the large discrepancies in density between the two point clouds, registration methods based on the classical point-to-point ICP metrics fail to provide a good pose estimate. The difficulty lies in the fact that there are no direct correspondences between the source and the target point clouds. Moreover, these figures also depict the dissimilarity between normals pertaining to the same surface, which theoretically should have the same orientations. This change is due mainly to the presence of wide amount of noise, pattern scanning, distortion and varying resolutions. This is the major problem of the methods that use geometric features according to [START_REF] Serafin | Fast and Robust 3D Feature Extraction from Sparse Point Clouds[END_REF][START_REF] Holz | Registration with the Point Cloud Library PCL[END_REF].

In this paper, a novel registration method is introduced, exploiting normals differently. It does not seek for each point its nearest neighbor sharing the same normal, nor introducing normals in the error function or in the minimization process. To remediate for the disturbances in the registration framework, hence inaccuracies in the final result, a voxelization is performed on both clouds to maintain the topological details of the scene. Then for each voxel, a normal-based classification of points is done. Thereafter, only one point of each local surface is maintained for the association step. This process results into points which offers better compatibility in terms of surface representation between the source and the target point clouds. This approach is not about points sampling, but rather an improved selection of points is achieved for the matching stage. Thereon, the process evolves in a classical ICP like framework for pose estimation.

The contributions of this paper are threefold. The first main contribution is to perform an efficient voxelization method operating on top of a clustering technique for electing one representative point for each local surface. This process aims to establish two sets of points which are most likely to be matched, hence providing good correspondences. It consists of three stages: voxelization, clustering and matching. As a result, it:

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j)
1. reduces the amount of data to be processed during the matching phase, 2. improves the matching robustness by allowing only the association of compatible points, 3. avoids wrong associations that decrease alignment accuracy, 4. does not use unnecessary points that do not provide further information more than the ones used, 5. improves convergence and accuracy simultaneously.

The second contribution is that the proposed method is totally independent of the density (number of points, scanning resolution) of the two clouds, scanning patterns (nature of sensors). It takes as input point clouds of different resolution, gathered by different sensors, or with the same sensor. It is also based only on the geometric characteristics of the points, which makes it independent of weather and illumination conditions. Thirdly, normals are computed once before starting the process and are used only to distinguish the different local surfaces. They are not used in the alignment process.

The rest of this paper is organized as follows: in Section 2, an overview of the state of the art of registration methods is given. Section 3 details the proposed method. This is followed by experiments and a comparison with the state-of-the art methods in order to evaluate the proposed approach. Finally,conclusions and suggestions of future works are discussed in Section 5.

Related Work

Registration algorithms assemble two representations of an environment in a single reference frame. The problem of registration has been dealt with extensively in several studies over the last 25 years. This started with geometric approaches leading to the appearance of the Iterative Closest Point (ICP) algorithm [START_REF] Besl | A Method for Registration of 3-D Shapes[END_REF][START_REF] Chen | Object modeling by registration of multiple range images[END_REF]. ICP is used to calculate the optimal transformation fitting two point clouds by a two-step process: matching of points and minimizing a metric describing the misalignment [START_REF] Marani | A modified iterative closest point algorithm for 3D point cloud registration[END_REF]. These two-steps iterate to minimize the matching error and thus improve alignment. In the literature, two main groups of registration methods are identified:

feature-based methods (approaches based on features extraction); dense methods (approaches exploit all the points in the cloud).

Feature-based Approaches

Feature-based methods are generally used in outdoor environments [START_REF] Maddern | Real-time probabilistic fusion of sparse 3d lidar and dense stereo[END_REF]. They are based on the use of features, which may be points that are easily identified by their apparent character (position, local information contents, mathematical definition, etc.) with respect to the other points. A good feature requires stability and distinctiveness [START_REF] Serafin | Fast and Robust 3D Feature Extraction from Sparse Point Clouds[END_REF]. In other words, detected features should be consistent in all the frames. They should be robust to noise and invariant to rotation, perspective distortion and changes of scale [28,9] [11].

Features extraction from point cloud representation is well documented literature. One can find the 3D Scale Invariant Feature Transform (3DSIFT), which is an extension of the 2D version proposed by Lowe in 1999 [START_REF] Lowe | Object recognition from local scale-invariant features[END_REF]. The 3D version was adapted by the PCL [START_REF] Rusu | 3D is here: point cloud library[END_REF] community using the curvature of points instead of the intensity of pixels [START_REF] Hänsch | Comparison of 3D interest point detectors and descriptors for point cloud fusion[END_REF]. The method uses a pyramidal approach to reach the scale invariance characteristic of features. To achieve invariance against rotation, it assigns orientations to keypoints. This adds to an incomplete list of features such as FPFH [START_REF] Rusu | Fast Point Feature Histograms (FPFH) for 3D registration[END_REF], VFH (Viewpoint Feature Histogram) [START_REF] Rusu | Fast 3d recognition and pose using the viewpoint feature histogram[END_REF], CVFH (clustered viewpoint feature histogram) [START_REF] Aldoma | CAD-model recognition and 6DOF pose estimation using 3D cues[END_REF] to name a few.

However, feature extraction techniques are often cumbersome to determine and pose a problem to real-time applications [START_REF] Costa | Robust 3/6 DoF self-localization system with selective map update for mobile robot platforms[END_REF], making them unsuitable for applications that require efficiency. Furthermore, the necessity of very dense clouds are required in order to obtain good features, which compromises with the use of sparse clouds [START_REF] Agamennoni | Point Clouds Registration with Probabilistic Data Association[END_REF][START_REF] Serafin | Fast and Robust 3D Feature Extraction from Sparse Point Clouds[END_REF][START_REF] Yang | Automatic registration of large-scale urban scene point clouds based on semantic feature points[END_REF][START_REF] Velas | Collar Line Segments for fast odometry estimation from Velodyne point clouds[END_REF]. More importantly, these methods are environment specific [START_REF] Cadena | Past, Present, and Future of Simultaneous Localization And Mapping: Towards the Robust-Perception Age[END_REF], which may result in the rejection of good data [START_REF] Nieto | Scan-SLAM: Combining EKF-SLAM and scan correlation[END_REF].

Dense Approaches

Dense approaches make use of all the points from both clouds, and require an initial guess (transformation) between the two clouds, which makes them sensitive to wrong initialization [START_REF] Yang | Go-icp: Solving 3d registration efficiently and globally optimally[END_REF][START_REF] Serafin | NICP: Dense normal based point cloud registration[END_REF][START_REF] Costa | Robust 3/6 DoF self-localization system with selective map update for mobile robot platforms[END_REF]. Despite the use of all the points, these methods are generally faster than feature-based approaches [START_REF] Serafin | NICP: Dense normal based point cloud registration[END_REF].

Dense techniques are however well adapted to a classical ICP framework. As pointed out by Pomerleau [START_REF] Pomerleau | A Review of Point Cloud Registration Algorithms for Mobile Robotics[END_REF], its easy implementation and simplicity, are both its strength and its weakness. This led to the emergence of many variants of the original solution, adapted in many ways, throughout the years. At the very outset, Chen [START_REF] Chen | Object modeling by registration of multiple range images[END_REF], improved the standard ICP by using point-to-plane metric instead of the Euclidean distance error. This approach takes advantage of surface normal information to reject wrong pairing. However, this approach fails when dealing with clouds of different densities, since normals computation are affected by the change in resolution, presence of noise and distortion [START_REF] Das | Mapping, planning, and sample detection strategies for autonomous exploration[END_REF][START_REF] Holz | Registration with the Point Cloud Library PCL[END_REF].

The Normal Distributions Transform (3D-NDT) [START_REF] Magnusson | Scan registration for autonomous mining vehicles using 3D-NDT[END_REF] discretizes the 3D points with their normals in cells, where each one is modeled by a matrix, representing the probability of occupation of its points (linear, planar and spherical). Then, a non-linear optimization is performed to calculate the transformation between the two clouds. However, according to [START_REF] Das | Mapping, planning, and sample detection strategies for autonomous exploration[END_REF], the NDT is not suitable for systems with low computing power capability.

An efficient approach for dense 3D data registration was presented in [START_REF] Segal | Generalized-ICP[END_REF]. This probabilistic version of ICP called Generalized ICP (GICP) is based on a Maximum Likelihood Estimation (MLE) probabilistic model. It exploits local planar patches in both point clouds, which leads to plane-to-plane concept. Since this algorithm is pointto-plane variant of ICP, it has similar drawbacks, especially those related to normals computation. For instance in [START_REF] Holz | Registration with the Point Cloud Library PCL[END_REF], it is shown that the non-uniform point densities cause inaccurate estimates, which degrade the performance of the algorithm. Moreover, Agamennoni [START_REF] Agamennoni | Point Clouds Registration with Probabilistic Data Association[END_REF] affirmed that the GICP does not work well in outdoor unstructured environment.

Serafin [START_REF] Serafin | NICP: Dense normal based point cloud registration[END_REF] extended the GICP algorithm by using the normals in the error function and in the selection of correspondences, which according to the authors, increases the robustness of the registration.

Our approach, called CICP for Cluster Iterative Closest Point, uses an (NDT and NICP)-like representation, however, it is different from the NDT in the way it uses the points of each voxel to determine local surfaces and get one representative point from each local surface to the matching process. In contrast, NDT computes a Gaussian distribution in points of each voxel using the vicinity of each point. Whereas, NICP uses an image projection of the voxel grid representation to compute statistics, and considers each point with the local features of the surrounding surface. These features, namely normal and curvature, are calculated for each point from its neighboring points, with a computational complexity of O(K × N), where K is the number of the neighboring points used to compute each normal and N the total number of points. Additionally, these features are used later in the process of point matching between the two clouds, as opposed to our method, that does not use normals in the matching process. Because of the difference in density, pattern scanning, and presence of noise, will lead to noisy normals and, hence, inaccurate results.

Proposed Method

In this paper, a novel registration method exploiting normals is introduced. We adopt the Rusinkiewicz [START_REF] Rusinkiewicz | Efficient variants of the ICP algorithm[END_REF] decomposition and propose a new selection strategy, which aims to improve the pairing process. Figure 2 illustrates the pipeline of the proposed method.

Surface Normal Segmentation

CICP starts with the estimation of normals of the source and target point clouds using Principal Component Analysis (PCA) [START_REF] Jolliffe | Principal Component Analysis for Special Types of Data[END_REF] follows:

C = 1 k Σ k i=1 (p i -p) T (p i -p), (1)
where, C is the covariance matrix of the nearest neighbors (NNs), k is the number of considered nearest neighbours p i , and p is its corresponding centroid along the tangent plane.

Thereafter, the target cloud is subdivided into small voxels. Points belonging to each voxel are subjected to a classification process based on their normals, giving rise to different groups of points, according to the geometric variation of each voxel. Each group of points represents a local surface since they share the same normal vector. A single point is chosen from a local surface extract to be used for the matching process. The closest point to the centroid of each local surface is elected a winner.

Similarly, the source point cloud is first transformed into the reference frame of the target cloud using a pose estimate before undergoing a similar process; voxelization, normals-based classification, designation of point's representatives. At the end of these steps, the method results into two improved sets of points from the corresponding clouds. Each set contains the most probable points to match with the points of the second set (this is more particularly in the overlapping area of the two clouds, as it reflects the same geometry seen from two different viewpoints).

CICP Matching Pipeline

The main contribution of this paper is the proposal of a new selection strategy. As mentioned above, instead of matching point-to-point as the classical ICP variants, points pass through an election process, which gives rise to one representative point for each small region. These representatives appear as the most likely points to be matched between each other. These good matches ultimately result in an accurate motion between the two clouds (shown in the results section). This election process is based on 3D position of points and their normals. It consists of three sub-tasks: (1) voxelization, (2) clustering and (3) representative election. The first task performs a spatial grouping which attempts to preserve the topological information based on the 3D position of the points.

A set of 3D cubic regions (voxels) is generated where all points within the voxel have very close spatial positions. The second task bundles all points of each voxel based on their normals. Once this grouping is done, we perform the last task, which selects one point for each cluster (local surface) in the voxel for the matching process. Algorithm 1 depicts the workflow of CICP.

Voxelization It is applied in order to maintain the topological details of the scrutinized surface. As normal computation depends on the number of neighbouring points and as the resolution of points of the two clouds is different, voxelization with the same voxel size aims to generate equivalent local regions in the two clouds. A common criterion of comparison now becomes feasible. Therefore, the voxel size parameter is of paramount importance for our technique and it should be chosen carefully in order to keep the fundamental characteristics of both point clouds; be it dense or sparse with topological details. A voxel grid with cell size d is generated, where the following set of rules are verified:

Definition 1 (Sparse Cloud). A sparse cloud is a cloud C = (V, P) in which:

|P| = O(|V |).
Definition 2 (Dense Cloud). A dense cloud is a cloud C = (V, P) in which:

|P| = O(k * |V |), with k > 2.
whereby, V : set of voxels, P: set of points, O: proportionality operator. Definitions 1 and 2 are proposed to frame the notions of sparsity and density of point clouds. The voxel size is set according to the number of points in the sparse cloud, so that each voxel contains at least one point. This choice ensures a significant difference in density between the two clouds. A dense cloud, in our case, contains at least twice as many points as the sparse cloud. Otherwise, they are considered as equivalent.

At the beginning, the procedure applies a bounding box to the entire sparse cloud by finding the minimum and maximum positions of points along the three axes X, Y and Z. The number of voxels for this bounding box is determined by the number of points and the voxel size is deduced. The same procedure is applied to the dense cloud.

1. Voxel assignment: each voxel is identified by a unique linear index. If i, j, k represent the voxel indices in the X, Y , Z dimensions, respectively, numDivX, numDivY are the number of voxels along X and Y axes, the formula to encode the linear index is [START_REF] Tazir | Color-based 3D Point Cloud Reduction[END_REF]:

idx = i + j × numDivX + k × numDivX × numDivY (2)
According to (2), we assign an index idx to each point. This relationship allows direct access to the desired voxel, thereby avoiding a linear search as in [START_REF] Wiemann | Data Handling in Large-Scale Surface Reconstruction[END_REF].

2. Voxel suppression: as the shape of the point cloud is arbitrary, the step of delimiting points by a bounding box creates many empty voxels which are later pruned out. Eventually, voxelization helps to filter noise from voxels where there is insufficient occupational evidence. An illustration of the described approach is given in Fig. 3. Clustering The process of electing one point from each local surface makes them good candidates for point correspondence searching, thereby rejecting wrong matches impacting alignment accuracy. At first, all the "voxelized" points are taken and a classification method is applied to identify points belonging to the same surface. In our work, k-means clustering [START_REF] Arthur | K-means++: The Advantages of Careful Seeding[END_REF] is used the classification technique based on the normal of each point. The appropriate number of clusters (local surfaces) within each voxel is determined using the Elbow method [START_REF] Tazir | Color-based 3D Point Cloud Reduction[END_REF][START_REF] Tibshirani | Estimating the number of clusters in a data set via the gap statistic[END_REF]. An illustration of the described approach is given in Figure 3.

Grouping the point clouds using their normal aims at: improving the robustness of the matching step by only allowing the pairing of compatible points, reducing the amount of data to be processed during the matching stage.

Matching The clustering process generates a reduced, but different number of points in both clouds. These two sets of points are used for matching. To boost up the matching process, an off shelf PCL [START_REF] Rusu | 3D is here: point cloud library[END_REF] implementation of the k-d trees is used, whereby matching is achieved using L 2 norm. Outliers between the dense and sparse sets are handled using a suitable threshold.

Optimization Framework

In the case of a point-to-point metric, the error function to be minimized is given by:

E (x) = N ∑ i=1 T ( x) p i -q i 2 (3)
The localization problem of a sparse to a dense point cloud (or vice-versa) resolves to estimating the relative transformation T ( x) between point clouds {p, q} : ∀{p i , q i } ∈ R 3 . The principle of rigid body motion is applied where the transformation of a point tethered to a coordinate frame represent the whole compact body motion. For any point pair lying on the body, metric properties such as distances and orientation are preserved. This kind of body motion, discussed subsequently forms part of the special euclidean group SE(3).

Inter-frame incremental displacement is further defined as an element of the Lie groups applied on the smooth differential manifold of SE(3) [START_REF] Blanco | A tutorial on se(3) transformation parameterizations and on-manifold optimization[END_REF], also known as the group of direct affine isometries. Motion is parametrized as a twist (a velocity screw motion around an axis in space), denoted as

x = {[ω, υ]|υ ∈ R 3 , ω ∈ so(3)} ∈ se(3): ω = ω x ω y ω z , υ = υ x υ y υ z , with so(3) = { ω ∈ R 3×3 | ω = -ω ,
where ω and υ are the angular and linear velocities respectively. The reconstruction of a group action T ∈ SE(3) from the twist consists of applying the exponential map using Rodriguez formula [START_REF] Ma | An Invitation to 3-D Vision[END_REF].

Equation ( 3) is solved iteratively in a Gauss Newton fashion, where at each iteration, a new error E and a new Jacobian matrix J(0) are computed in order to obtain the update x by:

x = -J(0) T J(0) -1 J(0) T e(x) (4)
and the rigid transformation is updated as follows:

T ←-T T (x) (5)
Minimization is stopped when the error: e 2 < α occurs, or when the calculated increment becomes too small: x 2 < ε, where α and ε are predefined stop criteria.

Results

Our CICP approach is implemented in C++ without code optimization and our algorithm is thoroughly evaluated by conducting multiple experiments. The computational efficiency of the algorithm is beyond the scope of this paper. We rather focus on the methodology. The proposed method does not require any knowledge about the external orientation of the sensors at the time of acquisition, their position and orientation are estimated by the algorithm, the only requirement is that the two clouds share a tolerable overlap.

The experimental is set up as shown in Fig. 4. The centre of the two sensors; Velodyne HDL32 and that of the Leica P20 are perfectly superimposed with the help of the STANLEY Cubix cross line laser. The velodyne is then physically displaced and rotated by known translations and rotations from the graduated set up in order to perturb the 6 degrees of freedom (dof) transformation. Data acquisition is then performed under different scenarios in order to exert our CICP algorithm. Table 1 below summarizes the various experiments performed in a controlled environment. For each experiment, CICP is initialized at Identity, i.e. x = [0, 0, 0, 0, 0, 0].

Dense-Sparse Registration with CICP

Two clouds are acquired with different sensors; the denser cloud produced by a 3D LiDAR Leica P20 laser scanner and the sparser cloud with an HDL-32E Velodyne LiDAR sensor. A Leica P20 generates very detailed and dense point clouds. Depending on the resolution chosen during the scanning process, these clouds can exceed 100 millions of points for a single scan. All computations are performed on a laptop with the following specifications; Intel Core i74800MQ processor, 2.7 GHz, and 32 GB of RAM. For reasons of computational resources, we perform a sampling process using [START_REF] Tazir | Color-based 3D Point Cloud Reduction[END_REF] in order to reduce the number of points to the order of few millions without losing use-ful information. Figures 5(a), 5(d) illustrate the output of the sampling process with 986 344 and 2 732 783 points for the office and PAVIN5 scenes, respectively.On the other hand, the HDL-32E Velodyne sensor generates sparse point clouds that do not exceed 70 000 points. This represents a ratio of 14 times between the two clouds from the first environment and a ratio of 40 times, for clouds of the second environment.

Fig. 5 shows the registration process of such point clouds using the CICP method. On the left, the green cloud is from a LiDAR Leica P20 and the blue cloud is from an HDL32-E Velodyne. The corresponding results are shown on the right. In order to verify the convergence of the optimization, a comparison between the two clouds at the start and the end of the registration process is recorded together with the convergence profile obtained from the evolution of the RMSE error as a function of the number of iterations to convergence (see Fig. 6). An exit loop condition is imposed on the translation (10 -3 ) and rotation rates (10 -4 ). We run the algorithm for several indoor and outdoor scenes, with different viewpoints, as depicted in Table 1 and shown in Fig. 6. Each experiment is performed more than 20 times. For instance, for the experiment 1 (Expt 1), the displacement between the two clouds is [150, 170, 35, 5, 0, 0], where the first three values correspond to the translation in millimeters and the last three to the rotation in degrees. As for the fourth experiment, the displacement is [65, 45, 200, 0, 7, 5], which took 34 iterations for the algorithm to converge. From this table, three elements are identified that influence the registration results; inter-frame displacement, difference in density between the two clouds, as well as the nature of Table 1. CICP registration applied to mainly two compiled data sets; OFFICE and PAVIN. The resolutions of corresponding dense and sparse cloud are given along with the initial physical measured transformation from our set up given by the first row of each experiment, whilst the second row depicts the results output by our algorithm. Convergence is evaluated from the RMSE and the number of iterations required for full registration. environment (indoor or outdoor). Overall, we would like to highlight the fact that displacement between two viewing angles are quite consequent keeping in mind that dense techniques generally require an inter frame displacement since the cost function is linearized around x = 0.

Expt Envir # dense # sparse t x t y t z θ x θ y θ z RMSE # iter. -onment cloud cloud (mm) (mm) (mm) ( • ) ( • ) ( • ) (m) 1
Continuing with our discussion on the influence of initial displacement, and let us take the case of experiments Expt 3 and Expt 6, which represent a pure translation and a pure rotation, respectively. These two experiments as a sample of several experiments that we carry out, show that generally, the pure rotation requires more energy to reach the convergence with respect to the case of pure translation. Regarding the influence of density, a quick look shows that the denser the clouds becomes, the more the RMSE decreases, leading to better registration. Finally, we observe that CICP performance depends on the scene. In fact, the impact of scene affects the performance of registration as observed by the difference of the number of iterations required to reach the convergence domain between PAVIN's and office data sets. It is clear that the indoor environment performs better registration than the outdoor scene. This is possibly caused by the richness in planar regions of the former. It should not be overlooked that the outdoor environment contains a large amount of noise and outliers. This can be seen on Expt 8 and Expt 9, in which the initial displacement is the same in both experiments. However, the alignment for the office dataset requires 34 iterations to converge instead of 56 for the PAVIN data set.

For the sake of illustration, we take four experiments arbitrarily (Expt 3, Expt 4, Expt 7, Expt 8), and showe their state before and after registration with their convergence profile in Fig. 6. A closer look to the RMSE curves in the third row of this figure reveals that the residues which are far away are successfully minimized. However, the convergence begins very quickly and then stabilizes for a while before it reaches its minimum. This is mainly due to the fact that there is not a perfect point-to-point equivalence in the two pairing sets. This is quite logical and it can be explained by the large difference in density between the two clouds, the noise, and the clustering defects on the two clouds.

Comparison with Existing Methods

In order to compare our method with the existing state-of-the-art methods, we use implemented routines of PCL [START_REF] Rusu | 3D is here: point cloud library[END_REF] library for the NDT algorithm, GICP, point-to-plane ICP and simple ICP for dense methods. For the case of feature-based methods (methods based on features extraction) we also use PCL implementations of SIFT3D and FPFH to extract characteristic points from the two clouds, and use simple ICP to perform matching. The performance of each method is evaluated using three metrics: the accuracy, the relative translational error and the relative rotational error. The former describes the evolution of the root-mean-square point-to-point distance; this can be expressed mathematically as:

RMSE = 1 n Σ n i=1 E i 2 (6)
where n is the number of points and E i is the distance error between the source points and its correspondent in the target cloud in each iteration. This can be expressed as follows:

E i = Σ m i=0 p i -q i ( 7
)
where m is the total number of points in the sparse cloud. p i and q i which represent two points of the source and target cloud, respectively. The second metric is the Relative Translational Error (RTE), which measures the translation gap between the ground truth (t GT ) and the estimated (t E ) translation vectors.

RT E = t GT -t E 2 (8)
The Relative Rotational Error (RRE) is the sum of the absolute differences of the three Euler angles, calculated from the two rotation matrices R GT (ground truth rotation matrix) and R E (estimated rotation matrix). RRE is calculated by the equation 9 as:

RRE = |Roll R -1 GT R E | + |Pitch R -1 GT R E | + |Yaw R -1 GT R E | (9)
Table 2 presents the results gathered in processing two indoor and outdoor scenes with the state-of-the-art methods. Bold values show the best result. Quantitatively, the RMSE value of the indoor scene reaches 6 cm in the case of point-to-point, 6.2 cm point-to-plane ICP, 6.3 cm for the NDT and the GICP, more than 5 cm for SIFT3D and less than 3 cm for the proposed method. The maximum number of iterations for each test is fixed at 500 beyond which the algorithm is considered as not having converged if it reaches that ceiling, as is the case of the FPFH method.

Figure 7 shows comparison of convergences between different registration methods. Again CICP outperforms the state of the art. In addition to that, it is shown that CICP is robust against scene variation.

Dense-to-Dense Data

Figure 8 shows the state of the two dense clouds before and after the registration. Despite the large number of points, the final result is correctly aligned. Here is another benefit of our approach, the fact of not considering the entire set of points for matching, but only a collected set of points from each local surface, which improves the convergence speed.

Conclusion

In this paper, a novel selection technique is introduced on the forefront of an ICP framework. Throughout the experimental section, the performance of the algorithm has been demonstrated where otherwise classical state-of-the art techniques fails or performs poorly. The main highlight of this work is the way 3D surfaces are segmented in a point cloud representation using voxelization and clustering approaches. The advantages are multifold; both sparse and dense point clouds are subsampled by maintaining the geometry of the surface. Moreover, better point estimates obtained are used for later stages of matching and registration. In so doing, the convergence domain of the cost function is greatly improved leading to faster convergence of the algorithm where again, classical techniques fails. Additionally, our method of normal based segmentation not only improves on the weakness of the heterogeneous problem of sparse to dense registration but also deals with sensor noise leading to noisy normal extraction. Finally, our CICP approach improves on the precision of registration and outperforms the state of the art. This work can also be viewed as a direct solution to localization problems of a mobile robot in an a priori map constructed with different depth sensors. Further work includes an extension of our approach to vast scale outdoor mapping and localization as well as optimizing the computational time on the CPU for real time applications.

Figure 1 .

 1 Figure 1. Dense to sparse registration: (a) point cloud obtained from the Leica P20 LiDAR with 88 556 380 points; (f) point cloud obtained from an HDL-32E Velodyne with 69 984 points; (b) and (g) are their corresponding point clouds with normal vectors; (c) registration result of point cloud (a) and (f) using our proposed method; (h) normal vectors corresponding to (c); (d), (e) and (i) are exploded views of places indicated in (a); (j) is a close up view of (i).

Figure 2 .

 2 Figure 2. Overview of the CICP pipeline. Given two point clouds, CICP starts by computing the surface normals of the two clouds. It looks for points sharing the same local properties, and then elects one representative point from each local cluster. This election process is based on 3D position of points and their normals. It consists of three sub-tasks: (1) Voxelization: a set of 3D cubic regions (voxels) is generated where all voxel points have very close spatial positions.(2) Clustering: classify all points of each voxel according to their normals. (3) Matching: once this grouping step is completed, the last task consists in selecting one point from each cluster (local surface) in each voxel. Representative points serve as candidates for correspondence process. As a result, few points are used in the matching process, but which are most likely to be associated, thereby, improving on convergence and accuracy simultaneously.

  (a) Voxelized and clusterized point cloud (b) Electing one point from each cluster for the matching phase

Figure 3 .

 3 Figure 3. Voxelized/normal-based clustering for matching process.

Algorithm 1 : 3 T 4 S 5 T 8 SFigure 4 .

 134584 Figure 4. Experimental set up for data collection from Velodyne HDL32 (left) and Leica P20 (right) sensors.

Figure 5 .

 5 Figure 5. Registration results with CICP algorithm.

Figure 6 .

 6 Figure 6. CICP results applied to different data sets from various environments.

Figure 7 .

 7 Figure 7. Convergence comparison between different registration methods.

Figure 8 .

 8 Figure 8. Registration of two dense clouds of indoor scene captured by Leica P20 sensor.

Table 2 .

 2 Comparison with the state-of-the-art methods.

					Dense			Feature-based
			ICP	pt2pl	NDT	GICP	CICP	SIFT	FPFH
				ICP				3D+ICP	+ICP
	Office	RMSE (m)	0.0602	0.0620	0.0636	0.0636	0.0299	0.0516	failed
	[m, m, m, • , • , • ]	RTE (m)	0.2811	0.2482	0.2019	0.2178	0.0169	0.0191	failed
	[0, 0.5, 0.5, 20, 0, 10]	RRE ( • )	1.8476	0.5243	1.0517	0.0978	0.0144	0.7660	failed
	PAVIN	RMSE (m)	0.0836	0.0804	0.0824	0.0860	0.0347	0.0678	failed
	[m, m, m, • , • , • ]	RTE (m)	0.0365	0.0253	0.0315	0.0642	0.0092	0.021	failed
	[0, 0.5, 0.3, 0, 0, 10]	RRE ( • )	0.1919	0.1747	0.2048	0.3211	0.1198	0.1989	failed

Leica P20: http://leica-geosystems.com/

Velodyne LiDAR: http://velodynelidar.com/

PAVIN: http://www.institutpascal.uca.fr/index.php/en/the-institut-pascal/equipments


Guillaume Stempfel

Franc Liva Ralaivola

¸ois Denis

Learning from Noisy Data using Hyperplane Sampling and Sample Averages

We present a new classification algorithm capable of learning from data corrupted by a class dependent uniform classification noise. The produced classifier is a linear classifier, and the algorithm works seamlessly when using kernels. The algorithm relies on the sampling of random hyperplanes that help the building of new training examples of which the correct classes are known; a linear classifier (e.g. an SVM) is learned from these examples and output by the algorithm. The produced examples are sample averages computed from the data at hand with respect to areas of the space defined by the random hyperplanes and the target hyperplane. A statistical analysis of the properties of these sample averages is provided as well as results from numerical simulations conducted on synthetic datasets. These simulations show that the linear and kernelized versions of our algorithm are effective for learning from both noise-free and noisy data.

Introduction

Learning from noisy data is a problem of interest both from the practical and theoretical points of view. In this paper, we focus on a particular noise setting in a binary classification framework where the noise process uniformly flips the label of an example to the opposite label with a probability that depends on each class. An instance of this classification setting might be that of automatic spam filtering where a user might erroneously label regular mails as spam and conversely; it is obvious in this example that the probability of mislabelling an email depends on its true nature (spam or non spam).

From the theoretical and algorithmic point of views, there are very few simple learning strategies that are guaranteed to output a reliable classifier from data corrupted by the classification noise process depicted above. It must be noted that despite soft-margin Support Vector Machines seem to be a viable strategy to learn from noisy data and that generalization bounds for SVM expressed in terms of margin and the values of slack variables do exist, there is no result, to our knowledge, about the characteristics of the solution to the SVM quadratic program when noisy data are involved.

Here, we propose a strategy to learn a large margin classifier from noisy data. The algorithm proposed relies on the sampling of random hyperplanes that help the building of new training examples of which the correct classes are known; a linear classifier (e.g. a perceptron) is learned from these examples and output by the algorithm. The produced examples are sample averages computed from the data at hand with respect to areas of the space defined by the random hyperplanes and the target hyperplane.

The paper is organized as follows. Section 2 formalizes the problem and introduces notation. Section 3 presents the different parts of the algorithms, i.e. the computation of the sample averages and their use as inputs to an SVM classifier. Numerical simulations are presented in Section 4: they show the behavior of our algorithm on linearly separable distributions and nonlinearly separable distributions in the noise-free and noisy contexts.

2 Notation and Definition of the Problem X denotes the input space, assumed to be an Hilbert space, equipped with an inner product denoted by •. We restrict our study to the binary classification problem and the target space Y is Y {-1, +1}. Throughout the analysis, we additionally make the assumption of the existence of zero bias separating hyperplanes (i.e hyperplanes defined as w • x = 0 that pass through the origin of the space). These assumptions make our analysis seemlessly applicable when using kernels.

In order to simplify the definition and the writing of the proofs, we will consider normalized labeled examples, that is, for a pair (x, y), we will consider ( x x , y). Note that the transformation does not change the difficulty of the problem.

Definition 1 (γ-separable distributions). For γ > 0, D γ is the set of distributions on X such that for any D ∈ D γ , there exists a unit vector w * ∈ X such that P (x,y)∼D [y(w * • x) < γ] = 0. This means that given a vector w * , we consider a deterministic labelling y(x) of x according to the sign of w * • x, i.e., y(x) = sign (w * • x).

For a distribution D defined over the labelled space X × Y, D X denotes the distribution marginalized over Y.

Definition 2 (Class-conditional classification noise process y η ). Let η + , η -∈ [0, 1) such that η + + η -< 1 and define η = [η + η -]. Let the random process y η : Y → Y maps y as follows:

y η (y) = 8 > > < > > : +1 with prob. 1 -η + if y = +1 -1 with prob. η + if y = +1 +1 with prob. η -if y = -1 -1 with prob. 1 -η -if y = -1
Let γ > 0 and η + , η -∈ [0, 1). For a distribution D ∈ D γ , D η is the distribution over X × Y from which a labelled example is drawn according to the following process: (a) draw (x, y) according to D and (b) return (x, y η (y)). The noisy version S η of a random set S = {(x 1 , y 1 ), . . . , (x m , y m )} is defined likewise, i.e. S η = {(x 1 , y η (y 1 )), . . . , (x m , y η (y m )}, which, from here on, will be shorthanded as S = {(x 1 , y η 1 ), . . . , (x m , y η m )}. With the previous definitions at hand, the problem that we tackle is the learning from samples drawn from D η , for a noise-free distribution D ∈ D γ . We exhibit a learning algorithm that, given a finite sample from D η (altered by class-conditional classification noise), outputs a linear classifier w that is targetted to a low error rate on the noise free distribution D. If learning algorithms exist to learn from distributions that undergo a uniform classification noise process (with the same noise rate on both classes), as, for instance the perceptron-like approaches described in [START_REF] Bylander | Learning Linear Threshold Functions in the Presence of Classification Noise[END_REF][START_REF] Blum | A Polynomial-Time Algorithm for Learning Noisy Linear Threshold Functions[END_REF][START_REF] Cohen | Learning Noisy Perceptrons by a Perceptron in Polynomial Time[END_REF], or coresetbased strategies [START_REF] Har-Peled | Maximum margin coresets for active and noise tolerant learning[END_REF], they have not been actually specialized to the handling of class dependent classification noise. According to [START_REF] Ralaivola | CN=CPCN[END_REF], these algorithms could be extended in a straightforward way to this setting, but we do think that it is of major importance to have a dedicated learning algorithm for this problem. In addition, we truly think that the strategy presented in this paper is amenable to a full theoretical analysis showing it is a Probably Approximately Correct strategy. This is the topic of future researches.

3 The Sample Average Machine

High Level Description

The new algorithm that we propose, called Sample Average Machine (SAM), implements a very simple yet effective two-step learning strategy. Given a noisy training sample S η = {(x 1 , y η 1 ), . . . , (x m , y η m )} drawn according to D η (the noisy version of D) the algorithm, depicted in Algorithm 1, works as follows.

1. SAM creates a new sample S µ of examples of which the correct labels, according to the target hyperplane w * , are known with high probability. To produce these examples, SAM uniformly picks random vectors on the unit hypersphere defined in the subspace spanned by the vectors of the training set (this space is therefore of dimension at most m). For each random hyperplane w, at least one new point µ is added to S µ : µ is a sample estimate of the mean vector of D X restricted to one of the four subspaces delimited by w and w * (see Fig. 1). Provided m is large enough, the correct label of µ is known with high probability. 2. SAM learns a support vector classifier from S µ . The hyperplane returned by SAM is the one that will be used to classify new data.

The remaining of this section is largely devoted to the computation of the new examples the SVM is built upon. We give their definitions and provide their statistical properties. The proof heavily relies on concentration results, which makes it possible to have dimension independent results.

Generation of the Training Sample S µ

In the following, we consider a linearly separable distribution (with margin γ > 0) D ∈ D γ , of target hyperplane w * . The noise vector associated to the noise process is

η = [η + η -]. Let S = {(x 1 , y 1 ), . . . , (x m , y m )} a sample of m labelled examples drawn from D and S η = {(x 1 , y η 1 ), . . . , (x m , y η m )} is the noisy version of S.
The first step to generate the new training points is the sampling of a vector w uniformly picked on the unit hypersphere defined in the space spanned by the x 1 , . . . , x m . Given w and w * , X can be decomposed into four disjoint subspaces (as illustrated on Figure 1):

A(w) = {x ∈ X |w • x ≥ 0, y(x) = +1} B(w) = {x ∈ X |w • x < 0, y(x) = +1} E(w) = {x ∈ X |w • x < 0, y(x) = -1} F (w) = {x ∈ X |w • x ≥ 0, y(x) = -1} . (1)
Algorithm 1 Sample Average Machine (see text for details) Likewise, the samples S and S η are divided into four subsets each:

Input: m0, n, η = [η + η -], S η = {(x1,
y A(w) = y B(w) = +1, y E(w) = y F (w) = -1 for Z ∈ {A(w), B(w), E(w), F (w)} do compute μη (Z) according to equation (3) compute Pη (Z) according to equation (4) if Pη (Z) ≥ m 0 m then S µ = S µ ∪ {( μη (Z), yZ )} end if end for /*
Â(η) (w) = n xi ∈ S (η) |w • xi ≥ 0, y (η) i = +1 o B(η) (w) = n xi ∈ S (η) |w • xi < 0, y (η) i = +1 o Ê(η) (w) = n xi ∈ S (η) |w • xi < 0, y (η) i = -1 o F (η) (w) = n xi ∈ S (η) |w • xi ≥ 0, y (η) i = -1 o (2)
where the superscript (η) denotes an optional η in the definitions.

A few observations can be readily made: (a) the mean vectors of D restricted to each of the subspaces A(w), B(w), E(w), F (w) are of class +1, +1, -1, -1, respectively; (b) those mean vectors could be estimated by sample averages computed on the subsets Â(w), B(w), Ê(w), F (w) if these sets were known; (c) in the framework of classification from noisy dataset however, the only set that is accessible is S η and its corresponding subsets and these sample averages cannot be performed directly.

It nevertheless turns out that it is possible to derive unbiased estimates of the mean vectors of each subspace from S η . For a given w, it is therefore possible to compute (at most four) vectors of which the correct classes are known with high probability, provided the size of S η is large enough.

Sample Estimates of the Mean Vectors. The four vectors of interest that we would like to approximate (respectively of class +1, +1, -1, and -1) are defined for Z in Let β = 1 1-η + -η -. Consider the following vectors, which can be computed from S η :

{A(w), B(w), E(w), F (w)} by µ(Z) = E x∼D [x1l Z (x)] where 1l Z (z) = 1 if z ∈ Z
μη (A(w)) = 1 m β " (1 -η -) P x i ∈ Âη (w) xi -η -P x i ∈ F η (w) xi " μη (B(w)) = 1 m β " (1 -η -) P x i ∈ Bη (w) xi -η -P x i ∈ Êη (w) xi " μη (E(w)) = 1 m β " (1 -η + ) P x i ∈ Êη (w) xi -η + P x i ∈ Bη (w) xi " μη (F (w)) = 1 m β " (1 -η + ) P x i ∈ F η (w) xi -η + P x i ∈ Âη (w) xi " . (3)
Then, we have the following lemma.

Lemma 1. ∀D ∈ D 0 , ∀w ∈ X , ∀η + , η -∈ [0, 1) such that η + + η -< 1, ∀m ∈ N, ∀S η = {(x 1 , y η 1 ), . . . , (x m , y η m )} drawn from D η , for Z in {A(w), B(w),E(w), F (w)} the following holds: E S η ∼D η ( μη (Z)) = µ(Z).
Proof. Let us focus on the case Z = A(w) (the proof works similarly for the other three subspaces). In order to work out the proof of the equality, it suffices to observe that the random variable S η is driven by two random process. The first one is that of drawing random vectors x 1 , . . . , x m from D X while the second is the class conditional random noise process (see Definition 2), which depends upon the first one through the labels of the x i 's. We therefore have:

Focusing on the innermost expectation, droping , for sake of space, the fact that the x 1 , . . . , x m are fixed, and using the shorthand y η for y η (y 1 ) . . . y η (y m ), we have:

Eyη [ μη (A(w))] = Eyη 2 4 1 m β 0 @ (1 -η -) X x i ∈ Âη (w) xi -η -X x i ∈ F η (w) xi 1 A 3 5 = β m Eyη 2 4 (1 -η -) X w•x i ≥0 xi1 l {+1} (y η (yi)) -η -X w•x i ≥0 xi1 l {-1} (y η (yi)) 3 5 = β m Eyη 2 4 (1 -η -) 0 @ X x i ∈ Â(w) xi1 l {+1} (y η (+1)) + X x i ∈ F (w) xi1 l {+1} (y η (-1)) 1 A -η - 0 @ X x i ∈ Â(w) xi1 l {-1} (y η (+1)) + X x i ∈ F (w)
xi1 l {-1} (y η (-1))

1 A 3 5 = 1 m β 0 @ (1 -η -) 0 @ X x i ∈ Â(w) xi(1 -η + ) + X x i ∈ F (w) xi(η -) 1 A -η - 0 @ X x i ∈ Â(w) xi(η + ) + X x i ∈ F (w) xi(1 -η -) 1 A 1 A = 1 m X x i ∈ Â(w) xi = μ(A(w)).
Getting back to the full expectation, it is straightforward to check that: We therefore have proved that, given a vector w, the points computed in (3) are unbiased estimates of the mean of the subspaces delimited by w and w * .

ESη∼Dη (µ η (A(w))) =
Concentration of the Sample Estimates. We show that the distances between the sample estimates and their expected values can be bounded with high probability. To this end, we make use of McDiarmid's inequality: Theorem 1 ( [START_REF] Mcdiarmid | On the method of bounded differences[END_REF]). Let x 1 , . . . , x n be independent random variables taking values in a set X , and assume that f :

X n → R satisfies sup x 1 ,...,xn,x ′ i ∈X ˛f (x1, . . . , xn) -f (x1, . . . , xi-1, x ′ i , xi+1, . . . , xn) ˛≤ ci
for every 1 ≤ i ≤ n. Then, for every t > 0,

P {|f (x1, . . . , xn) -Ef (x1, . . . , xn)| ≥ t} ≤ 2 exp " - 2t 2 P n i=1 c 2 i « .
We have the following lemma.

Lemma 2. ∀γ > 0, ∀δ > 0, ∀D ∈ D γ , ∀w ∈ X , ∀S sample drawn from D, if |S| > 128 2 ln ( 8 δ ) γ 2
then, for Z ∈ {A(w), B(w), E(w), F (w)},

P n µ(Z) -μ(Z) > γ 32 o < δ 4 .
Proof. Again let Z = A(w). Consider S = {(x 1 , y 1 ), . . . , (x i , y i ), . . . , (x m , y m )} and S (i) = {(x 1 , y 1 ), . . . , (x ′ i , y ′ i ), . . . , (x m , y m )} two samples of m examples drawn from D that only differs on their ith variable. It can be observed that for any i in 1, . . . , m:

˛ µ(A(w)) -μ(A(w)) - ' ' 'µ(A(w)) -μ(i) (A(w)) ' ' ' ˛≤ ' ' ' μ(A(w)) -μ(i) (A(w)) ' ' ' ≤ 2 m
Hence, by Theorem 1, Lemma 1 and according to the value of m, we have

P n | µ(A(w)) -μ(A(w)) -ES∼D [ µ(A(w)) -μ(A(w)) ]| > γ 64 o < 2 exp " - γ 2 m 8192 « < δ 4 .
If σ is a Rademacher vector of size m, i.e. a random vector whose entries independently take the values +1 or -1 with probability 0.5, we observe that:

ES∼D [ µ(A(w)) -μ(A(w)) ] = ES ˆ' ' E S ′ ˆμ ′ (A(w)) ˜-μ(A(w)) ' ' ≤ E SS ′ ˆ' ' μ′ (A(w)) -μ(A(w)) ' ' ˜(triangle ineq.) = E SS ′ »' ' ' ' 1 m X x ′ i ∈ Â′ (w)) x ′ i - 1 m X x i ∈ Â(w) xi ' ' ' ' - = E SS ′ σ » 1 m ' ' ' X m i=1 σi " x ′ i 1 l Â′ (w) (x ′ i ) -xi1 l Â(w) (xi) "' ' ' - (S and S ′ are iid samples) ≤ 2ESσ » 1 m ' ' ' X m i=1 σixi1 l Â(w) (xi) ' ' ' - = 2 m ESσ » "" X m i=1 σixi1 l Â′ (w) (xi) " • " X m j=1 σjxj1 l Â(w) (xj) "" 1 2 - ≤ 2 m " ESσ h X m i,j=1 σiσjxi • xj1 l Â(w) (xi)1 l Â(w) (xj) i" 1 2 (Jensen ineq.) = 2 m " ES h X m i=1 xi • xi1 l Â(w) (xi)1 l Â(w) (xi) i" 1 2 (Eσ [σiσj] = δij) = 2 m p mD(A(w)) ≤ 2 √ m
According to the lower bound on m, we can conclude that

P µ(A(w)) -μ(A(w)) > γ 32 < δ 4 .
The following lemma is very similar to the previous one except that the random process considered is the classification noise process.

Lemma 3. Let γ > 0. ∀η + , η -∈ [0, 1) such that η + + η -< 1, ∀δ > 0, ∀D ∈ D γ , ∀w ∈ X , ∀S η sample drawn from D η , if |S η | > 128 2 ln ( 8 δ ) γ 2 (1-η + -η -)
then for Z ∈ {A(w), B(w), C(w), D(w)}, we have:

P μ(Z) -μη (Z) > γ 32 < δ 4
Proof. The lines of the proof are exactly the same as for Lemma 2 (even though the complete proof is little bit more tedious because of the way μη (Z) is defined) and we leave it to the reader.

Finally, the next proposition readily follows.

Proposition 1. Let γ > 0. ∀η + , η -∈ [0, 1) such that η + + η -< 1, ∀δ > 0, ∀n ∈ N, ∀D ∈ D γ , ∀W = {w 1 , . . . , w n }, w i ∈ X and ∀S η drawn from D η , if |S η | > 128 2 ln ( 8n δ ) γ 2 (1-η + -η -) then the following holds: ∀Z ∈ {A(w), B(w), E(w), F (w)}, P n µ(Z) -μη (Z) > γ 16 o ≤ δ 2 , ∀w ∈ W.
Proof. Suppose that Z = A(w). By Lemma 2 and Lemma 3, we know that for a given w ∈ X :

P w : μη (A(w)) -µ(A(w)) > γ 16 ≤ δ 2n .
A simple union bound argument concludes the proof:

P n ∃w ∈ W : ' ' μ(S) A η (w i ) -µ A(w i ) ' ' > γ 16 o ≤ |W |P n w : μη (A(w)) -µ(A(w)) > γ 16 o ≤ n • δ 2n = δ 2 .

Margins for the Sample Estimates

In this subsection we show that it is possible to guarantee the (normalized) margins of the newly generated points (computed according to (3)), and therefore the correct labels, with a simple criterion based on the estimated sizes of the Ẑ(w) (for Z ∈ {A(w), B(w), C(w), D(w)}).

Proceeding in a way very similar to what we did previously for the estimation of the mean vectors, we introduce the following notation: for Z ∈ {A(w), B(w),E(w), F (w)},

P (Z) = E x∼D X 1l x∈Z and P (Z(w)) = 1 m xi∈ Ẑ(w) 1 = |Z(w)|
m , which is the sample estimate of P (Z(w)).

We define the 'probability estimates' counterparts of the mean vector estimates defined in (3):

Pη (A(w)) = 1 m β " (1 -η -)| Âη (w)| -η -| F η (w)| " Pη (B(w)) = 1 m β " (1 -η -)| Bη (w)| -η -| Êη (w)| " Pη (E(w)) = 1 m β " (1 -η + )| Êη (w)| -η + | Bη (w)| " Pη (F (w)) = 1 m β " (1 -η + )| F η (w)| -η + | Âη (w)| " (4)
With these definitions at hand it is easy to get the following proposition.

Proposition 2. Let γ > 0. ∀η + , η -∈ [0, 1) such that η + + η -< 1, ∀δ > 0, ∀n ∈ N, ∀D ∈ D γ , ∀W = {w 1 , . . . , w n }, w i ∈ X and ∀S η drawn from D η , if |S η | > 128 2 ln ( 8n δ ) γ 2 (1-η + -η -)
then the following holds:

∀Z ∈ {A(w), B(w), C(w), D(w)}, P  ˛P(A(w)) -P(A(w)) ˛> 1 16 ff ≤ δ 2 , ∀w ∈ W.
Proof. The steps of the proof are exactly the same as those leading to Proposition 1.

It suffices to observe that, for instance, | Âη (w)| = xi∈ Âη 1 is used here instead of xi∈ Âη x i (see [START_REF] Cohen | Learning Noisy Perceptrons by a Perceptron in Polynomial Time[END_REF]). The proof is left to the reader.

We finally have the following result, which shows that it is always possible with high probability to keep among the four points calculated in (3), at least one that has a positive margin with respect to the target hyperplane w * .

Theorem 2. Let γ > 0. ∀η + , η -∈ [0, 1) such that η + -η -< 1, ∀δ > 0, ∀n ∈ N, ∀D ∈ D γ , ∀W = {w 1 , . . . , w n }, w i ∈ X , w i = 1, ∀S η drawn from D η , if |S η | > 128 2 ln ( 8n δ ) γ 2 (1-η + -η -)
then with probability at least 1 -δ we can generate at least n points having margin γ 17 .

Proof. Given a target hyperplane w * , any unit vector w defines 4 subspaces, and there consequently is at least, one of these subspaces Z 0 ∈ {A(w), B(w), E(w), F (w)} such that P(Z 0 ) ≥ 1 4 . By Proposition 2

P P(Z 0 ) -Pη (Z 0 ) > 1 16 ≤ δ 2 , ∀w ∈ W,
and the computed estimate Pη (Z) of the probability for this area is at least 3 16 . By Proposition 2 again, a subspace with Pη (Z 0 ) ≥ 3 16 is such that P(Z 0 ) ≥ 1 8 . Let y 0 the label of the points in Z 0 (y 0 is known for each of the subspace A(w), B(w), E(w), F (w)). As for all x ∈ Z 0 , y 0 (w * • x) ≥ γ, we have y 0 w * • µ(Z 0 ) ≥ γ 8 and then, by Proposition 1, we can lower bound the normalized margin of μη (Z 0 ) with respect to w * for a given w (recall that Z 0 depends on w):

w * • μη (Z 0 ) w * μη (Z 0 ) ≥ w * • µ(Z 0 ) -w * µ(Z 0 ) -μη (Z 0 ) w * ( µ(Z 0 ) + µ(Z) -μη (Z 0 ) ) ≥ γ 8 -γ 16 1 + γ 16 ≥ γ 17 (with probability 1 -δ n )
Henceforth, the margin of the newly computed point μη (Z 0 ) is at least γ 17 , with probability 1 -δ n , provided Pη (Z 0 ) ≥ 3 16 . A union bound argument gives that we can generate, with probability 1 -δ, a set of n hyperplanes makes it possible to generate at least n new points with a margin greater than γ 17 .

Remark 1. Theorem 2 tells us that given a noisy training sample S η and using (3), we can generate a new linearly separable training set of size n with high probability provided that S η is large enough; in addition, if a margin of γ/17 is wanted for the points of the new set, it suffices to select the generated points µ η (Z) computed on the subspaces Z such that Pη (Z) ≥ 3 16 . In Algorithm 1, the lower bound on this probability estimate is an input of the learning procedure through parameter m 0 .

Numerical Simulations

Hyperplanes in Various Dimensions For this set of experiments, we randomly generate hyperplanes in R d , with d = 2, 10, 100. For a given hyperplane w * , the data to be learned from and to be classified are uniformly distributed in the hypercube [-1; +1] d , where a separation margin of γ = 0.02 is enforced (i.e. points such that |w * • x| < 0.02 are discarded). Noise is added to the training data according to several noise vectors η.

A cross-validation procedure is performed in order to choose the soft-margin parameter C for the SVM and the value of m 0 . The values tested for C are 1, 10, 100, 1000, 10000, 100000 and 0, 5, 10, 20, 40, 80 for m 0 . The cross-validation procedure works as follows: given a noisy training set of size m, and a noisy test set of the same size, n new points are generated from the training set according to (3) and a soft-margin SVM is learned on these data; then n new points are also generated from the test set using (3) as well and the accuracy of the learned classifier is measured on this new set; the parameters (C, m 0 ) giving the lowest test error are selected for a validation procedure made on an independent validation set. The results provided in Table 1 are the error rates measured on these validation sets (10 hyperplanes are learned for each d).

It is important to note that in the cross-validation procedure, new points must be generated as it is not possible to assess the actual error rate of a classifier on a noisy sample if the noise rates for the classes at hand are not equal.

For d = 2, the size m of the noisy training and test sets is 400 and the size of the validation set is 1000. For the learning procedure, 4m points are generated from which we remove those which do not comply with the value of m 0 ; 4m points are generated as well for the test set and they are all kept. For d = 10, we use the same setting using m = 1000, and for d = 1000, m is set to 3000.

The striking feature of the results is the very good ability of SAM to handle the classification noise. Indeed, in dimension 2 and 10, the achieved error rates are very low. For d = 100 however, the performances are a little worse but still very good with regard to the amount of noise. One point that is worth noting is that because the SVM is learned on points that are averages of areas of the space, and not the points themselves, the quality of the learning might not be optimal when no noise is present.

Banana Banana is 2-dimensional nonlinearly separable classification problem, which, in addition, has a Bayes error around 10%. The data we use are those made available by Gunnar Rätsch 1 , who provides 100 training of size m = 400 and test sets of size 4900. In order to tackle this nonlinearly separable problem, we make use of a Gaussian kernel of width 1.0: we pick this value because according to [START_REF] Rätsch | Soft Margins for AdaBoost[END_REF], this seems to be the parameter allowing for the best results when used with SVM. The cross-validation procedure implemented for Banana is a mix of the cross-validation procedure described for learning hyperplanes (see above) and the one implemented in [START_REF] Rätsch | Soft Margins for AdaBoost[END_REF]: it is a 2-fold validation process based on the first two training sets corrupted by classification noise that works by learning with SAM one of the two sets and testing on a set of averaged vectors computed according to (3) on the other set. The parameters that are selected are those that give the lowest test error. The number of points created to perform the learning is 4m = 1600. Table 2 summarizes the learning results. It must be noted that only 10 train/test sets are used to compute the mean error rates and the standard deviations. Fig. 2 depicts the learned concepts together with the corresponding noisy samples. Once again, it is striking how SAM is capable of handling the noise, even for relatively high rates. Again, a very good insensitivity of the procedure can be observed.

Conclusion

We have proposed a new method, the Sample Average Machine or SAM, for the learning of data altered by class conditional classification noise. Based on a relatively intuitive idea, SAM generates from a noisy sample another labelled sample of data whose correct classes are known with high probability, and learns an SVM on this newly generated sample. The simulation results are quite satisfactory, in particular from the noise tolerance perspective. Moreover, we provide theoretical results which prove the good statistical properties of the new computed set and justify the learning on this sample.

The next step of our research on this problem consists in giving a formal proof that the learning on the generated set is equivalent to the learning on the training sample, directly drawn from the distribution, and, finally showing that SAM (or a derived algorithm) can be fitted into the PAC framework.

  and 0 otherwise; we denote by μ(Z) = 1 m xi∈ Ẑ x i the sample estimate of µ Z (again, these sample estimates cannot be computed directly).

Fig. 1 .

 1 Fig.1. The subspaces A(w), B(w), E(w), F (w) defined by a random vector w and the target vector w * ; the corresponding subsets Â(w), B(w), Ê(w), F (w) of a noise free set S (white/black disks are positive/negative examples) are automatically deduced.

  Ex 1 ,...xm∼D X ˆEy η (y 1 ),...,y η (ym) [µ η (A(w))|x1, . . . , xm] = Ex 1 ,...xm∼D X [ μ(A(w))] = µ A(w) .

Fig. 2 .

 2 Fig. 2. Banana: the first row shows the noisy data (red/blue disks are positive/negative examples) with noises [0.0 0.0], [0.2 0.4] and [0.6 0], respectively; the second row shows the concept learned by SAM.

  SVM learning */ learn an SVM classifier wSVM on S µ return wSVM end while

Table 1 .

 1 Classification error rates (in %) together with standard deviation, obtained when learning noisy hyperplanes. The first table corresponds to d = 2 and the value for the pair (C, m0) obtained by a cross-validation procedure is provided for each noise rate. The first error rates in the second table correspond to d = 10 while the second to d = 100; for these experiments, it turned out that the cross-validation procedure output C = 1, m0 = 0.

	Noise Rate [0.0 , 0.0] [0.0 , 0.2] [0.0 , 0.4] [0.0 , 0.6] [0.2 , 0.0]
	(C, m0)	(10,20) (100,40) (100,20) (100,0)	(10,10)
	Error Rate 15.1 ± 1.1 14.0 ± 1.8 14.8 ± 1.3 16.6 ± 2.8 15.3 ± 1.4 Noise Rate [0.2 , 0.2] [0.2 , 0.4] [0.4 , 0.0] [0.4 , 0.2] [0.6 , 0.0]
	(C, m0) (100,40) (100,20) (1000,10) (100,10) (10000,5)
	Error Rate 15.3 ± 1.7 19.2 ± 4.9 13.1 ± 1.5 18.2 ± 3.7 14.3 ± 1.8

Table 2 .

 2 Banana error rates and standard deviation for various noise rates. The line labelled

(C, m0) specifies the values of (C, m0) selected by a 2-fold cross-validation (see text).

ESη∼Dη ( μη (A(w))) = Ex 1 ,...xm∼D X ˆEy η (y 1 ),...,y η (ym) [ μη (A(w))|x1, . . . , xm] ˜.

http://ida.first.fraunhofer.de/projects/bench/benchmarks.htm

Noise Rate [0.0 , 0.0] [0.0 , 0.2] [0.0 , 0.4] [0.0 , 0.6] [0.2 , 0.2] [0.2 , 0.4] (C, m0) (100,20) (10,20) (1000,10) (10,5) (1000,40) (100000,40) Error Rate 0.3 ± 0.4 0.6 ± 0.8 0.5 ± 0.8 1.4 ± 1.6 2.6 ± 1.4 4.4 ± 3.6 Noise Rate [0.0 , 0.0] [0.0 , 0.2] [0.0 , 0.4] [0.0 , 0.6] [0.2 , 0.2] [0.2 , 0.4]

Error Rate 2.0 ± 0.7 3.9 ± 0.9 3.8 ± 1.1 5.6 ± 1.6 6.2 ± 1.7 6.4 ± 1.8 Error Rate 5.1 ± 0.7 6.6 ± 1.0 9.1 ± 1.0 12.6 ± 1.0 10.0 ± 1.0 15.9 ± 1.3


O Debicki

N Mareau

L Ouvry

J Foucault

S Lesecq

G Dudnik

M Correvon

INSPEX: Integrated portable multi-sensor obstacle detection device Application to navigation for visually impaired people

According to the definition of visual impairement of the World Health Organisation:

 1.6 million people suffer of blindness in EU, only 5% are fully autonomous in their daily mobility.

[WHO12]  40% of the visually impaired suffer head level accidents at least once a month, and 30% suffer a fall accident at least once a month. [START_REF] Manduchi | Mobility-Related Accidents Experienced by People with Visual Impairment[END_REF] OptekSystems Sonar Glasses

• Existing Electronic Travel Aid (ETA):

• Some wearable solutions (sunglasses, gloves…) are sometimes considered as extra prosthesis, cumbersome and stigmatizing. • Rely on a single sensor technology: ultrasound which is sensitive to multi-echo leading to wrong detections. • Perception is limited to range sensing (of the nearest target). Current experiments on environmental perception require non autonomous and non embedded devices.

GoSense

• Context:

• Obstacle detection available in high-end vehicles • Multi-sensor (LiDAR, radar, ultrasonic, vision) solutions to compensate the drawbacks of each sensor technology • Data Fusion highly demanding in terms of memory and computation to construct a representation of the vehicle surroundings

Context and objectives

• Objectives: integrate obstacle detection capabilities in a portable device to offer a virtual safety cocoon to its "user". This requires:

• Integration of several sensor technologies to ensure detection of various obstacles (material, size, shape, colour) in 3D, at different heights and ranges, in various environment conditions • Low-power light-weight small-size range sensors • Highly efficient fusion technique to be embedded on a small computing platform (typ. µCont.)

• Modular system to tackle different application domains that might benefit of such functionality to improve the "user" navigation experience and safety (e.g. drone)

Occupancy grid

User needs

A portable multi-sensor obstacle detection device Radar LiDARs Ultrasound

• Integration of a single portable device in charge of environment perception.

• various weather conditions,

• weight and size constraints,

• connection to smart environment,

• detection of obstacles at head and waist level

• Different range sensing technologies (lidar, radar, ultrasound) , each one compensating the drawbacks of another one.

User needs User req.

Radar LiDARs Ultrasound

• Integration of a single portable device in charge of environment perception.

• various weather conditions,

• weight and size constraints,

• connection to smart environment,

• detection of obstacles at head and waist level

• Different range sensing technologies (lidar, radar, ultrasound) , each one compensating the drawbacks of another one.

User

A portable multi-sensor obstacle detection device Radar LiDARs Ultrasound

• Integration of a single portable device in charge of environment perception.

• various weather conditions,

• weight and size constraints,

• connection to smart environment,

• detection of obstacles at head and waist level

• Different range sensing technologies (lidar, radar, ultrasound) , each one compensating the drawbacks of another one.

User needs User req.

System requirements

• Portable device ( 200gr, 100cm3, 500mW)

• Environment sensors for Context awareness:

• Ultrasound sensors (range 2m),

• Radar (range 10m, FoV 60°),

• short/long ranges LiDAR (3-5m/10m),

• Inertial measurement unit

• Overlap of sensors Field-of-View

• Autonomous 10 hours lifetime • The software architecture has been designed in compliance with the system requirements, i.e. ensuring modularity, software reuse.

• The same runtime is running on both subsystems in order to ease future evolution of the hardware platform with the integration of the whole software on a unique computing platform (provided it offers enough memory and computational capability).

• Traditionally sensors and post-treatment are paired. There is no access to the raw data.

• In the white cane architecture, the post-treatment is done in the measurement front-end MCU: Focus on the fusion hub

• Model of environment is based on 3D occupancy grid.

• SigmaFusion TM is based on Bayesian sensor fusion.

• SigmaFusion TM 2D has been extended to the 3D.

• SigmaFusion TM has been extended to take into account odometry and filtering over time.

• A safety cocoon of 8000 cells is scanned at 20Hz. Information is compressed and sent to the user's smartphone for rendering. • The final prototype is a smart and miniaturized device issued from the user's needs. So the specification task started at the beginning of the project and preliminary investigations have been used as guidelines for the modules developments.

A incremental approach up to the final integration • The objective of Inspex project is to transpose automotive technologies for environment perception to portable devices.

• New system requirements imply development of innovative sensor technologies and algorithms while supporting the integration. • Integrating and sharing sensor post-treatment on the same CPU increases the entropy of data manipulated.

• 2D map fusion has been extended to 3D without additional computing resources.

• A roadmap of demonstrators allows confrontation of the results from the early prototype validation with the requirements of the final demo.

DAC 2018 |

 2018 San Francisco, CA | June 24-28 H2020 INSPEX project: Augmenting capabilities of the white cane • Our motivation: Help vision impaired persons on their daily commute.

•••

  Detect obstacles up to a distance of 4m, moving at a relative speed of ~1.4 m/s worst case Adapt or develop new sensing technologies : • Development of new sensors (LR lidar, UWB radar, Ultrasound piezo) • Mechanical and electronical integration of existing sensors (SR lidars) • Develop innovative post-treatment for obstacle detection, sensor positioning • Transfer data fusion techniques from the automotive domain into the white cane:• 3D data fusion based on 3D occupancy map of the device surrounding thanks to SigmaFusion® Algorithm (Bayesian fusion) architecture with an scalable ecosystem based on two high-end MCUs (CORTEX M7)• The split in two part is appropriate for separation of concerns (i.e. sensor providers do not interfere with the SigmaFusion algorithm/IP). The software architecture has been designed in compliance with the system requirements, i.e. ensuring modularity, software reuse.FramesDetailed architecture of the proof-of-concept demonstrator (software)

•Focus

  It enables accesses to sensors raw data • A post-treatment module can also consider information from other sensors • Enables development of more efficient algorithms based on the whole sensing ecosystem.

•

  The early prototype is used in sensor experiments and firmware development. It is based on standard development kits and out of the shelf sensors. It has enabled embedded software development during the specification of the mechanical and electrical constraints. • The intermediate prototype allows connectivity with the new sensors developed in the project. It embeds the two subsystems but one micro-controller can be bypassed to try to go one step further in the software integration.

•

  A basic architecture based on two high-end MCUs allows separation of concerns and helps mitigating the risks all along the project.

This work has been partly funded by the European Union's Horizon 2020 Research and Innovation programme under grant agreement no 730953. This work was supported in part by the Swiss secretariat for education, research and innovation (SERI) under grant 16.0136 730953


David S Dean
email: david.dean@u-bordeaux1.fr

Non-equilibrium fluctuation-induced interactions

Keywords: numbers: 05.40.-a, 05.10.Gg, 05.70.Ln

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

Introduction: fluctuation-induced interactions

The most famous example of a fluctuation-induced interaction is the Casimir effect [START_REF] Mostepanenko | The Casimir Effect and its Applications[END_REF]. Here, the presence of two conducting plates of area A modifies the zero point energy of the electromagnetic field and leads to an attractive force via the celebrated formula

f = - Aπ 2 240 hc L 4 , (1)
where L is the separation between the plates. At non-zero temperature thermal fluctuations come into play, giving rise to additional temperature-dependent interactions. The full theory of the electromagnetic Casimir effect is most compactly expressed via the Lifshitz theory, which is formulated in terms of the general dielectric and conductive properties of the interacting bodies [START_REF] Lifshitz | [END_REF]. At the level of pairwise, or two body, interaction, the Lifshitz theory recovers the general theory of van der Waals or dispersion forces [START_REF] Lifshitz | [END_REF].

In 1978 Fisher and de Gennes predicted that a critical Casimir effect would exist in critical (or near-critical systems) [START_REF] Fisher | [END_REF], a prediction that has been recently experimentally confirmed [4]. This critical Casimir force arises because objects or surfaces in critical systems modify the thermal fluctuations of the order parameter. In such systems the interaction can be tuned by, for instance, chemical treatment of the surfaces so they either favor the same or differing phases, for example, in a binary critical fluid. If L is the distance between two plates immersed in a critical fluid, then the free energy is a function of L and the average force between the two plates is given by

f = - ∂ F ∂ L , (2)
where F = -T ln(Z ) and Z is the partition function for the system (T being the temperature in units where the Boltzmann constant is 1). If we consider a simple Gaussian scalar field theory with Hamiltonian

H = 1 2 dx (∇φ) 2 + m 2 φ 2 , (3)
then we can compute the partition function Z = d[φ] exp(-β H ) by integrating over fields φ with the appropriate boundary conditions. For instance, the equilibrium force for parallel plates of area A imposing the Dirichlet-Dirichlet (DD) boundary conditions or Neumann-Neumann boundary conditions is attractive and given as [5]

f = - 2AT (4π ) d-1 2 ( d-1 2 ) × k d-2 dk √ k 2 + m 2 exp -2L √ k 2 + m 2 1 -exp(-2L √ k 2 + m 2 ) , (4)
where d denotes the dimension of space and (z) is the Euler gamma function. In agreement with the more general considerations of Fisher and de Gennes, we note that the force becomes long range when the Gaussian field theory is critical;

that is, when m = 0 we have

f = - AT (d)ζ (d) (16π) d-1 2 ( d-1 2 )L d , ( 5
)
where ζ is the Riemann zeta function. A wide range of Casimir-type interactions also occur in lipid membranes, where local membrane inclusions can couple to various physical parameters of the membrane, for instance, the curvature, lipid composition and lipid order parameters such as the lipid phase (solid-liquid-gel), the thickness of lipid bilayers or even lipid orientation. Similarly, regions of lipids of varying bending and elastic energies can experience effective interactions due to membrane height fluctuations; thus, membrane fluctuations can potentially modify the phase diagram of lipid bilayers composed of lipid mixtures [6]. Also, colloids and other objects in liquid crystals can couple to the local nematic order parameter and this leads to effective interactions between such objects [7]. Many of the examples of thermal fluctuation-induced interactions we have cited occur in soft matter systems, where time scales can be very long and dynamics slow. It is therefore natural to ask how fluctuation-induced interactions evolve with time toward their equilibrium value [8][9][10]. In addition, one can analyze a number of situations where one has non-equilibrium steady states: for instance, when the system is driven by external noise that does not obey the detailed balance necessary for thermal equilibrium [9][10][11][12]. Non-equilibrium Casimir effects have also been observed and predicted in a variety of driven granular systems [13,14], as well as in systems undergoing chemical reactions [15]. The imposition of temperature [16] or other thermodynamic gradients in systems will also lead to modifications of fluctuation-induced interactions with respect to the equilibrium situation and may provide a means of tuning fluctuation-induced interactions.

Non-equilibrium situations: the problem of computing the force

Out of equilibrium we cannot derive forces as the derivatives of a free energy, and in general, the force will depend on the dynamics, non-equilibrium driving and/or gradients in the system. In order to determine the non-equilibrium force, we thus need an expression for the force in terms of the field configuration. There are basically three approaches to this problem:

• The phenomenological approach. In non-thermal systems such as granular materials and systems of reaction diffusion equations, one can define a local pressure based on local kinetic arguments. For instance, the local order parameter φ can be associated with a local density ρ(φ) and the local pressure defined via a local equation of state, for instance P = Tρ, if one uses the ideal gas approximation [15]. In granular materials, the equation of state for hard spheres has also been used [13].

• The stress tensor. The instantaneous force on a surface S in the direction i is given by

f i = S T i j dS j , (6)
where T i j is the stress tensor; this approach has been widely used in both equilibrium and non-equilibrium situations [8,11,12,17].

• The principle of virtual work. Boundary conditions are effectively imposed by assuming that the surface has a coordinate L and a field φ effectively exists on the surface and interacts with the field via an interaction term [9,10,18], for instance,

H int = dx δ(z -L)V (φ), (7)
for a surface on the plane z = L. The total energy is thus a function H [φ, L] of both the field φ and the coordinate L. The force conjugate to the parameter L is thus given by the principle of virtual work as f = -∂ H ∂ L . As an example, taking V (φ) = λφ 2 /2 and then taking the limit λ → ∞ imposes Dirichlet boundary conditions on the surface at z = L.

It is straightforward to see that the principle of virtual work and stress tensor gives the same expressions for the average force in equilibrium. However, the computation of the force seems to require clarification of the field theory at a microscopic level; for instance, one needs to know how the field φ changes when the surface is displaced [10,19]. In electrostatics, there is no ambiguity as the force acting on a volume V is given by f = V dx ρE, where E is the local electric field and ρ is the charge density [START_REF] Landau | Electrodynamics of Continuous Media[END_REF]. From this expression, we can, in fact, derive the stress tensor expression for the force. In a binary mixture, the forces on the plates will be generated by the interaction between the molecules composing the plate and the molecules of the surrounding fluid; in a molecular dynamics simulation the force on the plates can thus be unambiguously computed. However, if we coarse grain the binary fluid and replace it with an Ising model, we immediately see that the numerical calculation of the instantaneous force for a given spin configuration is much less obvious.

The results for fluctuation-induced interactions governed by stochastic dissipative dynamics

Here we consider a free field theory with a general quadratic Hamiltonian

H = 1 2 dx dx φ(x) (x, x , L)φ(x ), (8)
and we will compute the force using the principle of virtual work discussed in the previous section. Here is a self-adjoint operator, i.e. (x, x ) = (x , x), and L represents any suitable free parameter in the problem, but for concreteness it could be the position of a plate that interacts with the field. For concreteness we consider the problem where the system is prepared in a state φ = 0 at the time t = 0 and then let it relax at some non-zero temperature T close to a critical point where it is massless. We also assume that the system is of infinite extent; that is to say the field exists outside the two plates. This means that the resulting force per unit area is a disjoining pressure; it has no bulk term and thus tends to zero when the plate separation becomes large.

We assume that the field obeys general over-damped stochastic dynamics [START_REF] Chaikin | Principles of Condensed Matter Physics[END_REF] ∂φ

(x) ∂t = -dx R(x, x ) δ H δφ(x ) + η(x, t) = -dx R (x, x )φ(x ) + η(x, t), ( 9
)
where R indicates the composed operator R (x, x ) = dy R(x, y) (y, x ). To satisfy detailed balance with noise that is uncorrelated in time, we chose the spatial noise correlator to be

η(x, t)η(x , t ) = 2T δ(t -t )R(x, x ). ( 10
)
When R(x, x ) = δ(xx ) we have non-conserved model A dynamics, and when R(x, x ) = -∇ 2 δ(xx ) we have conserved model B dynamics. It can be shown [9,10] that the Laplace transform of the time-average-dependent force is given by

L f (s) = T s ∂ ∂ L ln(Z ( s )), (11)
where the operator s is given by

s = + s 2 R -1 (12)
and Z ( s ) is the partition function for a Gaussian field theory defined via the operator s . This result means that the Laplace transform of the time-dependent Casimir force can be expressed in terms of the static Casimir force for another free field theory, and if this static partition function is known the time-dependent force can be extracted by inverting the Laplace transform. The late time equilibrium result is recovered from the pole at s = 0 in equation (11). Note that in the case of model B dynamics, the operator s is non-local due to the term R -1 = -∇ -2 . Thus, intriguingly, the study of the temporal evolution of fluctuation-induced forces can lead us to consider static Casimir forces arising from non-local field theories. In the case of model B dynamics, additional boundary conditions have to be specified for the flux of the field or its chemical potential at the interfaces [START_REF] Diehl | [END_REF]. This choice of boundary conditions still needs to be specified within the formalism described here, and the cases when it can be applied are yet to be worked out-the non-equilibrium evolution of the Casimir force under model B dynamics thus remains an open problem even for Gaussian field theories.

Let us consider the case of model A dynamics, i.e. where the dynamical operator R(xx ) = δ(xx ). This is the easiest case to analyze and it is the case that has been most studied in the literature via the other approaches mentioned above. We take H to be of the form of equation (3) in the critical case m = 0 and impose DD boundary conditions.

The Laplace transform can be almost inverted analytically, the temporal derivative of the out of equilibrium force is given by

d f dt = - 2AT (8π ) d 2 t d-1 2 ∂ ∂t ∞ n=1 1 √ t exp - L 2 n 2 2t , ( 13
)
where A is the area of the plates. Using the Poisson summation formula, we find that

d f dt = - AT 2(8π) d 2 t d+2 2 + AT 2(8π) d-1 2 t d+1 2 × ∂ ∂ L ∞ n=1 exp - 2π 2 n 2 t L 2 . ( 14
)
As model A dynamics is diffusive, we see that L 2 sets a time scale. The short-term behavior of the force as defined by the regime t/L 2 1 can be obtained from equation (13) and is given by

f (t) ∼ - 2AT (8π t) d 2 exp - L 2 2t . ( 15
)
The late time, defined via t/L 2 1, decay to the equilibrium force f eq can be extracted from equation ( 14) and is given by

f (t) ∼ f eq + AT d(8π ) d 2 t d 2 . ( 16
)
In equation ( 14), one sees the appearance of the eigenvalues for DD boundary conditions in the sum on the right-hand side. The first term can be seen to be due to the bulk on the exterior of the system. The results here also agree with those predicted using the stress tensor [8]; in the case of DD boundary conditions, the two approaches can be seen to be equivalent even out of equilibrium.

We now turn to a model of interacting dipoles which has the advantage of predicting a temporarily evolving Casimir, or van der Waals, force that can be, in principle, confronted with experiment and where, due to its electrostatic origin, the force can be precisely defined.

We consider a model dielectric medium [23] made up of local polarization fields p ν (x) at the point x of the medium. Here ν corresponds to a type or species of polarization field that notably has its own polarizability per unit volume denoted by χ ν (x). Next, we take two semi-infinite regions (slabs) V + and V -defined via the sign of the coordinate z, such that z > 0 in V + and z < 0 in V -. The two regions are separated in the z-direction by a distance L. The total energy for a given configuration of the dipole fields is

H = 1 2 dx dy νν p ν (x) • A νν (x, y, L)p ν (y), (17)
where

A νν (x, y, L) = δ(x -y)I δ νν χ (x) + D(x, y, L). ( 18
)
Here I is the identity matrix in three-dimensional space and the polarization energy equation ( 17) corresponds to the classical harmonic energy needed to generate a local polarization field. The second term D is the usual dipole-dipole interaction. We assume that the dipole field dynamics is of the model A type

∂ p iν (x) ∂t = -κ ν (x) δ H δ p iν (x) + ζ iν (x, t), ( 19
)
where κ ν (x) determines the local relaxation in the region x of space for dipoles of type ν and ζ iν (x, t) is Gaussian white noise in time with spatial correlation function satisfying detailed balance. The choice of model A type dynamics is based on the observation that polarization is not a conserved quantity (i.e. it can be generated locally at a given point without the need for charge transport to that point). The parameters in the above model can then be related to the frequency-dependent dielectric function of the slabs via

(ω) = 0 + ν χ ν 1 + iωχ ν κ ν ; (20)
this form of the dielectric function further justifies the choice of model A dynamics as it has the form of a superposition of the classical Debye expression for the dielectric function.

We now consider the situation where two such dielectric slabs are brought suddenly into proximity, parallel and at distance L along their axes. The development of the van der Waals force between the two slabs is physically due to the onset of the correlation between the dipoles in the two slabs. It be can be shown [23] that the time-dependent Casimir force is given by f (t) = -∂ F(t,L) ∂ L , where F is a time-dependent free energy whose Laplace transform is given by

LF = - T s ln Z -i s 2 , (21)
where Z [ (-i s 2 )] is the partition function for the static problem with

Z -i s 2 = d[φ] exp - β 2 dx x, - is 2 [∇φ(x)] 2 .
(22) The static result, which is identical to that predicted by Lifshitz theory, is recovered from the pole at s = 0. The Laplace transform is, in general, difficult to invert; however, in the case of the parallel slab geometry, it can be shown that the force is given by

f (L , t, T ) = - T A L 3 H (t), (23)
where H (t) is a time-dependent Hamaker constant that starts at zero and evolves to its static value as t → ∞.

We thus see that the scaling with distance L and time is completely different from that of model A dynamics for a massless scalar field, even though the static values have the same 1/L 3 scaling. The behavior of H (t) depends on the choice of dielectric function. For simple Debye-type dielectric functions, the decay to the equilibrium force is exponential; however, for certain dielectric functions [24], typical of glassy polymer-based materials, the decay can be of a power-law type [23].

Conclusions and open questions

We have seen that results on the out of equilibrium dynamics of fluctuation-induced forces are restricted to the simplest models of stochastic dynamics and that the field is replete with open questions and new directions for future study.

To conclude this paper, I list a number of questions and challenges I believe would be interesting to address.

• The definition of the out of equilibrium force. There seem to be a number of unresolved issues as to how one can define the instantaneous force in out of equilibrium situations for coarse-grained models. As well as permitting the determination of out of equilibrium forces, such methods, especially numerical ones, could lead to better schemes for measuring Casimir forces in simulations, even in equilibrium. Notably, such new methods will be useful to determine the fluctuations of fluctuation-induced interactions [17]. • Force dynamics for interacting field theories. Even within the framework of simple dissipative stochastic dynamics, the only existing results on dynamics are for free or Gaussian field theories. The treatment of interacting field theories remains an open challenge. A problem that seems attackable is the problem of an N components φ 4 field theory in the limit N → ∞ whose dynamics can be solved [26]. • Beyond overdamped stochastic dynamics. Clearly, the dynamics of binary liquid mixtures will be sensitive to hydrodynamic effects in most geometries and the coupling between hydrodynamics and the dynamics of the order parameter [START_REF] Chaikin | Principles of Condensed Matter Physics[END_REF]26] will change the temporal evolution of fluctuation-induced forces. Understanding the interplay between phase ordering in the presence of plates or inclusions and hydrodynamics is thus essential. Also it would be interesting, from a purely academic point of view, to see what effect the inclusion of inertial effects has on the evolution of fluctuationinduced forces; this should give the force propagation a wave-like component. Possibly, models for fluctuationinduced interactions in granular materials [13,14] could incorporate inertial effects. • Non-equilibrium steady states. One of the most active fields in the quantum electrodynamic Casimir effect is the study of forces between objects held at different temperatures [25]. In the realm of classical fluctuationinduced interactions, less work has been done. It can be shown that in the parallel plate geometry, if the three regions (outside left, inside and outside right) are held at different temperatures then the effective interaction is dominated by a bulk-like term [16] that, is mathematically similar in origin to the bulk difference in radiation pressure between regions of space at different temperatures. In soft matter systems, the absence of vacuum makes it natural to study the effect of temperature gradients, imposed for instance by holding the plates at different temperatures. The interacting dipole model presented in the previous section presents an interesting test bed for exploring the effects of temperature gradients and may indeed be realistic enough to make experimentally verifiable predictions.

Acknowledgments

I am grateful to V Démery, S Dietrich, J-B Fournier, A Gambassi, A J Gopinathan, V A Parsegian and R Podgornik for stimulating discussions and their collaboration on these topics and to the European Science Foundation and Institut Universitaire de France for financial support.