Loic Landrieu
email: loic.landrieu@ign.fr

Martin Simonovsky
email: martin.simonovsky@enpc.fr

Large-scale Point Cloud Semantic Segmentation with Superpoint Graphs

We propose a novel deep learning-based framework to tackle the challenge of semantic segmentation of largescale point clouds of millions of points. We argue that the organization of 3D point clouds can be efficiently captured by a structure called superpoint graph (SPG), derived from a partition of the scanned scene into geometrically homogeneous elements. SPGs offer a compact yet rich representation of contextual relationships between object parts, which is then exploited by a graph convolutional network. Our framework sets a new state of the art for segmenting outdoor LiDAR scans (+11.9 and +8.8 mIoU points for both Semantic3D test sets), as well as indoor scans (+12.4 mIoU points for the S3DIS dataset).

Introduction

Semantic segmentation of large 3D point clouds presents numerous challenges, the most obvious one being the scale of the data. Another hurdle is the lack of clear structure akin to the regular grid arrangement in images. These obstacles have likely prevented Convolutional Neural Networks (CNNs) from achieving on irregular data the impressive performances attained for speech processing or images.

Previous attempts at using deep learning for large 3D data were trying to replicate successful CNN architectures used for image segmentation. For example, SnapNet [5] converts a 3D point cloud into a set of virtual 2D RGBD snapshots, the semantic segmentation of which can then be projected on the original data. SegCloud [START_REF] Tchapmi | SEGCloud: Semantic segmentation of 3D point clouds[END_REF] uses 3D convolutions on a regular voxel grid. However, we argue that such methods do not capture the inherent structure of 3D point clouds, which results in limited discrimination performance. Indeed, converting point clouds to 2D format comes with loss of information and requires to perform surface reconstruction, a problem arguably as hard as semantic segmentation. Volumetric representation of point clouds is Both authors contributed equally to this work. inefficient and tends to discard small details.

Deep learning architectures specifically designed for 3D point clouds [START_REF] Qi | PointNet: Deep learning on point sets for 3D classification and segmentation[END_REF][START_REF] Simonovsky | Dynamic edgeconditioned filters in convolutional neural networks on graphs[END_REF]42,[START_REF] Qi | PointNet++: Deep hierarchical feature learning on point sets in a metric space[END_REF]10] display good results, but are limited by the size of inputs they can handle at once.

We propose a representation of large 3D point clouds as a collection of interconnected simple shapes coined superpoints, in spirit similar to superpixel methods for image segmentation [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]. As illustrated in Figure 1, this structure can be captured by an attributed directed graph called the superpoint graph (SPG). Its nodes represent simple shapes while edges describe their adjacency relationship characterized by rich edge features.

The SPG representation has several compelling advantages. First, instead of classifying individual points or voxels, it considers entire object parts as whole, which are easier to identify. Second, it is able to describe in detail the relationship between adjacent objects, which is crucial for contextual classification: cars are generally above roads, ceilings are surrounded by walls, etc. Third, the size of the SPG is defined by the number of simple structures in a scene rather than the total number of points, which is typically several order of magnitude smaller. This allows us to model long-range interaction which would be intractable otherwise without strong assumptions on the nature of the pairwise connections. Our contributions are as follows:

• We introduce superpoint graphs, a novel point cloud representation with rich edge features encoding the contextual relationship between object parts in 3D point clouds.

• Based on this representation, we are able to apply deep learning on large-scale point clouds without major sacrifice in fine details. Our architecture consists of Point-Nets [START_REF] Qi | PointNet: Deep learning on point sets for 3D classification and segmentation[END_REF] for superpoint embedding and graph convolutions for contextual segmentation. For the latter, we introduce a novel, more efficient version of Edge-Conditioned Convolutions [START_REF] Simonovsky | Dynamic edgeconditioned filters in convolutional neural networks on graphs[END_REF] as well as a new form of input gating in Gated Recurrent Units [START_REF] Cho | Learning phrase representations using RNN encoder-decoder for statistical machine translation[END_REF].

• We set a new state of the art on two publicly available datasets: Semantic3D [START_REF] Hackel | Semantic3d. net: A new large-scale point cloud classification benchmark[END_REF] and S3DIS [START_REF] Armeni | 3D semantic parsing of largescale indoor spaces[END_REF]. In particu-

Related Work

The classic approach to large-scale point cloud segmentation is to classify each point or voxel independently using handcrafted features derived from their local neighborhood [48]. The solution is then spatially regularized using graphical models [37,24,[START_REF] Lu | Simplified Markov random fields for efficient semantic labeling of 3d point clouds[END_REF]44,[START_REF] Kim | 3D scene understanding by voxel-CRF[END_REF]2,[START_REF] Niemeyer | Contextual classification of lidar data and building object detection in urban areas[END_REF]35,49] or structured optimization [START_REF] Landrieu | A structured regularization framework for spatially smoothing semantic labelings of 3D point clouds[END_REF]. Clustering as preprocessing [16,13] or postprocessing [START_REF] Weinmann | A hybrid semantic point cloud classification-segmentation framework based on geometric features and semantic rules[END_REF] have been used by several frameworks to improve the accuracy of the classification.

Deep Learning on Point Clouds. Several different approaches going beyond naive volumetric processing of point clouds have been proposed recently, notably setbased [START_REF] Qi | PointNet: Deep learning on point sets for 3D classification and segmentation[END_REF][START_REF] Qi | PointNet++: Deep hierarchical feature learning on point sets in a metric space[END_REF], tree-based [42,[START_REF] Klokov | Escape from cells: Deep Kd-networks for the recognition of 3D point cloud models[END_REF], and graph-based [START_REF] Simonovsky | Dynamic edgeconditioned filters in convolutional neural networks on graphs[END_REF]. However, very few methods with deep learning components have been demonstrated to be able to segment large-scale point clouds. PointNet [START_REF] Qi | PointNet: Deep learning on point sets for 3D classification and segmentation[END_REF] can segment large clouds with a sliding window approach, therefore constraining contextual information within a small area only. Engelmann et al. [10] improves on this by increasing the context scope with multi-scale windows or by considering directly neighboring window positions on a voxel grid. SEGCloud [START_REF] Tchapmi | SEGCloud: Semantic segmentation of 3D point clouds[END_REF] handles large clouds by voxelizing followed by interpolation back to the original resolution and post-processing with a conditional random field (CRF). None of these approaches is able to consider fine details and long-range contextual information simultaneously. In contrast, our pipeline partitions point clouds in an adaptive way according to their geometric complexity and allows deep learning architecture to use both fine detail and interactions over long distance.

Graph Convolutions. A key step of our approach is using graph convolutions to spread contextual information. Formulations that are able to deal with graphs of variable sizes can be seen as a form of message passing over graph edges [12]. Of particular interest are models supporting continuous edge attributes [START_REF] Simonovsky | Dynamic edgeconditioned filters in convolutional neural networks on graphs[END_REF][START_REF] Monti | Geometric deep learning on graphs and manifolds using mixture model CNNs[END_REF], which we use to represent interactions. In image segmentation, convolutions on graphs built over superpixels have been used for postprocessing: Liang et al. [START_REF] Liang | Semantic object parsing with graph LSTM[END_REF][START_REF] Liang | Interpretable structure-evolving LSTM[END_REF] traverses such graphs in a sequential node order based on unary confidences to improve the final labels. We update graph nodes in parallel and exploit edge attributes for informative context modeling. Xu et al. [START_REF] Xu | Scene graph generation by iterative message passing[END_REF] convolves information over graphs of object detections to infer their contextual relationships. Our work infers relationships implicitly to improve segmentation results. Qi et al. [START_REF] Qi | 3D graph neural networks for RGBD semantic segmentation[END_REF] also relies on graph convolutions on 3D point clouds. However, we process large point clouds instead of small RGBD images with nodes embedded in 3D instead of 2D in a novel, rich-attributed graph. Finally, we note that graph convolutions also bear functional similarity to deep learning formulations of CRFs [START_REF] Zheng | Conditional random fields as recurrent neural networks[END_REF], which we discuss more in Section 3.4.

Method

The main obstacle that our framework tries to overcome is the size of LiDAR scans. Indeed, they can reach hundreds of millions of points, making direct deep learning approaches intractable. ing to a geometrically simple primitive, which we assume to be semantically homogeneous. Such primitives can be reliably represented by downsampling small point clouds to at most hundreds of points. This small size allows us to utilize recent point cloud embedding methods such as PointNet [START_REF] Qi | PointNet: Deep learning on point sets for 3D classification and segmentation[END_REF].

3 Contextual segmentation: The graph of superpoints is by orders of magnitude smaller than any graph built on the original point cloud. Deep learning algorithms based on graph convolutions can then be used to classify its nodes using rich edge features facilitating longrange interactions.

The SPG representation allows us to perform end-to-end learning of the trainable two last steps. We will describe each step of our pipeline in the following subsections.

Geometric Partition with a Global Energy

In this subsection, we describe our method for partitioning the input point cloud into parts of simple shape. Our objective is not to retrieve individual objects such as cars or chairs, but rather to break down the objects into simple parts, as seen in Figure 3. However, the clusters being geometrically simple, one can expect them to be semantically homogeneous as well, i.e. not to cover objects of different classes. Note that this step of the pipeline is purely unsupervised and makes no use of class labels beyond validation.

We follow the global energy model described by [13] for its computational efficiency. Another advantage is that the segmentation is adaptive to the local geometric complexity. In other words, the segments obtained can be large simple shapes such as roads or walls, as well as much smaller components such as parts of a car or a chair.

Let us consider the input point cloud C as a set of n 3D points. Each point i ∈ C is defined by its 3D position p i , and, if available, other observations o i such as color or intensity. For each point, we compute a set of d g geometric features f i ∈ R dg characterizing the shape of its local neighborhood. In this paper, we use three dimensionality values proposed by [START_REF] Demantk | Dimensionality based scale selection in 3D lidar point clouds[END_REF]: linearity, planarity and scattering, as well as the verticality feature introduced by [13]. We also compute the elevation of each point, defined as the z coordinate of p i normalized over the whole input cloud.

The global energy proposed by [13] is defined with respect to the 10-nearest neighbor adjacency graph G nn = (C, E nn ) of the point cloud (note that this is not the SPG). The geometrically homogeneous partition is defined as the constant connected components of the solution of the following optimization problem: is chosen to be linearly decreasing with respect to the edge length. The factor µ is the regularization strength and determines the coarseness of the resulting partition.

arg min g∈R dg i∈C g i -f i 2 + µ (i,j)∈Enn w i,j [g i -g j = 0] , (1)
The problem defined in Equation 1 is known as generalized minimal partition problem, and can be seen as a continuous-space version of the Potts energy model, or an 0 variant of the graph total variation. The minimized functional being nonconvex and noncontinuous implies that the problem cannot realistically be solved exactly for large point clouds. However, the 0 -cut pursuit algorithm introduced by [26] is able to quickly find an approximate solution with a few graph-cut iterations. In contrast to other optimization methods such as α-expansion [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF], the 0 -cut pursuit algorithm does not require selecting the size of the partition in advance. The constant connected components

S = {S 1 , • • • , S k } of

Superpoint Graph Construction

In this subsection, we describe how we compute the SPG as well as its key features. The SPG is a structured representation of the point cloud, defined as an oriented attributed graph G = (S, E, F ) whose nodes are the set of superpoints S and edges E (referred to as superedges) represent the adjacency between superpoints. The superedges are annotated by a set of d f features: F ∈ R E×d f characterizing the adjacency relationship between superpoints.

We define G vor = (C, E vor ) as the symmetric Voronoi adjacency graph of the complete input point cloud as defined by [20]. Two superpoints S and T are adjacent if there is at least one edge in E vor with one end in S and one end in T :

E = (S, T ) ∈ S 2 | ∃ (i, j) ∈ E vor ∩ (S × T ) .
(2)

Important spatial features associated with a superedge (S, T ) are obtained from the set of offsets δ(S, T ) for edges in E vor linking both superpoints:

δ (S, T ) = {(p i -p j ) | (i, j) ∈ E vor ∩ (S × T )} . (3)
Superedge features can also be derived by comparing the shape and size of the adjacent superpoints. To this end, we compute |S| as the number of points comprised in a superpoint S, as well as shape features length (S) = λ 1 , surface (S) = λ 1 λ 2 , volume (S) = λ 1 λ 2 λ 3 derived from the eigenvalues λ 1 , λ 2 , λ 3 of the covariance of the positions of the points comprised in each superpoint, sorted by decreasing value. In Table 1, we describe a list of the different superedge features used in this paper. Note that the break of symmetry in the edge features makes the SPG a directed graph.

Superpoint Embedding

The goal of this stage is to compute a descriptor for every superpoint S i by embedding it into a vector z i of fixed-size dimensionality d z . Note that each superpoint is embedded in isolation; contextual information required for its reliable classification is provided only in the following stage by the means of graph convolutions.

Several deep learning-based methods have been proposed for this purpose recently. We choose PointNet [START_REF] Qi | PointNet: Deep learning on point sets for 3D classification and segmentation[END_REF] for its remarkable simplicity, efficiency, and robustness. In PointNet, input points are first aligned by a Spatial Transformer Network [START_REF] Jaderberg | Spatial transformer networks[END_REF], independently processed by multilayer perceptrons (MLPs), and finally max-pooled to summarize the shape.

In our case, input shapes are geometrically simple objects, which can be reliably represented by a small amount of points and embedded by a rather compact PointNet. This is important to limit the memory needed when evaluating many superpoints on current GPUs. In particular, we subsample superpoints on-the-fly down to n p = 128 points to maintain efficient computation in batches and facilitate data augmentation. Superpoints of less than n p points are sampled with replacement, which in principle does not affect the evaluation of PointNet due to its max-pooling. However, we observed that including very small superpoints of less than n minp = 40 points in training harms the overall performance. Thus, embedding of such superpoints is set to zero so that their classification relies solely on contextual information.

In order for PointNet to learn spatial distribution of different shapes, each superpoint is rescaled to unit sphere before embedding. Points are represented by their normalized position p i , observations o i , and geometric features f i (since these are already available precomputed from the partitioning step). Furthermore, the original metric diameter of the superpoint is concatenated as an additional feature after PointNet max-pooling in order to stay covariant with shape sizes.

Contextual Segmentation

The final stage of the pipeline is to classify each superpoint S i based on its embedding z i and its local surroundings within the SPG. Graph convolutions are naturally suited to this task. In this section, we explain the propagation model of our system.

Our approach builds on the ideas from Gated Graph Neural Networks [30] and Edge-Conditioned Convolutions (ECC) [START_REF] Simonovsky | Dynamic edgeconditioned filters in convolutional neural networks on graphs[END_REF]. The general idea is that superpoints refine their embedding according to pieces of information passed along superedges. Concretely, each superpoint S i maintains its state hidden in a Gated Recurrent Unit (GRU) [START_REF] Cho | Learning phrase representations using RNN encoder-decoder for statistical machine translation[END_REF]. The hidden state is initialized with embedding z i and is then processed over several iterations (time steps) t = 1 . . . T . At each iteration t, a GRU takes its hidden state h i to superpoint i is computed as a weighted sum of hidden states h (t) j of neighboring superpoints j. The actual weighting for a superedge (j, i) depends on its attributes F ji,• , listed in Table 1. In particular, it is computed from the attributes by a multi-layer perceptron Θ, so-called Filter Generating Network. Formally:

h (t+1) i = (1 -u (t) i ) q (t) i + u (t) i h (t) i q (t) i = tanh(x (t) 1,i + r (t) i h (t) 1,i ) u (t) i = σ(x (t) 2,i + h (t) 2,i ), r (t) i = σ(x (t) 3,i + h (t) 3,i ) (4)
(h

(t) 1,i , h (t)
2,i , h

(t) 3,i ) T = ρ(W h h (t) i + b h ) (x (t) 1,i , x (t) 2,i , x (t) 3,i ) T = ρ(W x x (t) i + b x ) (5) x (t) i = σ(W g h (t) i + b g ) m (t) i (6) m (t) i = mean j|(j,i)∈E Θ(F ji,• ; W e ) h (t) j (7) h (1) i = z i , y i = W o (h (1)
i , . . . , h

(T +1) i ) T , (8)
where is element-wise multiplication, σ(•) sigmoid function, and W • and b • are trainable parameters shared among all GRUs. Equation 4 lists the standard GRU rules [START_REF] Cho | Learning phrase representations using RNN encoder-decoder for statistical machine translation[END_REF] with its update gate u

i and reset gate r i , with being a small constant. Finally, the model includes three interesting extensions in Equations 6-8, which we detail below.

Input Gating. We argue that GRU should possess the ability to down-weight (parts of) an input vector based on its hidden state. For example, GRU might learn to ignore its context if its class state is highly certain or to direct its attention to only specific feature channels. Equation 6 achieves this by gating message m (t) i by the hidden state before using it as input x (t) i .

Edge-Conditioned Convolution. ECC plays a crucial role in our model as it can dynamically generate filtering weights for any value of continuous attributes F ji,• by processing them with a multi-layer perceptron Θ. In the original formulation [START_REF] Simonovsky | Dynamic edgeconditioned filters in convolutional neural networks on graphs[END_REF] (ECC-MV), Θ regresses a weight matrix to perform matrix-vector multiplication Θ(F ji,• ; W e )h (t) j for each edge. In this work, we propose a lightweight variant with lower memory requirements and fewer parameters, which is beneficial for datasets with few but large point clouds. Specifically, we regress only an edge-specific weight vector and perform element-wise multiplication as in Equation 7 (ECC-VV). Channel mixing, albeit in an edge-unspecific fashion, is postponed to Equation 5. Finally, let us remark that Θ is shared over time iterations and that self-loops as proposed in [START_REF] Simonovsky | Dynamic edgeconditioned filters in convolutional neural networks on graphs[END_REF] are not necessary due to the existence of hidden states in GRUs.

State Concatenation. Inspired by DenseNet [START_REF] Huang | Densely connected convolutional networks[END_REF], we concatenate hidden states over all time steps and linearly transform them to produce segmentation logits y i in Equation 8. This allows to exploit the dynamics of hidden states due to increasing receptive field for the final classification.

Relation to CRFs. In image segmentation, postprocessing of convolutional outputs using Conditional Random Fields (CRFs) is widely popular. Several inference algorithms can be formulated as (recurrent) network layers amendable to end-to-end learning [START_REF] Zheng | Conditional random fields as recurrent neural networks[END_REF][START_REF] Schwing | Fully connected deep structured networks[END_REF], possibly with general pairwise potentials [33, 7, 28]. While our method of information propagation shares both these characteristics, our GRUs operate on d z -dimensional intermediate feature space, which is richer and less constrained than low-dimensional vectors representing beliefs over classes, as also discussed in [START_REF] Gadde | Superpixel convolutional networks using bilateral inceptions[END_REF]. Such enhanced access to information is motivated by the desire to learn a powerful representation of context, which goes beyond belief compatibilities, as well as the desire to be able to discriminate our often relatively weak unaries (superpixel embeddings). We empirically evaluate these claims in Section 4.3.

Further Details

Adjacency Graphs. In this paper, we use two different adjacency graphs between points of the input clouds: G nn in Section 3.1 and G vor in Section 3.2. Indeed, different definitions of adjacency have different advantages. Voronoi adjacency is more suited to capture long-range relationships between superpoints, which is beneficial for the SPG. Nearest neighbors adjacency tends not to connect objects separated by a small gap. This is desirable for the global energy but tends to produce a SPG with many small connected components, decreasing embedding quality. Fixed radius adjacency should be avoided in general as it handles the variable density of LiDAR scans poorly.

Training. While the geometric partitioning step is unsupervised, superpoint embedding and contextual segmentation are trained jointly in a supervised way with cross entropy loss. Superpoints are assumed to be semantically homogeneous and, consequently, assigned a hard ground truth label corresponding to the majority label among their contained points. We also considered using soft labels corresponding to normalized histograms of point labels and training with Kullback-Leibler [25] divergence loss. It performed slightly worse in our initial experiments, though.

Naive training on large SPGs may approach memory limits of current GPUs. We circumvent this issue by randomly subsampling the sets of superpoints at each iteration and training on induced subgraphs, i.e. graphs composed of subsets of nodes and the original edges connecting them. Specifically, graph neighborhoods of order 3 are sampled to select at most 512 superpoints per SPG with more than n minp points, as smaller superpoints are not embedded. Note that as the induced graph is a union of small neighborhoods, relationships over many hops may still be formed and learned. This strategy also doubles as data augmentation and a strong regularization, together with randomized sampling of point clouds described in Section 3.3. Additional data augmentation is performed by randomly rotating superpoints around the vertical axis and jittering point features by Gaussian noise N (0, 0.01) truncated to [-0.05, 0.05].

Testing. In modern deep learning frameworks, testing can be made very memory-efficient by discarding layer activations as soon as the follow-up layers have been computed. In practice, we were able to label full SPGs at once. To compensate for randomness due to subsampling of point clouds in PointNets, we average logits obtained over 10 runs with different seeds.

Experiments

We evaluate our pipeline on the two currently largest point cloud segmentation benchmarks, Semantic3D [START_REF] Hackel | Semantic3d. net: A new large-scale point cloud classification benchmark[END_REF] and Stanford Large-Scale 3D Indoor Spaces (S3DIS) [START_REF] Armeni | 3D semantic parsing of largescale indoor spaces[END_REF], on both of which we set the new state of the art. Furthermore, we perform an ablation study of our pipeline in Section 4.3.

Even though the two data sets are quite different in nature (large outdoor scenes for Semantic3D, smaller indoor scanning for S3DIS), we use nearly the same model for both. The deep model is rather compact and 6 GB of GPU memory is enough for both testing and training. We refer to Appendix A for precise details on hyperparameter selection, architecture configuration, and training procedure.

Performance is evaluated using three metrics: per-class intersection over union (IoU), per-class accuracy (Acc), and overall accuracy (OA), defined as the proportion of correctly classified points. We stress that the metrics are computed on the original point clouds, not on superpoints.

Semantic3D

Semantic3D [START_REF] Hackel | Semantic3d. net: A new large-scale point cloud classification benchmark[END_REF] is the largest available LiDAR dataset with over 3 billion points from a variety of urban and rural scenes. Each point has RGB and intensity values (the latter of which we do not use). The dataset consists of 15 training scans and 15 test scans with withheld labels. We also evaluate on the reduced set of 4 subsampled scans, as common in past work.

In Table 2, we provide the results of our algorithm compared to other state of the art recent algorithms and in Figure 3, we provide qualitative results of our framework. Our framework improves significantly on the state of the art of semantic segmentation for this data set, i.e. by nearly 12 mIoU points on the reduced set and by nearly 9 mIoU points on the full set. In particular, we observe a steep gain on the "artefact" class. This can be explained by the ability of the partitioning algorithm to detect artifacts due to their singular shape, while they are hard to capture using snapshots, as suggested by [5]. Furthermore, these small object are often merged with the road when performing spatial regularization.

Stanford Large-Scale 3D Indoor Spaces

The S3DIS dataset [START_REF] Armeni | 3D semantic parsing of largescale indoor spaces[END_REF] consists of 3D RGB point clouds of six floors from three different buildings split into individual rooms. We evaluate our framework following two dominant strategies found in previous works. As advocated by [START_REF] Qi | PointNet: Deep learning on point sets for 3D classification and segmentation[END_REF]10], we perform 6-fold cross validation with microaveraging, i.e. computing metrics once over the merged predictions of all test folds. Following [START_REF] Tchapmi | SEGCloud: Semantic segmentation of 3D point clouds[END_REF], we also report the performance on the fifth fold only (Area 5), corresponding to a building not present in the other folds. Since some classes in this data set cannot be partitioned purely using geometric features (such as boards or paintings on walls), we concatenate the color information o to the geometric features f for the partitioning step.

The quantitative results are displayed in Table 3, with qualitative results in Figure 3 and in Appendix D. S3DIS is a difficult dataset with hard to retrieve classes such as white boards on white walls and columns within walls. From the quantitative results we can see that our framework performs better than other methods on average. Notably, doors are able to be correctly classified at a higher rate than other approaches, as long as they are open, as illustrated in Figure 3. Indeed, doors are geometrically similar to walls, but their position with respect to the door frame allows our network to retrieve them correctly. On the other hand, the partition merges white boards with walls, depriving the network from the opportunity to even learn to classify them: the IoU of boards for theoretical perfect classification of superpoints (as in Section 4.3) is only 51.3.

Computation Time. In Table 4, we report computation time over the different steps of our pipeline for the inference on Area 5 measured on a 4 GHz CPU and GTX 1080 Ti GPU. While the bulk of time is spent on the CPU for partitioning and SPG computation, we show that voxelization as pre-processing, detailed in Appendix A, leads to a significant speed-up as well as improved accuracy.

Ablation Studies

To better understand the influence of various design choices made in our framework, we compare it to several baselines and perform an ablation study. Due to the lack of public ground truth for test sets of Semantic3D, we evaluate on S3DIS with 6-fold cross validation and show comparison of different models to our Best model in Table 5.

Performance Limits. The contribution of contextual segmentation can be bounded both from below and above.

The lower bound (Unary) is estimated by training PointNet with d z = 13 but otherwise the same architecture, denoted as PointNet13, to directly predict class logits, without SPG and GRUs. The upper bound (Perfect) corresponds to assigning each superpoint its ground truth label, and thus sets the limit of performance due to the geometric partition. We can see that contextual segmentation is able to win roughly 22 mIoU points over unaries, confirming its importance. Nevertheless, the learned model still has room of up to 26 mIoU points for improvement, while about 12 mIoU points are forfeited to the semantic inhomogeneity of superpoints. CRFs. We compare the effect of our GRU+ECCbased network to CRF-based regularization. As a baseline (iCRF), we post-process Unary outputs by CRF inference over SPG connectivity with scalar transition matrix, as described by [13]. Next (CRF -ECC), we adapt CRF-RNN framework of Zheng et al. [START_REF] Zheng | Conditional random fields as recurrent neural networks[END_REF] to general graphs with edgeconditioned convolutions (see Appendix B for details) and train it with PointNet13 end-to-end. Finally (GRU13), we modify Best to use PointNet13. We observe that iCRF barely improves accuracy (+1 mIoU), which is to be expected, since the partitioning step already encourages spatial regularity. CRF -ECC does better (+15 mIoU) due to end-to-end learning and use of edge attributes, though it is still below GRU13 (+18 mIoU), which performs more complex operations and does not enforce normalization of the embedding. Nevertheless, the 32 channels used in Best instead of the 13 used in GRU13 provide even more freedom for feature representation (+22 mIoU).

Ablation. We explore the advantages of several design choices by individually removing them from Best in order to compare the framework's performance with and without them. In NoInputGate we remove input gating in GRU; in NoConcat we only consider the last hidden state in GRU for output as in ECC -VV we use the proposed lightweight formulation of ECC. We can see that each of the first two choices accounts for about 5 mIoU points. Next, without edge features our method falls back even below iCRF to the level of Unary, which validates their design and overall motivation for SPG. ECC -VV decreases the performance on the S3DIS dataset by 3 mIoU points, whereas it has improved the performance on Semantic3D by 2 mIoU. Finally, we invite the reader to Appendix C for further ablations.

y i = W o h (T +1) i instead of

Conclusion

We presented a deep learning framework for performing semantic segmentation of large point clouds based on a partition into simple shapes. We showed that SPGs allow us to use effective deep learning tools, which would not be able to handle the data volume otherwise. Our method significantly improves on the state of the art on two publicly available datasets. Our experimental analysis suggested that future improvements can be made in both partitioning and learning deep contextual classifiers.

The source code in PyTorch as well as the trained models are available at https://github.com/loicland/ superpoint_graph. everywhere. Input points have d p =11 dimensional features for Semantic3D (position p i , color o i , geometric features f i ), with 3 additional ones for S3DIS (room-normalized spatial coordinates, as in past work [START_REF] Qi | PointNet: Deep learning on point sets for 3D classification and segmentation[END_REF]).

Segmentation Network. We use embedding dimensionality d z = 32 and T = 10 iterations. ECC-VV is used for Semantic3D (there are only 15 point clouds even though the amount of points is large), while ECC-MV is used for S3DIS (large number of point clouds). Filter-generating network Θ is a MLP with 4 layers (widths 32, 128, 64, and 32 or 32 2 for ECC-VV or ECC-MV) with ReLUs. Batch Normalization is used only after the third parametric layer. No bias is used in the last layer. Superedges have d f = 13 dimensional features, normalized by mean subtraction and scaling to unit variance based on the whole training set.

Training. We train using Adam [START_REF] Kingma | A method for stochastic optimization[END_REF] with initial learning rate 0.01 and batch size 2, i.e. effectively up to 1024 superpoints per batch. For Semantic3D, we train for 500 epochs with stepwise learning rate decay of 0.7 at epochs 350, 400, and 450. For S3DIS, we train for 250 epochs with steps at 200 and 230. We clip gradients within [-1, 1].

B. CRF-ECC

In this section, we describe our adaptation of CRF-RNN mean field inference by Zheng et al. [START_REF] Zheng | Conditional random fields as recurrent neural networks[END_REF] for postprocessing PointNet embeddings in SPG, denoted as unary potentials U i here.

The original work proposed a dense CRF with pairwise potentials Ψ defined to be a mixture of m Gaussian kernels as Ψ ij = µ m w m K m (F ij ), where µ is label compatibility matrix, w are parameters, and K are fixed Gaussian kernels applied on edge features.

We replace this definition of the pairwise term with a Filter generating network Θ [START_REF] Simonovsky | Dynamic edgeconditioned filters in convolutional neural networks on graphs[END_REF] parameterized with weights W e , which generalizes the message passing and compatibility transform steps of Zheng et al. . Furthermore, we use superedge connectivity E instead of assuming a complete graph. The pseudo-code is listed in Algorithm 1. Its output are marginal probability distributions Q. In practice we run the inference for T = 10 iterations.

Algorithm 1 CRF-ECC Q i ← softmax(U i ) while not converged do Qi ← j|(j,i)∈E Θ(F ji,• ; W e )Q j Qi ← U i -Qi Q i ← softmax( Qi ) end while

C. Extended Ablation Studies

In this section, we present additional set of experiments to validate our design choices and present their results in Table 6.

a) Spatial Transformer Network. While STN makes superpoint embedding orientation invariant, the relationship with surrounding objects are still captured by superedges, which are orientation variant. In practice, STN helps by 4 mIoU points.

b) Geometric Features. Geometric features f i are computed in the geometric partition step and can therefore be used in the following learning step for free. While Point-Nets could be expected to learn similar features from the data, this is hampered by superpoint subsampling, and therefore their explicit use helps (+4 mIoU).

c) Sampling Superpoints. The main effect of subsampling SPG is regularization by data augmentation. Too small a sample size leads to disregarding contextual information (-4 mIoU) while too large a size leads to overfitting (-2 mIoU). Lower memory requirements at training is an extra benefit. There is no subsampling at test time.

d) Long-range Context. We observe that limiting the range of context information in SPG harms the performance. Specifically, capping distances in G vor to 1 m (as used in PointNet [START_REF] Qi | PointNet: Deep learning on point sets for 3D classification and segmentation[END_REF]) or 5 m (as used in SegCloud [START_REF] Tchapmi | SEGCloud: Semantic segmentation of 3D point clouds[END_REF]) worsens the performance of our method (even more on our Semantic 3D validation set).

e) Input Gate. We evaluate the effect of input gating (IG) for GRUs as well as LSTM units. While a LSTM unit achieves higher score than a GRU (-3 mIoU), the proposed IG reverses this situation in favor of GRU (+1 mIoU). Unlike the standard input gate of LSTM, which controls the information flow from the hidden state and input to the cell, our IG controls the input even before it is used to compute all other gates. f) Regularization Strength µ. We investigate the balance between superpoints' discriminative potential and their homogeneity controlled by parameter µ . We observe that the system is able to perform reasonably over a range of SPG sizes.

g) Superpoint Sizes. We include a breakdown of superpoint sizes for µ = 0.03 in relation to hyperparameters n minp = 40 and n p = 128, showing that 93% of points are in embedded superpoints, and 79% in superpoints that are subsampled.

Superedge Features. Finally, in Table 7 we evaluate empirical importance of individual superedge features by removing them from Best. Although no single feature is crucial, the most being offset deviation (+3 mIoU), we remind the reader than without any superedge features the net-Furthermore, SegCloud divides the inference into cubes without overlap, possibly causing inconsistencies across boundaries.

D. Video Illustration

We provide a video illustrating our method and qualitative results on S3DIS dataset, which can be viewed at https://youtu.be/Ijr3kGSU_tU.

1 arXiv: 1711 Figure 1 :

 117111 Figure 1: Visualization of individual steps in our pipeline. An input point cloud (a) is partitioned into geometrically simple shapes, called superpoints (b). Based on this preprocessing, a superpoints graph (SPG) is constructed by linking nearby superpoints by superedges with rich attributes (c). Finally, superpoints are transformed into compact embeddings, processed with graph convolutions to make use of contextual information, and classified into semantic labels.

Figure 2 :

 2 Figure2: Illustration of our framework on a toy scan of a table and a chair. We perform geometric partitioning on the point cloud (a), which allows us to build the superpoint graph (b). Each superpoint is embedded by a PointNet network. The embeddings are then refined in GRUs by message passing along superedges to produce the final labeling (c).

  where [•] is the Iverson bracket. The edge weight w ∈ R |E| +

  (t)i and an incoming message m (t) i as input, and computes its new hidden state h (t+1) i . The incoming message m (t)

  improve stability during training, in Equation 5 we apply Layer Normalization [4] defined as ρ(a) := (a-mean(a))/(std(a)+ ) separately to linearly transformed input x

Figure 3 :

 3 Figure 3: Example visualizations on both datasets. The colors in (b) are chosen randomly for each element of the partition.

Figure 5 :

 5 Figure 5: Histogram of points contained in superpoints of different size (in log scale) on the full S3DIS dataset. The embedding threshold n minp and subsampling threshold n p are marked in red.

  lar, we improve mean per-class intersection over union (mIoU) by 11.9 points for the Semantic3D reduced test set, by 8.8 points for the Semantic3D full test set, and by up to 12.4 points for the S3DIS dataset.

Table 1 :

 1 the solution of Equation 1 define our geometrically simple elements, and are referred as super-List of d f = 13 superedge features characterizing the adjacency between two superpoints S and T . points (i.e. set of points) in the rest of this paper.

	Feature name	Size	Description
	mean offset	3	mean m∈δ(S,T ) δ m
	offset deviation	3	std m∈δ(S,T ) δ m
	centroid offset	3	mean i∈S p i -mean j∈T p j
	length ratio	1	log length (S) /length (T )
	surface ratio	1	log surface (S) /surface (T )
	volume ratio	1	log volume (S) /volume (T )
	point count ratio	1	log |S|/|T |

Table 2 :

 2 Intersection over union metric for the different classes of the Semantic3D dataset. OA is the global accuracy, while mIoU refers to the unweighted average of IoU of each class.

	(a) RGB point cloud	(b) Geometric partitioning	(c) Prediction	(d) Ground truth

Table 3 :

 3 Results on the S3DIS dataset on fold "Area 5" (top) and micro-averaged over all 6 folds (bottom). Intersection over union is shown split per class.

	Step	Full cloud	2 cm	3 cm	4 cm
	Voxelization	0	40	24	16
	Feature computation	439	194	88	43
	Geometric partition	3428	1013	447	238
	SPG computation	3800	958	436	252
	Inference	10 × 24	10 × 11 10 × 6 10 × 5
	Total	7907	2315	1055	599
	mIoU 6-fold	54.1	60.2	62.1	57.1

Table 4 :

 4 Computation time in seconds for the inference on S3DIS Area 5 (68 rooms, 78 649 682 points) for different voxel sizes.

Table 5 :

 5 concatenation of all steps; in NoEdgeFeat we perform homogeneous regularization by setting all superedge features to scalar 1; and Ablation study and comparison to various baselines on S3DIS (6-fold cross validation).

	Model	mAcc mIoU
	Best	73.0	62.1
	Perfect	92.7	88.2
	Unary	50.8	40.0
	iCRF	51.5	40.7
	CRF -ECC	65.6	55.3
	GRU13	69.1	58.5
	NoInputGate	68.6	57.5
	NoConcat	69.3	57.7
	NoEdgeFeat	50.1	39.9
	ECC -VV	70.2	59.4

Table 6 :

 6 Ablation study of design decisions on S3DIS (6fold cross validation). Our choices in bold.

	Model	mAcc mIoU
	Best	73.0	62.1
	no mean offset	72.5	61.8
	no offset deviation	71.7	59.3
	no centroid offset	74.5	61.2
	no len/surf/vol ratios	71.2	60.7
	no point count ratio	72.7	61.7

Table 7 :

 7 Ablation study of superedge features on S3DIS (6fold cross validation).

		•10 6
		3
	number of points	1 2
		40 128	1,000	10,000
		size of superpoints

Superpoint embedding: Each node of the SPG corresponds to a small part of the point cloud correspond-

Appendix

A. Model Details

Voxelization. We pre-process input point clouds with voxelization subsampling by computing per-voxel mean positions and observations over a regular 3D grid (5 cm bins for Semantic3D and 3 cm bins for S3DIS dataset). The resulting semantic segmentation is interpolated back to the original point cloud in a nearest neighbor fashion. Voxelization helps decreasing the computation time and memory requirement, and improves the accuracy of the semantic segmentation by acting as a form of geometric and radiometric denoising as well (Table 4 in the main paper). The quality of further steps is practically not affected, as superpoints are usually strongly subsampled for embedding during learning and inference anyway (Section 3.3 in the main paper). Geometric Partition. We set regularization strength µ = 0.8 for Semantic3D and µ = 0.03 for S3DIS, which strikes a balance between semantic homogeneity of superpoints and the potential for their successful discrimination (S3DIS is composed of smaller semantic parts than Semantic3D). In addition to five geometric features f (linearity, planarity, scattering, verticality, elevation), we use color information o for clustering in S3DIS due to some classes being geometrically indistinguishable, such as boards or doors.


Nicolas Coppo

Mathieu Darnet

François Bretaudeau

Pierre Wawrzyniak

Julien Porte

F Rembert

Thomas Jacob

Claudine Lamotte

Characterization of intermediate depth reservoirs in complex coastal areas with CSEM method

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

Characterization of intermediate depth reservoirs in complex coastal areas

with CSEM method 1 BRGM, n.coppo@brgm.fr 1 BRGM, m.darnet@brgm.fr 1 BRGM, f.bretaudeau@brgm.fr 1 BRGM, p.wawrzyniak@brgm.fr 1 BRGM, j.porte@brgm.fr 2 METIS-UPMC, flore.rembert@upmc.fr 1 BRGM, t.jacob@brgm.fr 1 BRGM, c.lamotte@brgm.fr

SUMMARY

The characterization of water resources at intermediate depths (200-1500 m) in coastal peri-urban environments requires the implementation of specific and complementary geophysical techniques to overcome as much as possible the influence of anthropogenic noise and problems related to human infrastructures. Within the frame of the DEM'EAUX THAU project aiming at understanding and ultimately better manage the interactions between meteoritic, marine and deep thermal water in a karstic carbonate reservoir, we undertook to acquire 3D electromagnetic data around Balaruc-les-Bains (Hérault, France) to complete the hydrogeological and geological investigations. 91 controlled-source electromagnetic (CSEM) stations were collected both on-shore (56) and off-shore (35) on the Etang de Thau over ~25 km 2 in the frequency range 0.125 -8192 Hz using 2 powerful 22kVA transmitters running simultaneously at opposite side of the survey area.

This paper mainly focuses on the analysis of the CSEM data in obvious coastal 3D context. As a first step, apparent resistivity maps show a good agreement with both gravity data and expected geology of the area. Then, 1-D inversion (based on OCCAM1D, Key (2009)) was tried but failed due to the complexity of the area, mainly caused by the highly conductive Etang de Thau lake (~0.16 Ωm) and its rugged geometry. In a third step, 2.5-D inversions along profiles crossing the Etang de Thau have been tried using MARE2DEM (Key and Ovall, 2011) to deal with the shallow highly conductive layer but only partially reconstructed the resistivity distribution at depth due to the too simplistic 2.5D assumption.

Finally 3-D modelling considering the water layer and its complete bathymetry has been performed using the POLYEM3D software (Bretaudeau et al. 2018this workshop). The simulations highlight that even though the depth of the water does not exceed 15m, significant charge accumulations all along the coast are visible and depend on the depth and shape of the water layer. These effects are of the order of magnitude of the anomalies sought and can clearly not be neglected during the inversion. Current work consists in performing 3D inversion using the POLYEM3D inversion code by including all the complexity of the water layer.

Keywords: Controlled-Source Electro-Magnetics, On and Off-shore CSEM acquisition, Groundwater management, 3-D EM modelling and inversion


Laurent Bako

Guillaume Mercère

Stéphane Lecoeuche

Marco Lovera

RECURSIVE SUBSPACE IDENTIFICATION OF HAMMERSTEIN MODELS BASED ON LS-SVM

Keywords: LS-SVM, identification, subspace methods, recursive algorithms, nonlinear systems

This paper presents a recursive scheme for the identification of Hammerstein MIMO models. The Markov parameters of the system are determined first by a Least Squares Support Vector Machines (LS-SVM) regression through an over-parameterization technique. Then, a state space realization of the system is retrieved using an adapted online subspace identification method. Simulation results are provided to demonstrate the effectiveness of the algorithm in the presence of white output noise.

INTRODUCTION

System identification based on input-output data measurements is a necessary step in most applications such as control design, plant diagnosis and monitoring. A number of methods have been developed for the identification of linear systems as well in input-output form as in state space representation. Some of the proposed techniques have even been extended to slowly time varying systems [START_REF] Lovera | Recursive subspace identification of linear and non-linear Wiener state-space models[END_REF], [START_REF] Mercère | A new recursive method for subspace identification of noisy systems: EIVPM[END_REF], [START_REF] Oku | Recursive 4SID algorithms using gradient type subspace tracking[END_REF]. But nonlinear non-stationary systems identification in state space form still remains a very challenging research area. To that purpose the identification of a Hammerstein system is considered in this paper. This model consists of a nonlinear memoryless function f (•) followed by a linear dynamic system. The difficulty in separating the nonlinear part from the linear one is related to the fact that their respective terms are present only in crossproduct forms in the system equations. A typical way to proceed is to treat each of these crossproduct terms as an unknown. That results in an increased number of unknowns and the procedure is referred to as the over-parameterization method [START_REF] Bai | An optimal two-stage identification algorithm for Hammerstein-Wiener nonlinear systems[END_REF], (Goethals et al., 2005b). Generally, the static nonlinearity f (•) is expanded on a basis of known elementary functions (Radial Basis Functions, polynomials, B-splines, . . . ) and then, the coordinates associated to it are estimated. [START_REF] Dempsey | Identification of Hammerstein models with cubic spline nonlinearities[END_REF], [START_REF] Ramos | Identification of nonlinear systems using a b-splines parametric subspace approach[END_REF], [START_REF] Vörös | Modeling and identification of wiener systems with two-segment nonlinearities[END_REF]. Many proposed methods are based on input-output models and very few are based on state space models which may be sometimes more convenient for describing systems and especially multivariable systems. The features shared by most of these methods are that they often identify an input-output behavior of the system and use either iterative optimization or overparameterization or basis expansion technique to deal with the non-linearity. The present paper is inspired by (Goethals et al., 2005a) which firstly demonstrated the efficiency of componentwise LS-SVM in Hammerstein systems identification. Note that the proposed approach is very similar in the results to the well-known technique of expansion the nonlinearity on gaussian functions basis. Reliable results have been obtained for Linear Time Invariant systems. However, in actual problems there exist only a few number of such systems; almost any real life system is subject to some parameters variation during its operating time. In order to make up for this evolution, we extend in this paper the work presented in (Goethals et al., 2005a) to time-varying systems by using LS-SVM to estimate the nonlinear part of the system and ordinary least squares for recovering the linear part in state space form.

PROBLEM STATEMENT

Consider a discrete-time SIMO Hammerstein state space system represented as follows

x t+1 = Ax t + Bf (u t ) + w t y t = Cx t + Df (u t ) + v t (1)
where y t ∈ R ny and x t ∈ R nx are the output and state vectors respectively. The input u t ∈ R is a white gaussian sequence and w t ∈ R nx , v t ∈ R ny are some zero-mean noises. Moreover, the system under consideration is assumed to be stable and of known order n x . f (•) is a scalar-valued smooth nonlinear function. Given N consecutive inputoutput data measurements {u t , y t } N t=1 , the following work aims at determining the parameters (A, B, C, D) and f of the system in (1) up to a similarity transformation.

IDENTIFICATION

Let us write y t as a d-step ahead predictive model using equation (1). The equation obtained is very similar to a Finite Impulse Response (FIR) model

y t = CA d x t-d + d i=0 H i f (u t-i )+ d i=1 CA i-1 w t-i + v t et (2) ∀t > d,
where d is a positive integer, e t is the equation error and H i , defined as

H i = D if i = 0 CA i-1 B if i ≥ 1 (3)
are the so-called Markov parameters and represent the impulse responses of the system. These parameters do not depend on the basis of the state space representation in (1). Equation (2) amounts to

y t (s) = C s A d x t-d + d i=0 H i (s)f (u t-i ) + e t (s) (4) s = 1, • • • , n y
, where y t (s) is the sth output, C s the s-th row of the observation matrix C.

In (Goethals et al., 2005a), a method for off-line identification of Hammerstein models using LS-SVM regression techniques was proposed. This method produces reliable results but is computationally heavy since it requires costly singular value decomposition (SVD) steps (for computing the states sequences and for extracting the nonlinearity). Here, we would like to investigate a way of avoiding the SVD computation in this method for online application purposes. Although it is feasible to estimate a fully parameterized MIMO system by applying twice the LS-SVM regression procedure as in (Goethals et al., 2005a), only the case of SIMO system is considered in this paper for the sake of simplicity. Anyway, the system may be multi-input; the only requirement concerns the nonlinearity f (•) which has to be a scalar-valued function.

We will first estimate the terms H i f (•) and then the extended observability matrix. In equation ( 4), we choose the integer d at least larger than the order of the system and sufficiently large so that the term CA d x t-d may be considered negligible with respect to the noise term. This approximation is all the more acceptable as the system has rapid dynamics or d is relatively large. Then, equation (4) becomes

y t (s) ≈ d i=0 H i (s)f (u t-i ) + e t (s) (5)
for s = 1, . . . , n y . In the case that the system does not present rapid dynamics it may be safer to use the observer-augmenting procedure described in [START_REF] Juang | Identification of observer/Kalman filter Markov parameters : theory and experiments[END_REF] instead. This latter algorithm introduces an observer gain to set the poles of the resulting system as close to zero as possible so that the accuracy of the previous equation is noticeably improved even for small values of d.

Recall that neither the parameters H i nor the nonlinear function f (•) are known. The LS-SVM regression procedure developed by Suykens [START_REF] Suykens | Least Squares Support Vector Machines[END_REF] consists in modeling the nonlinearity as f (•) = ω ϕ(•) + δ where ϕ : R nu → R n h is a nonlinear mapping from the input space toward a possibly high dimensional feature space R n h , δ is a bias compensator and ω ∈ R n h . However, direct application of this idea results here in a hard non-convex optimization problem (Goethals et al., 2005b). Fortunately, this difficulty can be circumvented by combin-ing over-parameterization and componentwise LS-SVM [START_REF] Pelckmans | Componentwise least squares support vector machines[END_REF]. Letting

g i,s (•) = H i (s)f (•) = ω i,s ϕ(•) + δ i,s (6)
we get from (5)

y t (s) = d i=0 g i,s (u t-i ) + e t (s) = (7) = d i=0 ω T i,s ϕ(u t-i ) + δ s + e t (s) with δ s = d i=0 δ i,s , s = 1, • • • , n y , t = d + 1, • • • , N .
Since only input-output data (u t , y t ) are available, we would like to compute the nonlinear components g i,s from the sample values of their sum along i. That corresponds to the so-called component-wise LS-SVM problem. One can then notice that the solution of equation ( 7) is not unique. In fact, for any solution

{g i,s (•)} d i=0 of (7), {g i,s (•) + ε i,s (•)} d i=0 is also so- lution for all arbitrary functions {ε i,s (•)} d i=0 , sat- isfying d i=0 ε i,s (u t-i ) = 0, ∀ t = d + 1, • • • , N .
Hence, each function g i,s can only be estimated up to an additive undetermined function term. A typical way to overcome that difficulty is to introduce some additional constraint on the ω i,s ϕ(•) to be obtained by solving (7) so that the redundancy is removed. In (Goethals et al., 2005a), this constraint is chosen as a centering relation over the time window of the measurement points. An alternative way, which is quite more practical for recursive estimation, is used in the following to deal with this problem.

The LS-SVM constrained optimization problem is based on the cost function

J (w i,s , e t (s)) = 1 2 ny s=1 d i=0 ω i,s ω i,s + γ 2 ny s=1 N t=d+1 e t (s) 2 (8)
subject to (7), where γ is a regularization parameter and serves as a weight in the cost function to be minimized. In order to apply the Lagrange multipliers method, we formulate the associated Lagrangian as

L(ω i,s , δ s , α t,s , e t,s ) = J (ω i,s , e t,s ) - ny s=1 N t=d+1 α t,s d i=0 ω i,s ϕ(u t-i ) + δ s + e t (s) -y t (s)
where the α t,s stand for the Lagrange multipliers.

The associated Karush-Kuhn-Tucker optimality conditions are given by

∂L ∂ω i,s = 0 → ω i,s = N t=d+1 α t,s ϕ(u t-i ) (9) ∀i = 0, . . . , d ∀s = 1, . . . , n y ∂L ∂δ s = 0 → N t=d+1 α t,s = 0 ∀s = 1, . . . , n y (10) ∂L ∂α t,s = 0, ∀t = d + 1, . . . , N ∀s = 1, . . . , n y → y t (s) = d i=0 ω i,s ϕ(u t-i ) + δ s + e t (s) (11) ∂L ∂e t (s) = 0 → e t (s) = γ -1 α t,s , (12)
∀t = d + 1, . . . , N ∀s = 1, . . . , n y

Taking into account ( 9) and ( 12) in ( 11) gives

y t (s) = d i=0   N q=d+1 α q,s ϕ(u q-i ) ϕ(u t-i )   (13) + δ s + γ -1 α t,s
The kernel function is defined from R nu ×R nu to R as K(u q , u t ) = ϕ(u q ) ϕ(u t ) and selected here to be a gaussian Radial Basis Function K(u q , u t ) = exp(-u q -u t 2 2 /σ 2 ), with σ a constant. Equations ( 10), ( 13) can be summarized in the following linear equation with as unknowns the dual parameters of the optimization problem

0 1 M 1 M Ω + γ -1 I M δ s ᾱs = 0 Ȳs ∀s = 1, . . . , n y (14)
where M = N -d, I M is the identity matrix of order M and Ω = (Ω ij ) ∈ R M×M , 1 M , ᾱs and Ȳs ∈ R M are defined as

Ω ij = d m=0 K(u i+d-m , u j+d-m ) , 1 M = 1 . . . 1 , ᾱs = α d+1,s . . . α N,s , Ȳs = y d+1 (s) . . . y N (s) .
The solution of equation ( 14) is obtained as

δ s = 1 M T -1 Ȳs /(1 M T -1 1 M ) (15) ᾱs = T -1 Ȳs -1 M δ s , s = 1, . . . , n y
where T = Ω + γ -1 I M . Note that matrix T is the same for all the outputs elements s. It depends only on the input u through the kernel function K. Obviously, the properties of that matrix (e.g., its conditioning) are directly related to those of the input.

Once the dual parameters α t,s and δ s are known, the searched nonlinear components from (6) can be expressed as:

ω i,s ϕ(•) = N t=d+1 α t,s K(u t-i , •) = ᾱ s ki (•), (16) s = 1, . . . , n y with ki (•) = K(u d+1-i , •) . . . K(u N -i , •) .
It should be noticed that, depending on the selected kernel, the LS-SVM solution to the nonlinear regression problem is very similar to that obtained using the standard techniques of expanding the nonlinearity on polynomials or RBF bases. The solution ( 16) is known to be nonsparse, that is, some of the terms α t,s K(u t-i , •) have a too small contribution in the general error reduction and can therefore be removed without a great impact. Some pruning techniques have been proposed to render it sparse (see, e.g., (de Kruif and de Vries, 2003)) but our presentation does not address this problem.

It is known that the over-parameterization technique may result in an estimate for g i,s which does not satisfy any longer the relation ( 6). That means that the estimates for the components g i,s may not be collinear as they are expected to be. More precisely, the estimate may be corrupted by an unknown additive component ε i,s (•) which obeys to

d i=0 ε i,s (•) = 0 1 . Hence, we have ĝi,s = ω i,s ϕ(•) + δi,s = H i (s)f (•) + ε i,s (•).
Summing the above equation along the subscript i results in

d i=0 ω i,s ϕ(•) + δs = µ s f (•)
where µ s = i H i (s) is supposed to be nonzero. Note that for any nonzero scalar λ, H i (s)f (•) = (H i (s)λ) 1 λ f (•) . In other words f (•) can only be estimated up to a scalar factor just as the matrices (A, B, C, D) can be determined only up to a similar transformation. In the light of this remark, the static nonlinearity f can arbitrarily set to be f (

•) ← µ s f (•) f (•) = d i=0 ω i,s ϕ(•) + δs , s = 1, . . . , n y (17)
Since the nonlinearity f (•) does not depend on a particular output, it is not worth repeating this estimation for every output as that may increase the computational cost. The nonlinear function f (•) can be extracted from the measurements of the input and the ones of a single arbitrary output s or a combination of all the outputs. ¿From now on, given that the nonlinearity is known, the system matrices can be computed using any subspace identification method. The new linear system to be considered is that of input f (u) and of output the whole output vector y (instead of one component as above). We propose to

1 Indeed, ε(ζt) = ε 0 (ζ 1 t ) + • • • + ε d (ζ d+1 t
) = 0, {ζt} being an i.i.d. vector sequence. We assume then that ε 0 (u)+• • •+ ε d (u) = 0 for any scalar u estimate first the Markov parameters H i from the FIR model (4) using least squares regression and then compute the system matrices. The regression equation is:

y t = H 0 • • • H d    f (u t ) . . . f (u t-d )    + e t , (18)
t = d + 1, • • • , N
When this step is completed, all the terms H i i.e., D, CB, CAB, . . . , CA d-1 B and f are known. From these parameters the system matrices can be directly computed using an SVD but we are interested here in developing a recursive scheme which avoids this SVD process (which is known to be computationally very demanding).

RECURSIVE ALGORITHM

In the following, a recursive version of the identification algorithm presented in the previous section will be derived. While most recursive identification algorithms are developed by using a forgetting factor, in the present case we opt for a sliding window approach. For, assume that at a given instant N observations {u t , y t } k+N -1 t=k are available. When a new pair of data is acquired, the window is slid to be {u t , y t } k+N t=k+1 , that is, the oldest data is removed whenever a new sample becomes available. Then, matrix T changes in such a way that its dimensions remain fixed. A similar approach has been used in [START_REF] Liu | Online LS-SVM for function estimation and classification[END_REF] to implement an online LS-SVM based classifier.

We assume that the nonlinearity f (•) may undergo some changes as well in its parameters as in its form during time. It is then necessary to update its estimate. At times k and k + 1 the matrix T in ( 15) is formally given by

T (k) = η k ϕ k ϕ k ∆ k , T (k + 1) = ∆ k θ k+1 θ k+1 υ k+1
where the underlined elements indicate column vectors. Making use of the matrix block inversion identity we get:

T (k) -1 = g -1 k 1 -ϕ k ∆ -1 k -∆ -1 k ϕ k g k I + ∆ -1 k ϕ k ϕ k ∆ -1 k (19) with g k = η k -ϕ k ∆ -1 k ϕ k ∈ R. Assume now that T (k) -1 is known. Then, if we partition it as T (k) -1 = a k q k q k Q k ∆ -1
k can be deduced with respect to (19) as follows

∆ -1 k = Q k - 1 a k q k q k .
After that T (k + 1) -1 may be computed as:

T (k + 1) -1 = h -1 k+1 (h k+1 I + ∆ -1 k θ k+1 θ k+1 )∆ -1 k -∆ -1 k θ k+1 -θ k+1 ∆ -1 k 1 (20) with h k+1 = υ k+1 -θ k+1 ∆ -1
k θ k+1 ∈ R The dual parameters ᾱs and δ s follow immediately from T (k + 1) -1 and then, the nonlinear function f (•). Using the obtained estimate of the nonlinearity, the parameters H i are computed in a recursive least squares fashion from equation ( 18).

In order to recursively update the estimates of the system matrices, results from the field of recursive subspace identification (see, e.g., [START_REF] Lovera | Recursive subspace identification of linear and non-linear Wiener state-space models[END_REF][START_REF] Mercère | A new recursive method for subspace identification of noisy systems: EIVPM[END_REF] and the references therein) will be exploited. In particular, the socalled Propagator algorithm for system identification will be adopted [START_REF] Mercère | A new recursive method for subspace identification of noisy systems: EIVPM[END_REF], which is briefly summarized in the following. Consider the data vectors

y n(t) = y(t) • • • y(t + n -1) u n(t) = u(t) • • • u(t + n -1)
,

where n is chosen such that n x ≤ n ≤ d. From

(1), it is possible to derive this following relation

z n(t) = Γ nx(t) + n n(t) (21)
where z n = y n(t)-H nf (u n(t)), Γ n is the extended observability matrix, H n is the block Toeplitz matrix of the impulse responses from u to y and n n a sum of all the noise terms.

H n =      D 0 • • • 0 CB D • • • 0 . . . . . . . . . . . . CA n-2 B • • • CB D      Γ n = (C) (CA) • • • (CA n-1 )
The so-called observation vector is explicitly computed as z n(t) = y n(t) -Ĥn f (u n(t)) thanks to the estimates of H n and f (•) obtained from the previous subsection. Under the assumption that the system in (1) is observable, Γ n has at least n x linearly independent rows which can be gathered in a submatrix Γ 1 . Then, based on rows permutation, Γ n can be partitioned as

Γ n = Γ 1 P nΓ 1 = I nx P n Γ 1 ( 22
)
where the matrix P n is the so-called propagator [START_REF] Mercère | A new recursive method for subspace identification of noisy systems: EIVPM[END_REF]. Thus, since Γ 1 is nonsingular, a basis of the observability matrix can be estimated by estimating the propagator. To that purpose the following criterion is considered

J (P n) = R z2η -P nR z1η 2 F ( 23
)
where η is an instrumental variable supposed to be asymptotically uncorrelated with the noise but sufficiently correlated with the state vector,

R z2η = E(z 2 η ), R z1η = E(z 1 η ) are cross- correlation matrices. z 1 z 2
is a partition of z n consistent with the one for Γ n. The propagator is then obtained as P n(t) = R z2η (t)R -1 z1η (t) (see [START_REF] Mercère | A new recursive method for subspace identification of noisy systems: EIVPM[END_REF] for details in the recursive estimation of P n(t)) and the system matrices are finally recovered from

C = Γ n(1 : l) A = Γ n(1 : l(n -1), :) † Γ n(l + 1 : ln, :) D = H 0 , B = Γ n(1 : l(n -1), :) †    H 1 . . . H n-1   
where † denotes the Moore-Penrose pseudo-inverse.

Finally, the comlete algorithm can be summarised as follows. Whenever a new pair (u t , y t ) of data is measured:

(1) Compute the nonlinear function by computing the parameters ᾱ and δ. That is realized in a recursive way by using ( 20) and (15) (2) Update the Markov parameters from (18) using recursive least squares (3) Update the extended observability matrix Γ n by adapting recursively the propagator (4) Compute the system matrices A, B, C, D

SIMULATION RESULTS

To evaluate the proposed algorithm the following SISO example is considered:

x t+1 = -0.15 0.325 1.0 0 x t + 2 0 f (u t ) (24)
y t = 0.1750 1.1625

x t + f (u t ) + υ t
The static nonlinearity is chosen as f (u) = sinc(u)u 2 and the excitation as a white gaussian noise sequence u(t) ∼ N (0, 1). The output noise υ t also white and gaussian is chosen such that the signal to noise ratio is 30 dB. The algorithm parameters are chosen as σ = 1, d = 10, γ -1 = 0.003, n = 6. The simulation is driven with N = 300 over a time interval of 3000 points. Figure 1 illustrates the products between the nonlinearity and the two first Markov parameters on the interval [-5, 5]. The actual functions and their respective estimates are almost undistinguishable, at least on the density interval of the input sequence.

Figure 2 shows the poles of the system estimated recursively. Finally, Figure 3 shows the estimated poles of a time varying system derived from (24) by carrying out a slow and regular variation in the dynamics matrix as follows A(t) = 1 -10 -2 √ t -0.15 0.325 1.0 0 (25)

The results presented show a quite satisfactory tracking of the system parameters.

CONCLUDING REMARKS

A recursive algorithm combining LS-SVM regression and subspace identification techniques has been presented for the identification of SIMO Hammerstein systems. The whole algorithm works in a recursive fashion and is designed for slowly time-varying systems identification. Numerical simulations show quite good results. Future work will focus on the recursive identification of more general nonlinear state space models.

Fig. 1 .

 1 Fig. 1. (a): true (solid line) and estimated (dotted line) Markov parameter Df(). (b): true (solid line) and estimated (dotted line) Markov parameter CBf().

Fig. 2 .

 2 Fig. 2. Evolution of the eigenvalues of the estimated model (dotted) and of the actual system (24) (solid).

Fig. 3 .

 3 Fig. 3. Recursive estimates of the eigenvalues for the slowly time-varying system (25).